Skip to content

GitLab

Commit 140820bc authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

More stuff from dune-tnnmg-examples 

I am not really sure what this stuff is all about..
parent 23a38ca9
Hide whitespace changes
Inline Side-by-side
doc/manual/examples.cc 0 → 100755
View file @ 140820bc
#include <config.h>
#define DIMENSION 2
#define ALUGRID 1
//#define UGGRID 1
#define EMBEDDED_PYTHON 1
#undef EMBEDDED_PYTHON
#define FE_VERBOSE
// disable embedded python if python was not found
#ifndef HAVE_PYTHON
#undef EMBEDDED_PYTHON
#endif
#if EMBEDDED_PYTHON
#include <Python.h>
#endif
#include <fstream>
#include <iostream>
#include <iostream>
#include <vector>
#include <tr1/memory>
// dune includes ******************************************************
#ifdef UGGRID
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#endif
#ifdef ALUGRID
#include <dune/grid/alugrid.hh>
#include <dune/grid/io/file/dgfparser/dgfalu.hh>
#endif
#include <dune/common/configparser.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/grid/common/genericreferenceelements.hh>
// dune-fufem includes ******************************************************
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/prolongboundarypatch.hh>
#if EMBEDDED_PYTHON
#include <dune/fufem/functions/pythonfunction.hh>
#endif
#include <dune/fufem/functions/functions.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/functiontools/namedfunctionmap.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/localassemblers/lumpedmassassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#if HAVE_AMIRAMESH
#include <dune/fufem/functiontools/amirameshbasiswriter.hh>
#endif
#include <dune/fufem/functiontools/vtkbasiswriter.hh>
// dune-solvers includes ***************************************************
#include <dune/solvers/norms/norm.hh>
#include <dune/solvers/norms/sumnorm.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/diagnorm.hh>
#include <dune/solvers/norms/fullnorm.hh>
#include <dune/solvers/norms/blocknorm.hh>
// dune-tnnmg includes *******************************************************
#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
#include <dune/tnnmg/problem-classes/nonlinearity.hh>
#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
#include <dune/tnnmg/nonlinearities/shiftednonlinearity.hh>
#include <dune/tnnmg/iterationsteps/genericnonlinearjacobi.hh>
#include <dune/tnnmg/solvers/scalartnnmg.hh>
#include "laplace.hh"
template <int dim, int N>
class ConstantFunction :
public Dune::FunctionBase<double,double,dim,N>
{
public:
ConstantFunction(double c):
c_(c)
{}
double eval(int comp, const Dune::FieldVector<double,dim>& x) const
{
return c_;
}
void evalall(const Dune::FieldVector<double,dim>& x, Dune::FieldVector<double,N>& y) const
{
y = eval(0,x);
}
double c_;
};
template <int dim, int N>
class AbsFunction :
public Dune::FunctionBase<double,double,dim,N>
{
public:
AbsFunction()
{}
double eval(int comp, const Dune::FieldVector<double,dim>& x) const
{
return fabs(x[0]);
}
void evalall(const Dune::FieldVector<double,dim>& x, Dune::FieldVector<double,N>& y) const
{
y = eval(0,x);
}
};
int main (int argc, char *argv[]) try
{
#ifdef EMBEDDED_PYTHON
EmbeddedPython::start();
#endif
Dune::ConfigParser parset;
// parse command line once to obtain a possible parset argument
// parse parameter file
// parse command line a second time in order to allow overwriting parameters
parset.parseCmd(argc, argv);
parset.parseFile(parset.get("file", "adaptive_laplace.parset"));
parset.parseCmd(argc, argv);
// The grid dimension
const int dim = DIMENSION;
// define grid type
#ifdef UGGRID
typedef Dune::UGGrid<dim> GridType;
#endif
#ifdef ALUGRID
typedef Dune::ALUSimplexGrid<dim,dim> GridType;
#endif
typedef NamedFunctionMap<GridType::Codim<0>::Geometry::GlobalCoordinate,Dune::FieldVector<double,1> > Functions;
typedef Functions::Function Function;
Functions functions;
functions["one"] = new ConstantFunction<dim,1>(1.0);
functions["zero"] = new ConstantFunction<dim,1>(0.0);
functions["abs"] = new AbsFunction<dim,1>;
// Generate the grid
Dune::GridPtr<GridType> gridptr;
if (dim==2)
gridptr=Dune::GridPtr<GridType>(parset.get("grid.dgffile2d", "coarse.dgf").c_str());
else if (dim==3)
gridptr=Dune::GridPtr<GridType>(parset.get("grid.dgffile3d", "coarse.dgf").c_str());
GridType& grid = *gridptr;
// grid specific settings
#ifdef UGGRID
grid.setRefinementType(GridType::LOCAL);
grid.setClosureType(GridType::NONE);
#endif
#ifdef ALUGRID
if (parset.hasKey("grid.restore"))
{
double time;
grid.readGrid<Dune::xdr>(parset.get<std::string>("grid.backupfile"), time);
std::cout << "Grid restored. Time was " << time << std::endl;
}
else
#endif
{
//initial refinement of grid
for (int i=0; i<parset.get("grid.refine", 0) or i<parset.get("laplace.refinement.minlevel",0); ++i)
grid.globalRefine(1);
std::cout << "Grid refined." << std::endl;
}
#ifdef EMBEDDED_PYTHON
std::string importDictName = formatString(std::string("importAsDuneFunction%1dD"), dim);
if (parset.hasKey("python.inline"))
{
EmbeddedPython::run(parset.get("python.inline", ""));
EmbeddedPython::importFunctionsFromModule(functions, "__main__", importDictName);
}
if (parset.hasKey("python.module"))
EmbeddedPython::importFunctionsFromModule(functions, parset.get("python.module", "laplace"), importDictName);
#endif
LaplaceProblem<GridType> laplaceProblem(grid, parset.sub("laplace"), functions);
bool refined;
do
{
laplaceProblem.assemble();
laplaceProblem.solve();
refined = laplaceProblem.adapt();
}
while(refined);
#ifdef EMBEDDED_PYTHON
EmbeddedPython::stop();
#endif
// delete created functions
Functions::iterator it = functions.begin();
Functions::iterator end = functions.end();
for(; it!=end; ++it)
delete it->second;
return 0;
}
catch (Dune::Exception e)
{
std::cout << e << std::endl;
}
doc/manual/laplace.hh 0 → 100755
View file @ 140820bc
#ifndef TNNMG_EXAMPLE_LAPLACE_HH
#define TNNMG_EXAMPLE_LAPLACE_HH
#include <fstream>
#include <iostream>
#include <iostream>
#include <vector>
#include <tr1/memory>
// dune includes ******************************************************
#include <dune/common/parametertree.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/grid/utility/grapedataioformattypes.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/grid/common/genericreferenceelements.hh>
// dune-fufem includes ******************************************************
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/prolongboundarypatch.hh>
#if EMBEDDED_PYTHON
#include <dune/fufem/functions/pythonfunction.hh>
#endif
#include <dune/fufem/functions/functions.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/functiontools/namedfunctionmap.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/functiontools/gridfunctionadaptor.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
#include <dune/fufem/assemblers/localassemblers/lumpedmassassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#if HAVE_AMIRAMESH
#include <dune/fufem/functiontools/amirameshbasiswriter.hh>
#endif
#include <dune/fufem/functiontools/vtkbasiswriter.hh>
// dune-solvers includes ***************************************************
#include <dune/solvers/norms/norm.hh>
#include <dune/solvers/norms/sumnorm.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/diagnorm.hh>
#include <dune/solvers/norms/fullnorm.hh>
#include <dune/solvers/norms/blocknorm.hh>
#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/solvers/loopsolver.hh>
// dune-tnnmg includes *******************************************************
//#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
#include <dune/tnnmg/problem-classes/nonlinearity.hh>
#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
template<class T1>
std::string formatString(const std::string format, const T1& t1)
{
char buffer[1000];
sprintf(buffer, format.c_str(), t1);
return std::string(buffer);
}
template<class T1, class T2>
std::string formatString(const std::string format, const T1& t1, const T2& t2)
{
char buffer[1000];
sprintf(buffer, format.c_str(), t1, t2);
return std::string(buffer);
}
template <class GridType>
class LaplaceProblem
{
public:
static const int dim = GridType::dimension;
static const int block_size = 1;
typedef Dune::FieldMatrix<double,block_size,block_size> FineLocalMatrix;
typedef Dune::FieldVector<double,block_size> FineLocalVector;
typedef Dune::BCRSMatrix< FineLocalMatrix > FineMatrixType;
typedef Dune::BlockVector< FineLocalVector > FineVectorType;
/* in the linear Laplace problem there is no nonlinearity ;-) */
typedef ZeroNonlinearity<FineLocalVector,FineLocalMatrix> NonlinearityType;
/* we set the problemtypes:
the ConvexProblem depends on the NonlinearityType and the type of the fine grid matrix*/
typedef ConvexProblem<NonlinearityType, FineMatrixType> ConvexProblemType;
typedef BlockNonlinearTNNMGProblem<ConvexProblemType> TNNMGProblemType;
typedef typename TNNMGProblemType::Linearization::VectorType CoarseVectorType;
typedef typename TNNMGProblemType::Linearization::MatrixType CoarseMatrixType;
/* use a linear Gauss-Seidel as fine grid smoother */
typedef GenericNonlinearGS<TNNMGProblemType> FineSmoother;
typedef TruncatedNonsmoothNewtonMultigrid<TNNMGProblemType, FineSmoother> TNNMGType;
/* define the linear coarse grid solver type:
use linear Gauss-Seidel as smoother for linear multigrid */
typedef BlockGSStep<CoarseMatrixType,CoarseVectorType> CoarseSmoother;
typedef ScalarTNNMG::Transfer Transfer;
typedef NamedFunctionMap<typename GridType::template Codim<0>::Geometry::GlobalCoordinate,Dune::FieldVector<double,1> > FunctionMap;
typedef typename FunctionMap::Function Function;
typedef P1NodalBasis<typename GridType::LeafGridView> NonconformingBasis;
typedef ConformingBasis<NonconformingBasis> Basis;
LaplaceProblem(GridType& grid, Dune::ParameterTree& parset, const FunctionMap& functions) :
grid_(grid),
parset_(parset),
functions_(functions),
ncBasis_(grid.leafView()),
basis_(ncBasis_)
{
// create energy norm
norm_ = NormPointer(new EnergyNorm<FineMatrixType,FineVectorType>(stiffnessMatrix_));
// create a zero- nonlinearity, i.e. no nonlinearity ;-)
phi_ = NonlinearityPointer(new NonlinearityType());
// initialize solution vector
u_.resize(basis_.size());
u_ = 0.0;
refinementStep_ = 0;
typedef BoundaryPatchSegmentIndexProperty<typename GridType::LevelGridView> Property;
coarseBoundaryPatch_.setup(grid.levelView(0));
coarseBoundaryPatch_.insertFacesByProperty(Property(1));
}
~LaplaceProblem()
{
for(int l=0; l<transfer_.size(); ++l)
delete transfer_[l];
delete norm_;
}
template <class AssembleBasis>
void assembleRHS(const AssembleBasis& basis, const int rhsIntegrationOrder)
{
FunctionalAssembler<AssembleBasis> assembler(basis);
{
std::cout << "assembling right hand side" << std::endl;
const Function& rhsFunction = *functions_.getFunction(parset_.get("rhs.function", "one"));
L2FunctionalAssembler<GridType> l2Assembler(rhsFunction, rhsIntegrationOrder);
assembler.assemble(l2Assembler, rhs_);
}
{
std::cout << "assembling mass weight vector" << std::endl;
L2FunctionalAssembler<GridType> l2Assembler(*functions_.getFunction("one"), 1);
assembler.assemble(l2Assembler, massWeights_);
}
}
// assemble the whole problem
void assemble()
{
// setup boundary patch and boundary dofs
PatchProlongator<GridType>::prolong(coarseBoundaryPatch_, boundaryPatch_);
// setup boundary dofs
constructBoundaryDofs(boundaryPatch_, basis_, isBoundary_);
// assemble right hand side with given integration order
assembleRHS(basis_, parset_.get("rhs.integrationorder", 2));
// create a global operator assembler to be used for the matrices
OperatorAssembler<Basis, Basis> assembler(basis_, basis_);
typedef typename Basis::LocalFiniteElement FE;
{
// assemble stiffness matrix
std::cout << "assembling stiffness matrix" << std::endl;
LaplaceAssembler<GridType,FE,FE> laplaceAssembler;
assembler.assemble(laplaceAssembler, stiffnessMatrix_);
}
if (parset_.get("masslumping",true))
{
// assemble lumped mass matrix
std::cout << "assembling lumped mass matrix" << std::endl;
LumpedMassAssembler<GridType,FE,FE> lumpedMassAssembler;
assembler.assemble(lumpedMassAssembler, massMatrix_, true);
}
else
{
// assemble mass matrix
std::cout << "assembling mass matrix" << std::endl;
MassAssembler<GridType,FE,FE> massAssembler;
assembler.assemble(massAssembler, massMatrix_);
}
// multiply matrices by given coefficients and add them
stiffnessMatrix_ *= parset_.get("stiffnessfactor", 1.0);
massMatrix_ *= parset_.get("massfactor", 1.0);
stiffnessMatrix_ += massMatrix_;
// delete old transfer operators
for(int l=0; l<transfer_.size(); ++l)
delete transfer_[l];
// create new transfer operators
transfer_.resize(grid_.maxLevel());
for(int l=0; l<transfer_.size(); ++l)
transfer_[l] = new Transfer;
if (parset_.get("refinement.type", "local") == "global")
{
// assemble transfer operators on uniformly refined grid
// this is only for demonstration reasons
// the all at once aasembler below does also
// work for uniformly refined grids
for(int l=0; l<transfer_.size(); ++l)
transfer_[l]->setup(grid_, l,l+1);
}
else
{
// assemble transfer operator hierarchy on adaptively refined grid
TransferOperatorAssembler<GridType> transferOperatorAssembler(grid_);
transferOperatorAssembler.assembleOperatorPointerHierarchy(transfer_);
}
// compute entries of solution vector for dirichlet nodes
// as interpolation of given dirichlet function
Functions::interpolate(basis_, u_, *functions_.getFunction(parset_.get("dirichletfunction", "zero")),isBoundary_);
// GenericGridFunctionInterpolator<Basis,FineVectorType,Function>::interpolate(basis_, u_, *functions_.getFunction(parset_.get("dirichletfunction", "zero")), isBoundary_);
}
void solve()
{
std::cout << "u.size(): " << u_.size() << std::endl;
ConvexProblemType P(1.0, stiffnessMatrix_, 0.0, massWeights_, *phi_, rhs_, u_);
TNNMGProblemType tnnmgProblem(parset_.sub("tnnmg"), P);
/* the "nonlinear" or rather fine grid smoother */
FineSmoother fine_smoother;
/* the linear or coarse smoothers */
CoarseSmoother presmoother, postsmoother;
/* the linear (coarse grid) iterative solver step */
MultigridStep<CoarseMatrixType, CoarseVectorType> mgStep;
{
mgStep.setTransferOperators(transfer_);
mgStep.setNumberOfLevels(mgStep.mgTransfer_.size()+1);
mgStep.setSmoother(&presmoother, &postsmoother);
/* use one simple smoother iteration as base solver */
CoarseSmoother* baseSolverStep = new CoarseSmoother;
EnergyNorm<CoarseMatrixType, CoarseVectorType>* baseEnergyNorm = new EnergyNorm<CoarseMatrixType, CoarseVectorType>(*baseSolverStep);
mgStep.basesolver_ = new LoopSolver<CoarseVectorType>(baseSolverStep,
1,
parset_.get("mmg.tol", 1e-12),
baseEnergyNorm,
Solver::QUIET);
mgStep.setMGType(1, 3, 3);
mgStep.setNumberOfLevels(transfer_.size()+1);
}
Dune::BitSetVector<1> ignoreNodes(isBoundary_);
for (int i=0; i<ignoreNodes.size(); ++i)
{
if (basis_.isConstrained(i))
ignoreNodes[i].set();
}
TNNMGType tnnmgstep(mgStep, fine_smoother);
tnnmgstep.setProblem(P.u, tnnmgProblem);
tnnmgstep.ignoreNodes_ = &ignoreNodes;
tnnmgstep.setSmoothingSteps(1, 1, 0);
LoopSolver<VectorType> tnnmgsolver(&tnnmgstep,
parset_.get("tnnmg.maxiter", 100),
parset_.get("tnnmg.tol", 1e-12),
norm_,
Solver::FULL);
tnnmgsolver.preprocess();
tnnmgsolver.solve();
u_ = tnnmgstep.getSol();
}
void writeSolution()
{
#if HAVE_AMIRAMESH
if (parset_.get("solutionfile_am", "") != "")
{
AmiraMeshBasisWriter<Basis> writer(basis_);
writer.addP1Interpolation(u_);
writer.write(formatString(parset_.get("solutionfile_am", ""), refinementStep_));
}
#endif
if (parset_.get("solutionfile_vtk", "") != "")
{
VTKBasisWriter<Basis> writer(basis_);
writer.addP1Interpolation(u_);
writer.write(formatString(parset_.get("solutionfile_vtk", ""), refinementStep_));
}
}
bool adapt()
{
if (parset_.get("refinement.maxlevel",8) > grid_.maxLevel())
{
GridFunctionAdaptor<Basis> adaptor(basis_, true);
grid_.globalRefine(1);
ncBasis_.update(grid_.leafView());
basis_.update(grid_.leafView());
adaptor.adapt(u_);
return true;
}
return false;
}
private:
// const Dune::ParameterTree& parset_;
Dune::ParameterTree& parset_;
GridType& grid_;
const FunctionMap& functions_;
int refinementStep_;
NonconformingBasis ncBasis_;
Basis basis_;
BoundaryPatchBase< typename Basis::GridView > boundaryPatch_;
BoundaryPatchBase< typename GridType::LevelGridView > coarseBoundaryPatch_;
FineMatrixType stiffnessMatrix_;
FineMatrixType massMatrix_;
FineVectorType u_;
FineVectorType rhs_;
FineVectorType massWeights_;
typedef typename Dune::shared_ptr< NonlinearityType > NonlinearityPointer;
NonlinearityPointer phi_;
Dune::BitSetVector<1> isBoundary_;
typedef Norm<FineVectorType>* NormPointer;
NormPointer norm_;
std::vector<Transfer*> transfer_;
};
#endif
doc/manual/src/cube_centered.dgf 0 → 100755
View file @ 140820bc
DGF
Vertex % the verticies of the grid
-1 -1 1 % vertex 0
1 -1 1 % vertex 1
-1 -1 -1 % vertex 4