Skip to content
Snippets Groups Projects
Commit b698c106 authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Rewrite the solver. It now uses the material classes + trust-region solver and... 

Rewrite the solver. It now uses the material classes + trust-region solver and applies some adaptive refinement
parent c5dbb19f
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 269 additions and 207 deletions
nonlinelast.cc
+ 269
− 207
View file @ b698c106
#include <config.h>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/amirameshwriter.hh>
#include <dune/grid/io/file/amirameshreader.hh>
#include <dune/elasticity/assemblers/ogdenassembler.hh>
#include <dune/istl/io.hh>
#include <dune/fufem/boundarypatchprolongator.hh>
#include <dune/fufem/sampleonbitfield.hh>
#include <dune/fufem/readbitfield.hh>
#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/functiontools/gridfunctionadaptor.hh>
#include <dune/fufem/estimators/hierarchicalestimator.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/estimators/fractionalmarking.hh>
#include <dune/fufem/estimators/refinementindicator.hh>
#include <dune/fufem/utilities/dirichletbcassembler.hh>
#ifdef HAVE_IPOPT
#include <dune/solvers/solvers/quadraticipopt.hh>
#endif
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
// Choose a solver
#define IPOPT
//#define GAUSS_SEIDEL
//#define MULTIGRID
//#define IPOPT_BASE
#include <dune/elasticity/common/nonlinearelasticityproblem.hh>
#include <dune/elasticity/common/trustregionsolver.hh>
#define LAURSEN
#include <dune/elasticity/materials/neohookeanmaterial.hh>
#include <dune/elasticity/materials/neohookematerial.hh>
// The grid dimension
const int dim = 3;
typedef double field_type;
using namespace Dune;
int main (int argc, char *argv[]) try
{
// Some types that I need
typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
typedef BlockVector<FieldVector<double, dim> > VectorType;
// parse data file
......@@ -48,248 +55,303 @@ int main (int argc, char *argv[]) try
ParameterTreeParser::readINITree("nonlinelast.parset", parameterSet);
// read solver settings
const int numLevels = parameterSet.get<int>("numLevels");
const int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
const int numNewtonSteps = parameterSet.get<int>("numNewtonSteps");
const int numIt = parameterSet.get("numIt", int(0));
const int nu1 = parameterSet.get("nu1", int(0));
const int nu2 = parameterSet.get("nu2", int(0));
const int mu = parameterSet.get("mu", int(0));
const int baseIt = parameterSet.get("baseIt", int(0));
const double tolerance = parameterSet.get("tolerance", double(0));
const double baseTolerance = parameterSet.get("baseTolerance", double(0));
const double scaling = parameterSet.get<double>("scaling");
const int minLevel = parameterSet.get<int>("minLevel");
const int maxLevel = parameterSet.get<int>("maxLevel");
const bool useH1SemiNorm = parameterSet.get<bool>("useH1SemiNorm");
// read problem settings
std::string path = parameterSet.get<std::string>("path");
std::string resultPath = parameterSet.get<std::string>("resultPath");
std::string gridFile = parameterSet.get<std::string>("gridFile");
std::string dnFile = parameterSet.get<std::string>("dnFile");
std::string dvFile = parameterSet.get<std::string>("dvFile");
// /////////////////////////////
// Generate the grid
// /////////////////////////////
typedef UGGrid<dim> GridType;
typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;
typedef BoundaryPatch<GridType::LeafGridView> LeafBoundaryPatch;
GridType* grid = new GridType;
if (parameterSet.hasKey("parFile")) {
std::string parFile = parameterSet.get<std::string>("parFile");
grid = AmiraMeshReader<GridType>::read(path + gridFile, PSurfaceBoundary<dim-1>::read(path + parFile));
} else
AmiraMeshReader<GridType>::read(*grid, path + gridFile);
// Read coarse Dirichlet boundary values
BitSetVector<dim> coarseDirichletNodes;
VectorType coarseDirichletValues;
DirichletBCAssembler<GridType>::assembleDirichletBC(*grid, parameterSet,
coarseDirichletNodes, coarseDirichletValues,
path, "Amira");
const double scaling = parameterSet.get<double>("scaling");
coarseDirichletValues *= scaling;
GridType grid;
Dune::AmiraMeshReader<GridType>::read(grid, gridFile);
// //////////////////////
// Uniform refine grid
// /////////////////////
// Read Dirichlet boundary values
std::vector<LevelBoundaryPatch> dirichletBoundary(numLevels);
dirichletBoundary[0].setup(grid.levelGridView(0));
readBoundaryPatch<GridType>(dirichletBoundary[0], dnFile);
grid->setRefinementType(GridType::COPY);
for (int i=0; i<minLevel; i++)
grid->globalRefine(1);
std::vector<VectorType> dirichletValues(numLevels);
dirichletValues[0].resize(grid.size(0, dim));
AmiraMeshReader<GridType>::readFunction(dirichletValues[0], dvFile);
// Initial solution
VectorType x(grid->size(dim));
x = 0;
dirichletValues[0] *= scaling;
// /////////////////////////////////////
//setup the monotone multigrid step
// /////////////////////////////////////
for (int i=0; i<numLevels-1; i++)
grid.globalRefine(1);
BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary);
MonotoneMGStep<MatrixType,VectorType> mmgStep;
mmgStep.setMGType(parameterSet.get("mu",1),
parameterSet.get("preSmoothingSteps",3),
parameterSet.get("postSmoothingSteps",3));
std::vector<BitSetVector<dim> > dirichletNodes(numLevels);
for (int i=0; i<numLevels; i++) {
dirichletNodes[i].resize(grid.size(i, dim));
for (int j=0; j<grid.size(i, dim); j++) {
TrustRegionGSStep<MatrixType, VectorType> presmoother,postsmoother;
mmgStep.setSmoother(&presmoother, &postsmoother);
for (int k=0; k<dim; k++)
dirichletNodes[i][j][k] = dirichletBoundary[i].containsVertex(j);
}
MandelObstacleRestrictor<VectorType> obsRestrictor;
mmgStep.obstacleRestrictor_= &obsRestrictor;
mmgStep.verbosity_ = parameterSet.get<NumProc::VerbosityMode>("verbosity");
}
// Create a base solver
#ifdef HAVE_IPOPT
// hack Dirichlet data
dirichletNodes[numLevels-1][2*61 ] = false;
dirichletNodes[numLevels-1][2*61+1] = false;
dirichletNodes[numLevels-1][2*62 ] = false;
dirichletNodes[numLevels-1][2*62+1] = false;
dirichletNodes[numLevels-1][2*56 ] = false;
dirichletNodes[numLevels-1][2*56+1] = false;
dirichletNodes[numLevels-1][2*58 ] = false;
dirichletNodes[numLevels-1][2*58+1] = false;
dirichletNodes[numLevels-1][2*79 ] = false;
dirichletNodes[numLevels-1][2*79+1] = false;
dirichletNodes[numLevels-1][2*74 ] = false;
dirichletNodes[numLevels-1][2*74+1] = false;
dirichletNodes[numLevels-1][2*76 ] = false;
dirichletNodes[numLevels-1][2*76+1] = false;
dirichletNodes[numLevels-1][2*73 ] = false;
dirichletNodes[numLevels-1][2*73+1] = false;
dirichletNodes[numLevels-1][2*4 ] = false;
dirichletNodes[numLevels-1][2*4+1] = false;
dirichletNodes[numLevels-1][2*1 ] = false;
dirichletNodes[numLevels-1][2*1+1] = false;
dirichletNodes[numLevels-1][2*11 ] = false;
dirichletNodes[numLevels-1][2*11+1] = false;
dirichletNodes[numLevels-1][2*9 ] = false;
dirichletNodes[numLevels-1][2*9+1] = false;
dirichletNodes[numLevels-1][2*27 ] = false;
dirichletNodes[numLevels-1][2*27+1] = false;
dirichletNodes[numLevels-1][2*25 ] = false;
dirichletNodes[numLevels-1][2*25+1] = false;
dirichletNodes[numLevels-1][2*33 ] = false;
dirichletNodes[numLevels-1][2*33+1] = false;
sampleOnBitField(grid, dirichletValues, dirichletNodes);
dirichletValues[numLevels-1][0][1] = 0.2;
dirichletValues[numLevels-1][31][1] = 0.2;
int maxlevel = grid.maxLevel();
// Determine Dirichlet dofs
VectorType rhs(grid.size(grid.maxLevel(), dim));
VectorType x(grid.size(grid.maxLevel(), dim));
VectorType corr(grid.size(grid.maxLevel(), dim));
OperatorType problemMatrix;
QuadraticIPOptSolver<MatrixType,VectorType> baseSolver;
baseSolver.verbosity_ = NumProc::QUIET;
baseSolver.tolerance_ = parameterSet.get<field_type>("baseTolerance");
#endif
mmgStep.basesolver_ = &baseSolver;
// Initial solution
x = 0;
corr = 0;
rhs = 0;
EnergyNorm<MatrixType, VectorType> h1SemiNorm;
EnergyNorm<MatrixType, VectorType > energyNorm(mmgStep);
// Create a solver
#if defined IPOPT
::LoopSolver<VectorType> solver(&mmgStep,
parameterSet.get<int>("mgIterations"),
parameterSet.get<field_type>("mgTolerance"),
(useH1SemiNorm) ? &h1SemiNorm : &energyNorm,
Solver::FULL);
QuadraticIPOptSolver<OperatorType,VectorType> solver;
solver.verbosity_ = NumProc::QUIET;
solver.ignoreNodes_ = &dirichletNodes[maxlevel];
// /////////////////////////////////
// Setup the nonlinear material
// ////////////////////////////////
//
typedef P1NodalBasis<GridType::LeafGridView> P1Basis;
P1Basis p1Basis(grid->leafGridView());
typedef NeoHookeMaterial<P1Basis> MaterialType;
MaterialType material(p1Basis,
parameterSet.get<double>("E"),
parameterSet.get<double>("nu"),
parameterSet.get<bool>("vectorMode"));
#elif defined GAUSS_SEIDEL
L2Norm<ContactVector<dim> > l2Norm;
// ///////////////////////////////////////////////////
// Do a homotopy of the Dirichlet boundary data
// ///////////////////////////////////////////////////
typedef ProjectedBlockGSStep<OperatorType, VectorType> SmootherType;
SmootherType projectedBlockGSStep(*bilinearForm, x, rhs);
projectedBlockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];
field_type loadFactor = 0;
field_type loadIncrement = parameterSet.get<double>("loadIncrement");
::LoopSolver<OperatorType, ContactVector<dim> > solver;
solver.iterationStep = &projectedBlockGSStep;
solver.numIt = numIt;
solver.verbosity_ = Solver::FULL;
solver.errorNorm_ = &l2Norm;
solver.tolerance_ = tolerance;
int counter =0;
#elif defined MULTIGRID
L2Norm<ContactVector<dim> > l2Norm;
for (int level=minLevel; level<=maxLevel; level++) {
// First create a base solver
#ifdef IPOPT_BASE
LinearIPOptSolver baseSolver;
VectorType extForces(grid->size(dim));
extForces = 0;
#else // Gauss-Seidel is the base solver
// Dirichlet values on the fine grid
Dune::BitSetVector<dim> dirichletNodes;
VectorType dirichletValues;
ProjectedBlockGSStep<OperatorType, ContactVector<dim> > baseSolverStep;
DirichletBCAssembler<GridType>::assembleDirichletBC(*grid,
coarseDirichletNodes,coarseDirichletValues,
dirichletNodes, dirichletValues);
::LoopSolver<OperatorType, ContactVector<dim> > baseSolver;
baseSolver.iterationStep = &baseSolverStep;
baseSolver.numIt = baseIt;
baseSolver.verbosity_ = Solver::QUIET;
baseSolver.errorNorm_ = &l2Norm;
baseSolver.tolerance_ = baseTolerance;
#endif
mmgStep.ignoreNodes_ = &dirichletNodes;
// Make pre and postsmoothers
ProjectedBlockGSStep<OperatorType, ContactVector<dim> > presmoother;
ProjectedBlockGSStep<OperatorType, ContactVector<dim> > postsmoother;
//////////////////////////////////
// Create the transfer operator
//////////////////////////////////
std::vector<TruncatedCompressedMGTransfer<VectorType>* > mgTransfers(grid->maxLevel());
for (size_t i=0; i<mgTransfers.size(); i++) {
ContactMMGStep<OperatorType, ContactVector<dim>, FuncSpaceType > contactMMGStep(*bilinearForm, x, rhs, numLevels);
mgTransfers[i] = new TruncatedCompressedMGTransfer<VectorType>;
mgTransfers[i]->setup(*grid,i,i+1);
}
contactMMGStep.setMGType(1, nu1, nu2);
contactMMGStep.dirichletNodes_ = &dirichletNodes;
contactMMGStep.basesolver_ = &baseSolver;
contactMMGStep.presmoother_ = &presmoother;
contactMMGStep.postsmoother_ = &postsmoother;
contactMMGStep.funcSpace_ = funcSpace;
mmgStep.setTransferOperators(mgTransfers);
::LoopSolver<OperatorType, ContactVector<dim> > solver;
solver.iterationStep = &contactMMGStep;
solver.numIt = numIt;
solver.verbosity_ = Solver::FULL;
solver.errorNorm_ = &l2Norm;
solver.tolerance_ = tolerance;
// //////////////////////////////////////////////////////////
// Create obstacles
// //////////////////////////////////////////////////////////
#else
#warning You have to specify a solver!
#endif
typedef P1NodalBasis<GridType::LeafGridView,double> P1Basis;
typedef OgdenAssembler<P1Basis,P1Basis> Assembler;
BitSetVector<dim> hasObstacle(dirichletNodes.size(), true);
mmgStep.hasObstacle_ = &hasObstacle;
P1Basis p1Basis(grid.leafGridView());
Assembler ogdenAssembler(p1Basis,p1Basis);
LaplaceAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement,MatrixType::block_type> laplace;
OperatorAssembler<P1Basis, P1Basis> operatorAssembler(p1Basis, p1Basis);
for (int i=0; i<numHomotopySteps; i++) {
MatrixType laplaceMatrix;
operatorAssembler.assemble(laplace, laplaceMatrix);
h1SemiNorm.setMatrix(&laplaceMatrix);
// scale Dirichlet values
VectorType scaledDirichlet = dirichletValues[maxlevel];
scaledDirichlet *= (double(i)+1)/numHomotopySteps;
// Set new Dirichlet values in solution
for (int j=0; j<x.size(); j++)
for (int k=0; k<dim; k++)
if (dirichletNodes[maxlevel][j][k])
x[j][k] = scaledDirichlet[j][k];
//////////////////////////////////////////////////////////////
// Create the nonlinear elasticity problem
/////////////////////////////////////////////////////////////
// Make static contact problem
typedef NonlinearElasticityProblem<VectorType,MatrixType, P1Basis> ElasticityProblem;
ElasticityProblem elasticProblem(material, extForces);
// Create trust-region step
TrustRegionSolver<ElasticityProblem, VectorType, MatrixType> trustRegionSolver;
//TODO Allow a different norm here
const ParameterTree& trConfig = parameterSet.sub("trConfig");
trustRegionSolver.setup(trConfig, h1SemiNorm, elasticProblem, solver);
for (int j=0; j<numNewtonSteps; j++) {
std::cout << "iteration: " << j << std::endl;
typedef BasisGridFunction<P1Basis,VectorType> BasisGridFunction;
VectorType deformedX;
// Assemble the Jacobi matrix at the current solution
rhs = 0;
ogdenAssembler.assembleProblem(problemMatrix,x, rhs);
do {
corr = 0;
deformedX = x;
// Set new Dirichlet values in the right hand side
for (int k=0; k<rhs.size(); k++)
for (int l=0; l<dim; l++)
if (dirichletNodes[maxlevel][k][l])
rhs[k][l] = 0;
if (loadFactor < 1) {
solver.setProblem(problemMatrix, corr, rhs);
// ////////////////////////////////////////////////////
// Compute new feasible load increment
// ////////////////////////////////////////////////////
solver.preprocess();
#ifdef MULTIGRID
contactMMGStep.preprocess();
#endif
// Solve correction problem
solver.solve();
loadIncrement *= 2; // double it once, cause we're going to half it at least once, too!
bool isNan;
do {
// Add damped correction to the current solution
x += corr;
loadIncrement /= 2;
std::cout << "Correction: " << corr.infinity_norm() << std::endl;
// Set new Dirichlet values in solution
for (size_t k=0; k<dirichletNodes.size(); k++)
for (int j=0; j<dim; j++)
if (dirichletNodes[k][j])
deformedX[k][j] = dirichletValues[k][j] * (loadFactor + loadIncrement);
}
isNan = false;
std::shared_ptr<BasisGridFunction> disp = std::make_shared<BasisGridFunction>(p1Basis,deformedX);
field_type eng = material.energy(disp);
std::cout<<"Energy "<<eng<<std::endl;
isNan |= std::isnan(eng);
} while (isNan);
loadFactor += loadIncrement;
std::cout << "####################################################" << std::endl;
std::cout << "New load factor: " << loadFactor
<< " new load increment: " << loadIncrement << std::endl;
std::cout << "####################################################" << std::endl;
}
#ifdef MULTIGRID
x = contactMMGStep.getSol();
#endif
// add Dirichlet values
trustRegionSolver.setInitialIterate(deformedX);
trustRegionSolver.solve();
x = trustRegionSolver.getSol();
} while (loadFactor < 1);
// /////////////////////////////////////////////////////////////////////
// Refinement Loop
// /////////////////////////////////////////////////////////////////////
if (level==maxLevel)
break;
// ////////////////////////////////////////////////////////////////////////////
// Refine locally and transfer the current solution to the new leaf level
// ////////////////////////////////////////////////////////////////////////////
std::cout<<"Estimating error on refinement level "<<level<<std::endl;
// Create the materials
typedef P2NodalBasis<GridType::LeafGridView> P2Basis;
P2Basis p2Basis(grid->leafGridView());
typedef NeoHookeMaterial<P2Basis> MaterialTypeP2;
MaterialTypeP2 p2Material(p2Basis, parameterSet.get<field_type>("E"),
parameterSet.get<field_type>("nu"),
parameterSet.get<bool>("vectorMode"));
// P2 Forces
VectorType p2ExtForces(p2Basis.size());
p2ExtForces = 0;
typedef NonlinearElasticityProblem<VectorType,MatrixType, P2Basis> ElasticityProblemP2;
ElasticityProblemP2 elasticProblemP2(p2Material, p2ExtForces);
RefinementIndicator<GridType> refinementIndicator(*grid);
BitSetVector<1> scalarDirichletField;
DirichletBCAssembler<GridType>::inflateBitField(dirichletNodes, scalarDirichletField);
BoundaryPatch<GridType::LeafGridView> dirichletBoundary(grid->leafGridView(),scalarDirichletField);
//contactAssembler.postprocess(x);
HierarchicalEstimator<P2Basis,dim> estimator(*grid);
estimator.estimate(x, dirichletBoundary, refinementIndicator,elasticProblemP2);
//std::cout << x << std::endl;
// Output result
// ////////////////////////////////////////////////////
// Refine grids
// ////////////////////////////////////////////////////
std::cout<<"Mark elements"<<std::endl;
FractionalMarkingStrategy<GridType>::mark(refinementIndicator, *grid,
parameterSet.get<field_type>("refinementFraction"));
GridFunctionAdaptor<P1Basis> adaptor(p1Basis,true,true);
grid->preAdapt();
grid->adapt();
grid->postAdapt();
p1Basis.update();
adaptor.adapt(x);
std::cout << "########################################################" << std::endl;
std::cout << " Grid refined," <<" max level: " << grid->maxLevel()
<< " vertices: " << grid->size(dim)
<< " elements: " << grid->size(0)<<std::endl;
std::cout << "####################################################" << std::endl;
for (size_t j=0; j<mgTransfers.size(); j++)
delete(mgTransfers[j]);
Dune::LeafAmiraMeshWriter<GridType> amiramesh2;
amiramesh2.addLeafGrid(*grid,true);
amiramesh2.addVertexData(x,grid->leafGridView(),true);
std::ostringstream name;
name << "nonlinelastRefined" << std::setw(6) << std::setfill('0') <<level ;
amiramesh2.write(resultPath +name.str(),1);
}
LeafAmiraMeshWriter<GridType> amiramesh;
amiramesh.addLeafGrid(grid,true);
amiramesh.addVertexData(x, grid.leafGridView());
amiramesh.write("resultGrid",true);
amiramesh.addLeafGrid(*grid,true);
amiramesh.addVertexData(x, grid->leafGridView(),true);
amiramesh.write(resultPath +"nonlineast.result",1);
} catch (Exception e) {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment