diff --git a/dune/elasticity/materials/micromorphicallyrelaxedenergy.hh b/dune/elasticity/materials/micromorphicallyrelaxedenergy.hh index 6ecbb478c24a0a73b277f32ea5e5ec283e4b6423..67cd55bd3a131c39c0e334d018a433cd946d0687 100644 --- a/dune/elasticity/materials/micromorphicallyrelaxedenergy.hh +++ b/dune/elasticity/materials/micromorphicallyrelaxedenergy.hh @@ -4,6 +4,7 @@ #define DUNE_ELASTICITY_MATERIALS_MICROMORPHICALLYRELAXEDENERGY_HH #include <dune/common/fmatrix.hh> +#include <dune/common/transpose.hh> #include <dune/geometry/quadraturerules.hh> diff --git a/src/quasiconvexity-test-micromorphic.cc b/src/quasiconvexity-test-micromorphic.cc index f8ab4782d36d1a4e1e625fce9dda95765e8aa332..a17dc7dca49660e1bcd3e69e81333a5b0bc76e6a 100644 --- a/src/quasiconvexity-test-micromorphic.cc +++ b/src/quasiconvexity-test-micromorphic.cc @@ -1,5 +1,7 @@ #include <config.h> +#undef HAVE_DUNE_PARMG + // Includes for the ADOL-C automatic differentiation library // Need to come before (almost) all others. #include <adolc/adouble.h> @@ -8,12 +10,9 @@ #include <dune/fufem/utilities/adolcnamespaceinjections.hh> -#undef HAVE_DUNE_PARMG - #include <dune/common/bitsetvector.hh> #include <dune/common/parametertree.hh> #include <dune/common/parametertreeparser.hh> -#include <dune/common/transpose.hh> #include <dune/grid/uggrid.hh> #include <dune/grid/utility/structuredgridfactory.hh> @@ -109,7 +108,6 @@ void writeDisplacement(const Basis& basis, const Vector& deformation, const Fiel vtkWriter.write(filename); } - int main (int argc, char *argv[]) try { // initialize MPI, finalize is done automatically on exit @@ -172,16 +170,16 @@ int main (int argc, char *argv[]) try } else { std::string path = parameterSet.get<std::string>("path"); 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); #if HAVE_DUNE_PARMG - typedef Grid::LevelGridView GridView; + using GridView = Grid::LevelGridView; GridView gridView = grid->levelGridView(0); #else - typedef Grid::LeafGridView GridView; + using GridView = Grid::LeafGridView; GridView gridView = grid->leafGridView(); #endif @@ -210,7 +208,6 @@ int main (int argc, char *argv[]) try BitSetVector<1> dirichletVertices(gridView.size(dim), true); BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices); - // TODO: Do we still need to distinguish between dirichletNodes and dirichletDofs? BitSetVector<dim> dirichletDofs(feBasis.size(), false); constructBoundaryDofs(dirichletBoundary,feBasis,dirichletDofs); diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc index 173bf59f7d607773d4f85caca7fd64deb0edb9c8..d20e5a4e7bb091312f618b091ec925855a2075cb 100644 --- a/src/quasiconvexity-test.cc +++ b/src/quasiconvexity-test.cc @@ -23,9 +23,9 @@ #include <dune/functions/functionspacebases/interpolate.hh> #include <dune/functions/functionspacebases/lagrangebasis.hh> -#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh> #include <dune/functions/functionspacebases/periodicbasis.hh> #include <dune/functions/functionspacebases/powerbasis.hh> +#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh> #include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh> #include <dune/fufem/boundarypatch.hh> @@ -114,7 +114,7 @@ void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldM Vector displacement = nonperiodicX; 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()); std::vector<double> deformationDeterminant(p0Basis.size()); // K = 0.5 F.frobenius_norm2() / F.determinant(); @@ -200,7 +200,7 @@ int main (int argc, char *argv[]) try << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')" << std::endl; - typedef BlockVector<FieldVector<double,dim> > SolutionType; + using Vector = BlockVector<FieldVector<double,dim> >; // parse data file ParameterTree parameterSet; @@ -229,9 +229,9 @@ int main (int argc, char *argv[]) try // /////////////////////////////////////// // 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); @@ -241,12 +241,12 @@ int main (int argc, char *argv[]) try upper = parameterSet.get<FieldVector<double,dim> >("upper"); 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 { std::string path = parameterSet.get<std::string>("path"); 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); @@ -257,10 +257,10 @@ int main (int argc, char *argv[]) try std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl; #if HAVE_DUNE_PARMG - typedef GridType::LevelGridView GridView; + using GridView = Grid::LevelGridView; GridView gridView = grid->levelGridView(0); #else - typedef GridType::LeafGridView GridView; + using GridView = Grid::LeafGridView; GridView gridView = grid->leafGridView(); #endif @@ -348,11 +348,11 @@ int main (int argc, char *argv[]) try } #endif - // ////////////////////////// + //////////////////////////// // Initial iterate - // ////////////////////////// + //////////////////////////// - SolutionType x(feBasis.size()); + Vector x(feBasis.size()); auto powerBasis = makeBasis( gridView, @@ -522,7 +522,7 @@ int main (int argc, char *argv[]) try 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 @@ -532,43 +532,43 @@ int main (int argc, char *argv[]) try #ifdef HAVE_IPOPT // First create an IPOpt base solver - auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,SolutionType> >( + auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,Vector> >( baseTolerance, baseIterations, NumProc::QUIET, "mumps"); #else // 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! - TwoNorm<CorrectionType>* baseNorm = new TwoNorm<CorrectionType>; + TwoNorm<CorrectionType>* baseNorm = new TwoNorm<Vector>; - auto baseSolver = std::make_shared<::LoopSolver<CorrectionType>>(baseSolverStep, - baseIterations, - baseTolerance, - baseNorm, - Solver::QUIET); + auto baseSolver = std::make_shared<::LoopSolver<Vector>>(baseSolverStep, + baseIterations, + baseTolerance, + baseNorm, + Solver::QUIET); #endif // Make pre and postsmoothers - auto presmoother = std::make_shared< TrustRegionGSStep<Matrix, SolutionType> >(); - auto postsmoother = std::make_shared< TrustRegionGSStep<Matrix, SolutionType> >(); + auto presmoother = std::make_shared< TrustRegionGSStep<Matrix, Vector> >(); + 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->setMGType(mu, nu1, nu2); mmgStep->ignoreNodes_ = &dirichletDofs; mmgStep->setBaseSolver(baseSolver); mmgStep->setSmoother(presmoother, postsmoother); - mmgStep->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<SolutionType> >()); + mmgStep->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<Vector> >()); ////////////////////////////////////// // Create the transfer operators ////////////////////////////////////// - using TransferOperatorType = typename TruncatedCompressedMGTransfer<SolutionType>::TransferOperatorType; - std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<SolutionType>>> transferOperators(numLevels-1); + using TransferOperatorType = typename TruncatedCompressedMGTransfer<Vector>::TransferOperatorType; + std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<Vector>>> transferOperators(numLevels-1); // For the restriction operators: FE bases on all levels @@ -604,7 +604,7 @@ int main (int argc, char *argv[]) try assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases[j], *levelBases[j+1]); // Construct the local multigrid transfer matrix - transferOperators[j] = std::make_shared<TruncatedCompressedMGTransfer<SolutionType> >(); + transferOperators[j] = std::make_shared<TruncatedCompressedMGTransfer<Vector> >(); transferOperators[j]->setMatrix(transferMatrix); } @@ -613,13 +613,13 @@ int main (int argc, char *argv[]) try assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases.back(), *levelBases.back()); // 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); } mmgStep->setTransferOperators(transferOperators); - TrustRegionSolver<FEBasis,SolutionType> solver; + TrustRegionSolver<FEBasis,Vector> solver; solver.setup(*grid, assembler, x, @@ -708,7 +708,7 @@ int main (int argc, char *argv[]) try // Careful: The following code computes the initial iterate given by the Python file, // and *assumes* that that file describes the homogeneous deformation! - SolutionType homogeneous(feBasis.size()); + Vector homogeneous(feBasis.size()); Dune::Functions::interpolate(powerBasis, homogeneous, initialDeformation); std::shared_ptr<Elasticity::LocalEnergy<GridView, @@ -720,18 +720,18 @@ int main (int argc, char *argv[]) try LocalADOLCStiffness<GridView, 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 typedef DuneFunctionsBasis<FEBasis> FufemFEBasis; 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 << "Energy at current iterate: " << assembler.computeEnergy(x) << std::endl; - TrustRegionSolver<FEBasis,SolutionType> solver; + TrustRegionSolver<FEBasis,Vector> solver; solver.setup(*grid, &assembler, x,