-
Jonathan Youett authoredJonathan Youett authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nonlinelast.cc 14.53 KiB
#include <config.h>
#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/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>
#include <dune/fufem/utilities/gridconstruction.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>
#include <dune/elasticity/materials/adolcmaterial.hh>
#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
// The grid dimension
const int dim = 3;
typedef double field_type;
using namespace Dune;
int main (int argc, char *argv[]) try
{
Dune::MPIHelper::instance(argc, argv);
// 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");
// /////////////////////////////
// Generate the grid // /////////////////////////////
typedef UGGrid<dim> GridType;
std::unique_ptr<GridType> grid;
const auto& gridConfig = parameterSet.sub(parameterSet.get<std::string>("gridName"));
if (gridConfig.get<bool>("createGrid", false)) {
grid = GridConstruction<GridType,dim>::createGrid(gridConfig);
} else if (gridConfig.hasKey("amira_parFile")) {
#if HAVE_AMIRAMESH
auto gridFile = gridConfig.get<std::string>("gridFile");
auto parFile = gridConfig.get<std::string>("amira_parFile");
grid = AmiraMeshReader<GridType>::read(path + gridFile, PSurfaceBoundary<dim-1>::read(path + parFile));
#endif
} else {
#if HAVE_AMIRAMESH
grid = AmiraMeshReader<GridType>::read(path + gridConfig.get<std::string>("gridFile"));
#endif
}
// Read coarse Dirichlet boundary condition
BitSetVector<dim> coarseDirichletNodes(grid->size(dim),false);
VectorType coarseDirichletValues(grid->size(dim));
coarseDirichletValues = 0;
for (int i = 0; i < gridConfig.get<int>("nDirichletConditions"); ++i) {
const auto& dirConfig = gridConfig.sub(formatString("dirichlet%d",i));
BitSetVector<dim> dirichletNodes;
VectorType dirichletValues;
DirichletBCAssembler<GridType>::assembleDirichletBC(*grid, dirConfig,
dirichletNodes, dirichletValues, path);
auto scaling = dirConfig.get<double>("dvScaling", 1);
// combine with previous ones
coarseDirichletValues.axpy(scaling, dirichletValues);
for (size_t j = 0; j < dirichletNodes.size(); ++j)
for (int k = 0; k < dim; ++k)
if (dirichletNodes[j][k])
coarseDirichletNodes[j][k] = true;
}
// 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());
#if HAVE_ADOLC
MooneyRivlinMaterial<P1Basis> localEnergy(p1Basis,
parameterSet.get<double>("E"),
parameterSet.get<double>("nu"));
using MaterialType = AdolcMaterial<P1Basis>;
MaterialType material(p1Basis, localEnergy, parameterSet.get<bool>("vectorMode"));
#else
using MaterialType = NeoHookeanMaterial<P1Basis>;
MaterialType material(p1Basis,
parameterSet.get<double>("E"),
parameterSet.get<double>("nu"));
#endif
// ///////////////////////////////////////////////////
// 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();
LeafAmiraMeshWriter<GridType> amiramesh2;
amiramesh2.addLeafGrid(*grid,true);
amiramesh2.addVertexData(x, grid->leafGridView(),true);
std::string name = "loadingStep" + std::to_string(loadFactor);
amiramesh2.write(resultPath + name, 1);
} 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());
#if HAVE_ADOLC
typedef MooneyRivlinMaterial<P2Basis> MaterialType2;
MaterialType2 p2localEnergy(p2Basis, parameterSet.get<field_type>("E"),
parameterSet.get<field_type>("nu"));
AdolcMaterial<P2Basis> p2Material(p2Basis, p2localEnergy, false);
#else
using MaterialType2 = NeoHookeanMaterial<P2Basis>;
MaterialType material(p2Basis,
parameterSet.get<double>("E"),
parameterSet.get<double>("nu"));
#endif
// 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;
}