Skip to content
Snippets Groups Projects
Commit edff98e0 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Remove cosmetic differences between the quasiconvexity test programs

parent 51c88c81
Branches
No related tags found
No related merge requests found
Pipeline #36191 failed
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#define DUNE_ELASTICITY_MATERIALS_MICROMORPHICALLYRELAXEDENERGY_HH #define DUNE_ELASTICITY_MATERIALS_MICROMORPHICALLYRELAXEDENERGY_HH
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/transpose.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
......
#include <config.h> #include <config.h>
#undef HAVE_DUNE_PARMG
// Includes for the ADOL-C automatic differentiation library // Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others. // Need to come before (almost) all others.
#include <adolc/adouble.h> #include <adolc/adouble.h>
...@@ -8,12 +10,9 @@ ...@@ -8,12 +10,9 @@
#include <dune/fufem/utilities/adolcnamespaceinjections.hh> #include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#undef HAVE_DUNE_PARMG
#include <dune/common/bitsetvector.hh> #include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh> #include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh> #include <dune/common/parametertreeparser.hh>
#include <dune/common/transpose.hh>
#include <dune/grid/uggrid.hh> #include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh> #include <dune/grid/utility/structuredgridfactory.hh>
...@@ -109,7 +108,6 @@ void writeDisplacement(const Basis& basis, const Vector& deformation, const Fiel ...@@ -109,7 +108,6 @@ void writeDisplacement(const Basis& basis, const Vector& deformation, const Fiel
vtkWriter.write(filename); vtkWriter.write(filename);
} }
int main (int argc, char *argv[]) try int main (int argc, char *argv[]) try
{ {
// initialize MPI, finalize is done automatically on exit // initialize MPI, finalize is done automatically on exit
...@@ -172,16 +170,16 @@ int main (int argc, char *argv[]) try ...@@ -172,16 +170,16 @@ int main (int argc, char *argv[]) try
} else { } else {
std::string path = parameterSet.get<std::string>("path"); std::string path = parameterSet.get<std::string>("path");
std::string gridFile = parameterSet.get<std::string>("gridFile"); std::string gridFile = parameterSet.get<std::string>("gridFile");
grid = std::shared_ptr<Grid>(GmshReader<Grid>::read(path + "/" + gridFile)); grid = GmshReader<Grid>::read(path + "/" + gridFile);
} }
grid->globalRefine(numLevels-1); grid->globalRefine(numLevels-1);
#if HAVE_DUNE_PARMG #if HAVE_DUNE_PARMG
typedef Grid::LevelGridView GridView; using GridView = Grid::LevelGridView;
GridView gridView = grid->levelGridView(0); GridView gridView = grid->levelGridView(0);
#else #else
typedef Grid::LeafGridView GridView; using GridView = Grid::LeafGridView;
GridView gridView = grid->leafGridView(); GridView gridView = grid->leafGridView();
#endif #endif
...@@ -210,7 +208,6 @@ int main (int argc, char *argv[]) try ...@@ -210,7 +208,6 @@ int main (int argc, char *argv[]) try
BitSetVector<1> dirichletVertices(gridView.size(dim), true); BitSetVector<1> dirichletVertices(gridView.size(dim), true);
BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices); BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
// TODO: Do we still need to distinguish between dirichletNodes and dirichletDofs?
BitSetVector<dim> dirichletDofs(feBasis.size(), false); BitSetVector<dim> dirichletDofs(feBasis.size(), false);
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletDofs); constructBoundaryDofs(dirichletBoundary,feBasis,dirichletDofs);
......
...@@ -23,9 +23,9 @@ ...@@ -23,9 +23,9 @@
#include <dune/functions/functionspacebases/interpolate.hh> #include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh> #include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/functionspacebases/periodicbasis.hh> #include <dune/functions/functionspacebases/periodicbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh> #include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh> #include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
#include <dune/fufem/boundarypatch.hh> #include <dune/fufem/boundarypatch.hh>
...@@ -114,7 +114,7 @@ void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldM ...@@ -114,7 +114,7 @@ void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldM
Vector displacement = nonperiodicX; Vector displacement = nonperiodicX;
displacement += affineDisplacement; displacement += affineDisplacement;
// Compute the determinant per element, for better understanding of the SolutionType // Compute the determinant per element, for better understanding of the solution
Functions::LagrangeBasis<GridView,0> p0Basis(basis.gridView()); Functions::LagrangeBasis<GridView,0> p0Basis(basis.gridView());
std::vector<double> deformationDeterminant(p0Basis.size()); std::vector<double> deformationDeterminant(p0Basis.size());
// K = 0.5 F.frobenius_norm2() / F.determinant(); // K = 0.5 F.frobenius_norm2() / F.determinant();
...@@ -200,7 +200,7 @@ int main (int argc, char *argv[]) try ...@@ -200,7 +200,7 @@ int main (int argc, char *argv[]) try
<< std::endl << "sys.path.append(os.getcwd() + '/../../problems/')" << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
<< std::endl; << std::endl;
typedef BlockVector<FieldVector<double,dim> > SolutionType; using Vector = BlockVector<FieldVector<double,dim> >;
// parse data file // parse data file
ParameterTree parameterSet; ParameterTree parameterSet;
...@@ -229,9 +229,9 @@ int main (int argc, char *argv[]) try ...@@ -229,9 +229,9 @@ int main (int argc, char *argv[]) try
// /////////////////////////////////////// // ///////////////////////////////////////
// Create the grid // Create the grid
// /////////////////////////////////////// // ///////////////////////////////////////
typedef UGGrid<dim> GridType; using Grid = UGGrid<dim>;
std::shared_ptr<GridType> grid; std::shared_ptr<Grid> grid;
FieldVector<double,dim> lower(0), upper(1); FieldVector<double,dim> lower(0), upper(1);
...@@ -241,12 +241,12 @@ int main (int argc, char *argv[]) try ...@@ -241,12 +241,12 @@ int main (int argc, char *argv[]) try
upper = parameterSet.get<FieldVector<double,dim> >("upper"); upper = parameterSet.get<FieldVector<double,dim> >("upper");
std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements"); std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
grid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elements); grid = StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, elements);
} else { } else {
std::string path = parameterSet.get<std::string>("path"); std::string path = parameterSet.get<std::string>("path");
std::string gridFile = parameterSet.get<std::string>("gridFile"); std::string gridFile = parameterSet.get<std::string>("gridFile");
grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile)); grid = GmshReader<Grid>::read(path + "/" + gridFile);
} }
grid->globalRefine(numLevels-1); grid->globalRefine(numLevels-1);
...@@ -257,10 +257,10 @@ int main (int argc, char *argv[]) try ...@@ -257,10 +257,10 @@ int main (int argc, char *argv[]) try
std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl; std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
#if HAVE_DUNE_PARMG #if HAVE_DUNE_PARMG
typedef GridType::LevelGridView GridView; using GridView = Grid::LevelGridView;
GridView gridView = grid->levelGridView(0); GridView gridView = grid->levelGridView(0);
#else #else
typedef GridType::LeafGridView GridView; using GridView = Grid::LeafGridView;
GridView gridView = grid->leafGridView(); GridView gridView = grid->leafGridView();
#endif #endif
...@@ -352,7 +352,7 @@ int main (int argc, char *argv[]) try ...@@ -352,7 +352,7 @@ int main (int argc, char *argv[]) try
// Initial iterate // Initial iterate
//////////////////////////// ////////////////////////////
SolutionType x(feBasis.size()); Vector x(feBasis.size());
auto powerBasis = makeBasis( auto powerBasis = makeBasis(
gridView, gridView,
...@@ -522,7 +522,7 @@ int main (int argc, char *argv[]) try ...@@ -522,7 +522,7 @@ int main (int argc, char *argv[]) try
auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy); auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy);
auto assembler = std::make_shared<Elasticity::FEAssembler<FEBasis,SolutionType> >(feBasis, localADOLCStiffness); auto assembler = std::make_shared<Elasticity::FEAssembler<FEBasis,Vector> >(feBasis, localADOLCStiffness);
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// Create a trust-region solver // Create a trust-region solver
...@@ -532,43 +532,43 @@ int main (int argc, char *argv[]) try ...@@ -532,43 +532,43 @@ int main (int argc, char *argv[]) try
#ifdef HAVE_IPOPT #ifdef HAVE_IPOPT
// First create an IPOpt base solver // First create an IPOpt base solver
auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,SolutionType> >( auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,Vector> >(
baseTolerance, baseTolerance,
baseIterations, baseIterations,
NumProc::QUIET, NumProc::QUIET,
"mumps"); "mumps");
#else #else
// First create a Gauss-seidel base solver // First create a Gauss-seidel base solver
TrustRegionGSStep<Matrix, CorrectionType>* baseSolverStep = new TrustRegionGSStep<Matrix, SolutionType>; TrustRegionGSStep<Matrix, Vector>* baseSolverStep = new TrustRegionGSStep<Matrix, Vector>;
// Hack: the two-norm may not scale all that well, but it is fast! // Hack: the two-norm may not scale all that well, but it is fast!
TwoNorm<CorrectionType>* baseNorm = new TwoNorm<CorrectionType>; TwoNorm<CorrectionType>* baseNorm = new TwoNorm<Vector>;
auto baseSolver = std::make_shared<::LoopSolver<CorrectionType>>(baseSolverStep, auto baseSolver = std::make_shared<::LoopSolver<Vector>>(baseSolverStep,
baseIterations, baseIterations,
baseTolerance, baseTolerance,
baseNorm, baseNorm,
Solver::QUIET); Solver::QUIET);
#endif #endif
// Make pre and postsmoothers // Make pre and postsmoothers
auto presmoother = std::make_shared< TrustRegionGSStep<Matrix, SolutionType> >(); auto presmoother = std::make_shared< TrustRegionGSStep<Matrix, Vector> >();
auto postsmoother = std::make_shared< TrustRegionGSStep<Matrix, SolutionType> >(); auto postsmoother = std::make_shared< TrustRegionGSStep<Matrix, Vector> >();
auto mmgStep = std::make_shared<MonotoneMGStep<Matrix, SolutionType> >(); auto mmgStep = std::make_shared<MonotoneMGStep<Matrix, Vector> >();
mmgStep->setVerbosity(NumProc::FULL); mmgStep->setVerbosity(NumProc::FULL);
mmgStep->setMGType(mu, nu1, nu2); mmgStep->setMGType(mu, nu1, nu2);
mmgStep->ignoreNodes_ = &dirichletDofs; mmgStep->ignoreNodes_ = &dirichletDofs;
mmgStep->setBaseSolver(baseSolver); mmgStep->setBaseSolver(baseSolver);
mmgStep->setSmoother(presmoother, postsmoother); mmgStep->setSmoother(presmoother, postsmoother);
mmgStep->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<SolutionType> >()); mmgStep->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<Vector> >());
////////////////////////////////////// //////////////////////////////////////
// Create the transfer operators // Create the transfer operators
////////////////////////////////////// //////////////////////////////////////
using TransferOperatorType = typename TruncatedCompressedMGTransfer<SolutionType>::TransferOperatorType; using TransferOperatorType = typename TruncatedCompressedMGTransfer<Vector>::TransferOperatorType;
std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<SolutionType>>> transferOperators(numLevels-1); std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<Vector>>> transferOperators(numLevels-1);
// For the restriction operators: FE bases on all levels // For the restriction operators: FE bases on all levels
...@@ -604,7 +604,7 @@ int main (int argc, char *argv[]) try ...@@ -604,7 +604,7 @@ int main (int argc, char *argv[]) try
assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases[j], *levelBases[j+1]); assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases[j], *levelBases[j+1]);
// Construct the local multigrid transfer matrix // Construct the local multigrid transfer matrix
transferOperators[j] = std::make_shared<TruncatedCompressedMGTransfer<SolutionType> >(); transferOperators[j] = std::make_shared<TruncatedCompressedMGTransfer<Vector> >();
transferOperators[j]->setMatrix(transferMatrix); transferOperators[j]->setMatrix(transferMatrix);
} }
...@@ -613,13 +613,13 @@ int main (int argc, char *argv[]) try ...@@ -613,13 +613,13 @@ int main (int argc, char *argv[]) try
assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases.back(), *levelBases.back()); assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases.back(), *levelBases.back());
// Construct the local multigrid transfer matrix // Construct the local multigrid transfer matrix
transferOperators.push_back(std::make_shared<TruncatedCompressedMGTransfer<SolutionType> >()); transferOperators.push_back(std::make_shared<TruncatedCompressedMGTransfer<Vector> >());
transferOperators.back()->setMatrix(transferMatrix); transferOperators.back()->setMatrix(transferMatrix);
} }
mmgStep->setTransferOperators(transferOperators); mmgStep->setTransferOperators(transferOperators);
TrustRegionSolver<FEBasis,SolutionType> solver; TrustRegionSolver<FEBasis,Vector> solver;
solver.setup(*grid, solver.setup(*grid,
assembler, assembler,
x, x,
...@@ -708,7 +708,7 @@ int main (int argc, char *argv[]) try ...@@ -708,7 +708,7 @@ int main (int argc, char *argv[]) try
// Careful: The following code computes the initial iterate given by the Python file, // Careful: The following code computes the initial iterate given by the Python file,
// and *assumes* that that file describes the homogeneous deformation! // and *assumes* that that file describes the homogeneous deformation!
SolutionType homogeneous(feBasis.size()); Vector homogeneous(feBasis.size());
Dune::Functions::interpolate(powerBasis, homogeneous, initialDeformation); Dune::Functions::interpolate(powerBasis, homogeneous, initialDeformation);
std::shared_ptr<Elasticity::LocalEnergy<GridView, std::shared_ptr<Elasticity::LocalEnergy<GridView,
...@@ -720,18 +720,18 @@ int main (int argc, char *argv[]) try ...@@ -720,18 +720,18 @@ int main (int argc, char *argv[]) try
LocalADOLCStiffness<GridView, LocalADOLCStiffness<GridView,
FEBasis::LocalView::Tree::FiniteElement, FEBasis::LocalView::Tree::FiniteElement,
SolutionType> localADOLCStiffness(energy.get()); Vector> localADOLCStiffness(energy.get());
// dune-fufem-style FE basis for the transition from dune-fufem to dune-functions // dune-fufem-style FE basis for the transition from dune-fufem to dune-functions
typedef DuneFunctionsBasis<FEBasis> FufemFEBasis; typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
FufemFEBasis fufemFEBasis(feBasis); FufemFEBasis fufemFEBasis(feBasis);
FEAssembler<FufemFEBasis,SolutionType> assembler(fufemFEBasis, &localADOLCStiffness); FEAssembler<FufemFEBasis,Vector> assembler(fufemFEBasis, &localADOLCStiffness);
std::cout << "Homogeneous energy: " << assembler.computeEnergy(homogeneous) << std::endl; std::cout << "Homogeneous energy: " << assembler.computeEnergy(homogeneous) << std::endl;
std::cout << "Energy at current iterate: " << assembler.computeEnergy(x) << std::endl; std::cout << "Energy at current iterate: " << assembler.computeEnergy(x) << std::endl;
TrustRegionSolver<FEBasis,SolutionType> solver; TrustRegionSolver<FEBasis,Vector> solver;
solver.setup(*grid, solver.setup(*grid,
&assembler, &assembler,
x, x,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment