diff --git a/nonlinelast.cc b/nonlinelast.cc
index b2cd719e0224c8a5d88759edcac293b4395a232c..c29a256046b24704a8fd880a5900ac85bb62c0b5 100644
--- a/nonlinelast.cc
+++ b/nonlinelast.cc
@@ -1,43 +1,50 @@
#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,254 +55,309 @@ 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;
+
+
+ // //////////////////////
+ // Uniform refine grid
+ // /////////////////////
+
+ grid->setRefinementType(GridType::COPY);
+ for (int i=0; i<minLevel; i++)
+ grid->globalRefine(1);
+
+ // Initial solution
+ VectorType x(grid->size(dim));
+ x = 0;
+
+ // /////////////////////////////////////
+ //setup the monotone multigrid step
+ // /////////////////////////////////////
+
+
+ MonotoneMGStep<MatrixType,VectorType> mmgStep;
+ mmgStep.setMGType(parameterSet.get("mu",1),
+ parameterSet.get("preSmoothingSteps",3),
+ parameterSet.get("postSmoothingSteps",3));
+
+ TrustRegionGSStep<MatrixType, VectorType> presmoother,postsmoother;
+ mmgStep.setSmoother(&presmoother, &postsmoother);
+
+ MandelObstacleRestrictor<VectorType> obsRestrictor;
+ mmgStep.obstacleRestrictor_= &obsRestrictor;
+ mmgStep.verbosity_ = parameterSet.get<NumProc::VerbosityMode>("verbosity");
+
+ // Create a base solver
+#ifdef HAVE_IPOPT
+
+ QuadraticIPOptSolver<MatrixType,VectorType> baseSolver;
+ baseSolver.verbosity_ = NumProc::QUIET;
+ baseSolver.tolerance_ = parameterSet.get<field_type>("baseTolerance");
+#endif
+ mmgStep.basesolver_ = &baseSolver;
+
+ EnergyNorm<MatrixType, VectorType> h1SemiNorm;
+ EnergyNorm<MatrixType, VectorType > energyNorm(mmgStep);
+
+ ::LoopSolver<VectorType> solver(&mmgStep,
+ parameterSet.get<int>("mgIterations"),
+ parameterSet.get<field_type>("mgTolerance"),
+ (useH1SemiNorm) ? &h1SemiNorm : &energyNorm,
+ Solver::FULL);
+
+
+ // /////////////////////////////////
+ // 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"));
+
+
+ // ///////////////////////////////////////////////////
+ // Do a homotopy of the Dirichlet boundary data
+ // ///////////////////////////////////////////////////
- GridType grid;
+ field_type loadFactor = 0;
+ field_type loadIncrement = parameterSet.get<double>("loadIncrement");
- Dune::AmiraMeshReader<GridType>::read(grid, gridFile);
-
- // Read Dirichlet boundary values
- std::vector<LevelBoundaryPatch> dirichletBoundary(numLevels);
- dirichletBoundary[0].setup(grid.levelGridView(0));
- readBoundaryPatch<GridType>(dirichletBoundary[0], dnFile);
+ int counter =0;
- std::vector<VectorType> dirichletValues(numLevels);
- dirichletValues[0].resize(grid.size(0, dim));
- AmiraMeshReader<GridType>::readFunction(dirichletValues[0], dvFile);
- dirichletValues[0] *= scaling;
+ for (int level=minLevel; level<=maxLevel; level++) {
- for (int i=0; i<numLevels-1; i++)
- grid.globalRefine(1);
- BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary);
+ VectorType extForces(grid->size(dim));
+ extForces = 0;
- 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++) {
+ // Dirichlet values on the fine grid
+ Dune::BitSetVector<dim> dirichletNodes;
+ VectorType dirichletValues;
- for (int k=0; k<dim; k++)
- dirichletNodes[i][j][k] = dirichletBoundary[i].containsVertex(j);
+ DirichletBCAssembler<GridType>::assembleDirichletBC(*grid,
+ coarseDirichletNodes,coarseDirichletValues,
+ dirichletNodes, dirichletValues);
+
+ mmgStep.ignoreNodes_ = &dirichletNodes;
+
+ //////////////////////////////////
+ // Create the transfer operator
+ //////////////////////////////////
+
+ std::vector<TruncatedCompressedMGTransfer<VectorType>* > mgTransfers(grid->maxLevel());
+ for (size_t i=0; i<mgTransfers.size(); i++) {
+
+ mgTransfers[i] = new TruncatedCompressedMGTransfer<VectorType>;
+ mgTransfers[i]->setup(*grid,i,i+1);
}
- }
+ mmgStep.setTransferOperators(mgTransfers);
- // 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;
-
- // Initial solution
- x = 0;
- corr = 0;
- rhs = 0;
+ // //////////////////////////////////////////////////////////
+ // Create obstacles
+ // //////////////////////////////////////////////////////////
- // Create a solver
-#if defined IPOPT
+ BitSetVector<dim> hasObstacle(dirichletNodes.size(), true);
+ mmgStep.hasObstacle_ = &hasObstacle;
- QuadraticIPOptSolver<OperatorType,VectorType> solver;
- solver.verbosity_ = NumProc::QUIET;
+ LaplaceAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement,MatrixType::block_type> laplace;
+ OperatorAssembler<P1Basis, P1Basis> operatorAssembler(p1Basis, p1Basis);
- solver.ignoreNodes_ = &dirichletNodes[maxlevel];
+ MatrixType laplaceMatrix;
+ operatorAssembler.assemble(laplace, laplaceMatrix);
+ h1SemiNorm.setMatrix(&laplaceMatrix);
-#elif defined GAUSS_SEIDEL
- L2Norm<ContactVector<dim> > l2Norm;
+ //////////////////////////////////////////////////////////////
+ // Create the nonlinear elasticity problem
+ /////////////////////////////////////////////////////////////
- typedef ProjectedBlockGSStep<OperatorType, VectorType> SmootherType;
- SmootherType projectedBlockGSStep(*bilinearForm, x, rhs);
- projectedBlockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];
+ // Make static contact problem
+ typedef NonlinearElasticityProblem<VectorType,MatrixType, P1Basis> ElasticityProblem;
+ ElasticityProblem elasticProblem(material, extForces);
- ::LoopSolver<OperatorType, ContactVector<dim> > solver;
- solver.iterationStep = &projectedBlockGSStep;
- solver.numIt = numIt;
- solver.verbosity_ = Solver::FULL;
- solver.errorNorm_ = &l2Norm;
- solver.tolerance_ = tolerance;
-
-#elif defined MULTIGRID
+ // 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);
- L2Norm<ContactVector<dim> > l2Norm;
- // First create a base solver
-#ifdef IPOPT_BASE
+ typedef BasisGridFunction<P1Basis,VectorType> BasisGridFunction;
+ VectorType deformedX;
- LinearIPOptSolver baseSolver;
+ do {
-#else // Gauss-Seidel is the base solver
+ deformedX = x;
- ProjectedBlockGSStep<OperatorType, ContactVector<dim> > baseSolverStep;
+ if (loadFactor < 1) {
- ::LoopSolver<OperatorType, ContactVector<dim> > baseSolver;
- baseSolver.iterationStep = &baseSolverStep;
- baseSolver.numIt = baseIt;
- baseSolver.verbosity_ = Solver::QUIET;
- baseSolver.errorNorm_ = &l2Norm;
- baseSolver.tolerance_ = baseTolerance;
-#endif
+ // ////////////////////////////////////////////////////
+ // Compute new feasible load increment
+ // ////////////////////////////////////////////////////
- // Make pre and postsmoothers
- ProjectedBlockGSStep<OperatorType, ContactVector<dim> > presmoother;
- ProjectedBlockGSStep<OperatorType, ContactVector<dim> > postsmoother;
-
-
- ContactMMGStep<OperatorType, ContactVector<dim>, FuncSpaceType > contactMMGStep(*bilinearForm, x, rhs, numLevels);
-
- contactMMGStep.setMGType(1, nu1, nu2);
- contactMMGStep.dirichletNodes_ = &dirichletNodes;
- contactMMGStep.basesolver_ = &baseSolver;
- contactMMGStep.presmoother_ = &presmoother;
- contactMMGStep.postsmoother_ = &postsmoother;
- contactMMGStep.funcSpace_ = funcSpace;
-
- ::LoopSolver<OperatorType, ContactVector<dim> > solver;
- solver.iterationStep = &contactMMGStep;
- solver.numIt = numIt;
- solver.verbosity_ = Solver::FULL;
- solver.errorNorm_ = &l2Norm;
- solver.tolerance_ = tolerance;
-
-#else
- #warning You have to specify a solver!
-#endif
- typedef P1NodalBasis<GridType::LeafGridView,double> P1Basis;
- typedef OgdenAssembler<P1Basis,P1Basis> Assembler;
-
- P1Basis p1Basis(grid.leafGridView());
- Assembler ogdenAssembler(p1Basis,p1Basis);
-
- for (int i=0; i<numHomotopySteps; i++) {
- // scale Dirichlet values
- VectorType scaledDirichlet = dirichletValues[maxlevel];
- scaledDirichlet *= (double(i)+1)/numHomotopySteps;
+ loadIncrement *= 2; // double it once, cause we're going to half it at least once, too!
+ bool isNan;
+ do {
- // 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];
+ loadIncrement /= 2;
+ // 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;
- for (int j=0; j<numNewtonSteps; j++) {
+ 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);
- std::cout << "iteration: " << j << std::endl;
+ loadFactor += loadIncrement;
- // Assemble the Jacobi matrix at the current solution
- rhs = 0;
- ogdenAssembler.assembleProblem(problemMatrix,x, rhs);
+ std::cout << "####################################################" << std::endl;
+ std::cout << "New load factor: " << loadFactor
+ << " new load increment: " << loadIncrement << std::endl;
+ std::cout << "####################################################" << std::endl;
- corr = 0;
-
- // 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;
+ }
- solver.setProblem(problemMatrix, corr, rhs);
+ // add Dirichlet values
+ trustRegionSolver.setInitialIterate(deformedX);
+ trustRegionSolver.solve();
+ x = trustRegionSolver.getSol();
- solver.preprocess();
-#ifdef MULTIGRID
- contactMMGStep.preprocess();
-#endif
+ } while (loadFactor < 1);
- // Solve correction problem
- solver.solve();
- // Add damped correction to the current solution
- x += corr;
+ // /////////////////////////////////////////////////////////////////////
+ // Refinement Loop
+ // /////////////////////////////////////////////////////////////////////
- std::cout << "Correction: " << corr.infinity_norm() << std::endl;
+ 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);
-#ifdef MULTIGRID
- x = contactMMGStep.getSol();
-#endif
- //contactAssembler.postprocess(x);
+ RefinementIndicator<GridType> refinementIndicator(*grid);
- //std::cout << x << std::endl;
+ BitSetVector<1> scalarDirichletField;
+ DirichletBCAssembler<GridType>::inflateBitField(dirichletNodes, scalarDirichletField);
+ BoundaryPatch<GridType::LeafGridView> dirichletBoundary(grid->leafGridView(),scalarDirichletField);
+
+ HierarchicalEstimator<P2Basis,dim> estimator(*grid);
+ estimator.estimate(x, dirichletBoundary, refinementIndicator,elasticProblemP2);
+
+
+ // ////////////////////////////////////////////////////
+ // 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);
+
+
+ }
- // Output result
- LeafAmiraMeshWriter<GridType> amiramesh;
- amiramesh.addLeafGrid(grid,true);
- amiramesh.addVertexData(x, grid.leafGridView());
- amiramesh.write("resultGrid",true);
+ LeafAmiraMeshWriter<GridType> amiramesh;
+ amiramesh.addLeafGrid(*grid,true);
+ amiramesh.addVertexData(x, grid->leafGridView(),true);
+ amiramesh.write(resultPath +"nonlineast.result",1);
- } catch (Exception e) {
+} catch (Exception e) {
std::cout << e << std::endl;
- }
+}