Select Git revision
linelast.cc
Forked from
agnumpde / dune-elasticity
414 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
linelast.cc 10.83 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/istl/io.hh>
#include <dune/fufem/functiontools/gridfunctionadaptor.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/sampleonbitfield.hh>
#include <dune/fufem/boundarypatchprolongator.hh>
#include <dune/fufem/readbitfield.hh>
#include <dune/fufem/estimators/fractionalmarking.hh>
#include <dune/fufem/estimators/hierarchicalestimator.hh>
#ifdef HAVE_IPOPT
#include <dune/solvers/solvers/quadraticipopt.hh>
#endif
#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#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("linelast.parset", parameterSet);
// read solver settings
const int minLevel = parameterSet.get<int>("minLevel");
const int maxLevel = parameterSet.get<int>("maxLevel");
const int numIt = parameterSet.get<int>("numIt");
const int nu1 = parameterSet.get<int>("nu1");
const int nu2 = parameterSet.get<int>("nu2");
const int mu = parameterSet.get<int>("mu");
const int baseIt = parameterSet.get<int>("baseIt");
const double tolerance = parameterSet.get<double>("tolerance");
const double baseTolerance = parameterSet.get<double>("baseTolerance");
const bool nestedIteration = parameterSet.get<int>("nestedIteration");
const bool paramBoundaries = parameterSet.get<int>("paramBoundaries");
const double threshold = parameterSet.get<double>("threshold");
// read problem settings
std::string path = parameterSet.get<std::string>("path");
std::string gridFile = parameterSet.get<std::string>("gridFile");
std::string parFile = parameterSet.get("parFile", "");
std::string dirichletFile = parameterSet.get<std::string>("dirichletFile");
std::string dirichletValuesFile = parameterSet.get<std::string>("dirichletValuesFile");
// /////////////////////////////
// Generate the grid
// /////////////////////////////
typedef UGGrid<dim> GridType;
typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;
typedef BoundaryPatch<GridType::LeafGridView> LeafBoundaryPatch;
GridType grid;
grid.setRefinementType(GridType::COPY);
if (paramBoundaries)
Dune::AmiraMeshReader<GridType>::read(grid, path + gridFile, path + parFile);
else
Dune::AmiraMeshReader<GridType>::read(grid, path + gridFile);
// Read Dirichlet boundary values
std::vector<BitSetVector<dim> > dirichletNodes(maxLevel+1);
std::vector<LevelBoundaryPatch> dirichletBoundary(maxLevel+1);
dirichletBoundary[0].setup(grid.levelView(0));
readBoundaryPatch<GridType>(dirichletBoundary[0], path + dirichletFile);
std::vector<VectorType> dirichletValues(maxLevel+1);
dirichletValues[0].resize(grid.size(0, dim));
AmiraMeshReader<GridType>::readFunction(dirichletValues[0], path + dirichletValuesFile);
for (int i=0; i<minLevel; i++)
grid.globalRefine(1);
//int maxlevel = grid.maxLevel();
// initial solution
VectorType x(grid.size(grid.maxLevel(), dim));
x = 0;
// //////////////////////////////////////////////////
// Refinement loop
// //////////////////////////////////////////////////
while (true) {
// Determine Dirichlet dofs
BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary);
LeafBoundaryPatch leafDirichletBoundary(grid.leafView());
BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[0], leafDirichletBoundary);
for (int i=0; i<=grid.maxLevel(); i++) {
dirichletNodes[i].resize(grid.size(i,dim));
for (int j=0; j<grid.size(i,dim); j++)
dirichletNodes[i][j] = dirichletBoundary[i].containsVertex(j);
}
sampleOnBitField(grid, dirichletValues, dirichletNodes);
VectorType rhs(grid.size(grid.maxLevel(), dim));
rhs = 0;
for (int i=0; i<rhs.size(); i++)
for (int j=0; j<dim; j++) {
if (dirichletNodes[grid.maxLevel()][i][j])
x[i][j] = dirichletValues[grid.maxLevel()][i][j];
}
// Assemble elasticity problem#
typedef P1NodalBasis<GridType::LeafGridView,double> P1Basis;
P1Basis p1NodalBasis(grid.leafView());
OperatorAssembler<P1Basis,P1Basis> p1Assembler(p1NodalBasis, p1NodalBasis);
StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> p1LocalAssembler(2.5e5, 0.3);
OperatorType stiffnessMatrix;
p1Assembler.assemble(p1LocalAssembler, stiffnessMatrix);
// ///////////////////////////
// Create a solver
// ///////////////////////////
// First create a base solver
#ifdef IPOPT_BASE
#if !defined HAVE_IPOPT
#error You need to have IPOpt installed if you want to use it as the base solver!
#endif
QuadraticIPOptSolver<OperatorType,VectorType> baseSolver;
baseSolver.tolerance_ = baseTolerance;
baseSolver.verbosity_ = Solver::QUIET;
#else // Gauss-Seidel is the base solver
BlockGSStep<OperatorType, VectorType> baseSolverStep;
EnergyNorm<OperatorType,VectorType> baseEnergyNorm(baseSolverStep);
IterativeSolver<VectorType> baseSolver(&baseSolverStep,
baseIt,
baseTolerance,
&baseEnergyNorm,
Solver::QUIET);
#endif
// Make pre and postsmoothers
BlockGSStep<OperatorType, VectorType> presmoother;
BlockGSStep<OperatorType, VectorType> postsmoother;
MultigridStep<OperatorType, VectorType> multigridStep(stiffnessMatrix, x, rhs, grid.maxLevel()+1);
multigridStep.setMGType(1, nu1, nu2);
multigridStep.ignoreNodes_ = &dirichletNodes.back();
multigridStep.basesolver_ = &baseSolver;
multigridStep.setSmoother(&presmoother,&postsmoother);
multigridStep.mgTransfer_.resize(grid.maxLevel());
for (int i=0; i<multigridStep.mgTransfer_.size(); i++) {
CompressedMultigridTransfer<VectorType>* newTransferOp = new CompressedMultigridTransfer<VectorType>;
newTransferOp->setup(grid, i, i+1);
multigridStep.mgTransfer_[i] = newTransferOp;
}
EnergyNorm<OperatorType, VectorType> energyNorm(multigridStep);
LoopSolver<VectorType> solver(&multigridStep,
numIt,
tolerance,
&energyNorm,
Solver::FULL);
solver.preprocess(); multigridStep.preprocess();
// Compute solution
solver.solve();
x = multigridStep.getSol();
// ////////////////////////////////////////////////////////////////////////
// Leave adaptation loop if maximum number of levels has been reached
// ////////////////////////////////////////////////////////////////////////
if (grid.maxLevel() == maxLevel)
break;
// /////////////////////////////////////////////////////////////
// Estimate error and refine grid
// /////////////////////////////////////////////////////////////
HierarchicalEstimator<P2NodalBasis<GridType::LeafGridView>,dim> estimator(grid);
VectorType hierarchicError;
VirtualFunction<FieldVector<double,dim>,FieldVector<double,dim> >* volumeTermDummy=NULL;
VirtualFunction<FieldVector<double,dim>,FieldVector<double,dim> >* neumannTermDummy=NULL;
RefinementIndicator<GridType> indicator(grid);
typedef P2NodalBasis<GridType::LeafGridView>::LocalFiniteElement P2FiniteElement;
StVenantKirchhoffAssembler<GridType,P2FiniteElement,P2FiniteElement>* localStiffness = new StVenantKirchhoffAssembler<GridType,P2FiniteElement,P2FiniteElement>(2.5e5,0.3);
estimator.estimate(x,
volumeTermDummy,
neumannTermDummy,
leafDirichletBoundary,
indicator,
localStiffness);
//estimator.markForRefinement(grid, hierarchicError, threshold);
FractionalMarkingStrategy<GridType>::mark(indicator, grid, 0.2);
// ////////////////////////////////////////////////////
// Refine grid
// ////////////////////////////////////////////////////
GridFunctionAdaptor<P1Basis> adaptorP1(p1NodalBasis,true,true);
grid.preAdapt();
grid.adapt();
grid.postAdapt();
p1NodalBasis.update();
adaptorP1.adapt(x);
// if (paramBoundaries)
// improveGrid(grid,10);
std::cout << "########################################################" << std::endl;
std::cout << " Grid refined" << std::endl;
std::cout << " Level: " << grid.maxLevel()
<< " vertices: " << grid.size(grid.maxLevel(), dim)
<< " elements: " << grid.size(grid.maxLevel(), 0) << std::endl;
std::cout << "########################################################" << std::endl;
}
//std::cout << "Hierarchic Error: " << std::endl << hierarchicError << std::endl;
// Output result
LeafAmiraMeshWriter<GridType> amiramesh;
amiramesh.addLeafGrid(grid,true);
amiramesh.addVertexData(x, grid.leafView());
amiramesh.write("resultGrid");
} catch (Exception e) {
std::cout << e << std::endl;
}