Forked from
agnumpde / dune-elasticity
405 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nonlinelast.cc 9.46 KiB
#include <config.h>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.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/p1nodalbasis.hh>
#include <dune/solvers/solvers/quadraticipopt.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
// Choose a solver
#define IPOPT
//#define GAUSS_SEIDEL
//#define MULTIGRID
//#define IPOPT_BASE
// The grid dimension
const int dim = 3;
using namespace Dune;
int main (int argc, char *argv[]) try
{
// Some types that I need
typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
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 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");
// read problem settings
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;
GridType grid;
Dune::AmiraMeshReader<GridType>::read(grid, gridFile);
// Read Dirichlet boundary values
std::vector<LevelBoundaryPatch> dirichletBoundary(numLevels);
dirichletBoundary[0].setup(grid.levelView(0));
readBoundaryPatch<GridType>(dirichletBoundary[0], dnFile);
std::vector<VectorType> dirichletValues(numLevels);
dirichletValues[0].resize(grid.size(0, dim));
AmiraMeshReader<GridType>::readFunction(dirichletValues[0], dvFile);
dirichletValues[0] *= scaling;
for (int i=0; i<numLevels-1; i++)
grid.globalRefine(1);
BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary);
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++) {
for (int k=0; k<dim; k++)
dirichletNodes[i][j][k] = dirichletBoundary[i].containsVertex(j);
}
}
// 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 a solver
#if defined IPOPT
QuadraticIPOptSolver<OperatorType,VectorType> solver;
solver.verbosity_ = NumProc::QUIET;
solver.ignoreNodes_ = &dirichletNodes[maxlevel];
#elif defined GAUSS_SEIDEL
L2Norm<ContactVector<dim> > l2Norm;
typedef ProjectedBlockGSStep<OperatorType, VectorType> SmootherType;
SmootherType projectedBlockGSStep(*bilinearForm, x, rhs);
projectedBlockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];
::LoopSolver<OperatorType, ContactVector<dim> > solver;
solver.iterationStep = &projectedBlockGSStep;
solver.numIt = numIt;
solver.verbosity_ = Solver::FULL;
solver.errorNorm_ = &l2Norm;
solver.tolerance_ = tolerance;
#elif defined MULTIGRID
L2Norm<ContactVector<dim> > l2Norm;
// First create a base solver
#ifdef IPOPT_BASE
LinearIPOptSolver baseSolver;
#else // Gauss-Seidel is the base solver
ProjectedBlockGSStep<OperatorType, ContactVector<dim> > baseSolverStep;
::LoopSolver<OperatorType, ContactVector<dim> > baseSolver;
baseSolver.iterationStep = &baseSolverStep;
baseSolver.numIt = baseIt;
baseSolver.verbosity_ = Solver::QUIET;
baseSolver.errorNorm_ = &l2Norm;
baseSolver.tolerance_ = baseTolerance;
#endif
// 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.leafView());
Assembler ogdenAssembler(p1Basis,p1Basis);
for (int i=0; i<numHomotopySteps; i++) {
// 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];
for (int j=0; j<numNewtonSteps; j++) {
std::cout << "iteration: " << j << std::endl;
// Assemble the Jacobi matrix at the current solution
rhs = 0;
ogdenAssembler.assembleProblem(problemMatrix,x, rhs);
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);
solver.preprocess();
#ifdef MULTIGRID
contactMMGStep.preprocess();
#endif
// Solve correction problem
solver.solve();
// Add damped correction to the current solution
x += corr;
std::cout << "Correction: " << corr.infinity_norm() << std::endl;
}
}
#ifdef MULTIGRID
x = contactMMGStep.getSol();
#endif
//contactAssembler.postprocess(x);
//std::cout << x << std::endl;
// Output result
LeafAmiraMeshWriter<GridType> amiramesh;
amiramesh.addLeafGrid(grid,true);
amiramesh.addVertexData(x, grid.leafView());
amiramesh.write("resultGrid",true);
} catch (Exception e) {
std::cout << e << std::endl;
}