Select Git revision
Forked from
agnumpde / dune-elasticity
206 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nonlinelast.cc 12.61 KiB
#include <config.h>
#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/istl/io.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>
#include <dune/solvers/solvers/trustregionsolver.hh>
#include <dune/elasticity/common/nonlinearelasticityproblem.hh>
#define LAURSEN
#include <dune/elasticity/materials/neohookeanmaterial.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> > MatrixType;
typedef BlockVector<FieldVector<double, dim> > VectorType;
// parse data file
ParameterTree parameterSet;
if (argc==2)
ParameterTreeParser::readINITree(argv[1], parameterSet);
else
ParameterTreeParser::readINITree("nonlinelast.parset", parameterSet);
// read solver settings
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");
// /////////////////////////////
// Generate the grid
// /////////////////////////////
typedef UGGrid<dim> GridType;
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);
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);
mmgStep.setObstacleRestrictor(MandelObstacleRestrictor<VectorType>());
mmgStep.setVerbosity(parameterSet.get<NumProc::VerbosityMode>("verbosity"));
// Create a base solver
#ifdef HAVE_IPOPT
QuadraticIPOptSolver<MatrixType,VectorType> baseSolver(parameterSet.get<field_type>("baseTolerance"),
100, NumProc::QUIET);
#endif
mmgStep.setBaseSolver(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 NeoHookeanMaterial<P1Basis> MaterialType;
MaterialType material(p1Basis, parameterSet.get<double>("E"),
parameterSet.get<double>("nu"));
// ///////////////////////////////////////////////////
// Do a homotopy of the Dirichlet boundary data
// ///////////////////////////////////////////////////
field_type loadFactor = 0;
field_type loadIncrement = parameterSet.get<double>("loadIncrement");
for (int level=minLevel; level<=maxLevel; level++) {
VectorType extForces(grid->size(dim));
extForces = 0;
// Dirichlet values on the fine grid
Dune::BitSetVector<dim> dirichletNodes;
VectorType dirichletValues;
DirichletBCAssembler<GridType>::assembleDirichletBC(*grid,
coarseDirichletNodes,coarseDirichletValues,
dirichletNodes, dirichletValues);
mmgStep.setIgnore(dirichletNodes);
//////////////////////////////////
// Create the transfer operator
//////////////////////////////////
std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<VectorType> > > mgTransfers(grid->maxLevel());
for (size_t i=0; i<mgTransfers.size(); i++) {
mgTransfers[i] = std::make_shared<TruncatedCompressedMGTransfer<VectorType> >();
mgTransfers[i]->setup(*grid,i,i+1);
}
mmgStep.setTransferOperators(mgTransfers);
// //////////////////////////////////////////////////////////
// Create obstacles
// //////////////////////////////////////////////////////////
mmgStep.setHasObstacles(BitSetVector<dim>(dirichletNodes.size(), true));
LaplaceAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement,MatrixType::block_type> laplace;
OperatorAssembler<P1Basis, P1Basis> operatorAssembler(p1Basis, p1Basis);
MatrixType laplaceMatrix;
operatorAssembler.assemble(laplace, laplaceMatrix);
h1SemiNorm.setMatrix(&laplaceMatrix);
//////////////////////////////////////////////////////////////
// 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);
typedef BasisGridFunction<P1Basis,VectorType> BasisGridFunction; VectorType deformedX;
do {
deformedX = x;
if (loadFactor < 1) {
// ////////////////////////////////////////////////////
// Compute new feasible load increment
// ////////////////////////////////////////////////////
loadIncrement *= 2; // double it once, cause we're going to half it at least once, too!
bool isNan;
do {
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;
auto 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;
}
// 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 NeoHookeanMaterial<P2Basis> MaterialTypeP2;
MaterialTypeP2 p2Material(p2Basis, parameterSet.get<field_type>("E"),
parameterSet.get<field_type>("nu"));
// 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);
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;
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(),true);
amiramesh.write(resultPath +"nonlineast.result",1);
} catch (Exception e) {
std::cout << e << std::endl;
}