diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index b2dcaf24186480e39e3f542582d4e3f108f03929..335309b407391cccaf4149fbae220c21eb1d3b40 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -27,3 +27,8 @@ add_dependencies(examples navier_stokes.2d convection_diffusion.2d cahn_hilliard.2d) + +if (dune-foamgrid_FOUND) + add_amdis_executable(NAME surface.2d SOURCES surface.cc DIM 2 DOW 3) + add_dependencies(examples surface.2d) +endif () \ No newline at end of file diff --git a/examples/FactoryParametrization.hpp b/examples/FactoryParametrization.hpp new file mode 100644 index 0000000000000000000000000000000000000000..41290b3a98c88ae4d74169b74351696b496303f6 --- /dev/null +++ b/examples/FactoryParametrization.hpp @@ -0,0 +1,104 @@ +#pragma once + +#include <type_traits> +#include <vector> + +#include <dune/grid/common/gridfactory.hh> +#include <amdis/Output.hpp> + +namespace AMDiS +{ + template <class GridType, class Projection> + struct Parametrization + { + using Grid = GridType; + using ctype = typename Grid::ctype; + + static constexpr int dimension = Grid::dimension; + static constexpr int dimensionworld = Grid::dimensionworld; + + using LeafIntersection = typename Grid::LeafIntersection; + + template <int cd> + using Codim = typename Grid::template Codim<cd>; + }; + +} // end namespace AMDiS + + +namespace Dune +{ + + /// \brief A factory class for \ref SimpleGrid. + template <class Grid, class Projection> + class GridFactory<AMDiS::Parametrization<Grid, Projection>> + : public GridFactoryInterface<AMDiS::Parametrization<Grid, Projection>> + { + using GridWrapper = AMDiS::Parametrization<Grid, Projection>; + + /// Type used by the grid for coordinates + using ctype = typename Grid::ctype; + + static constexpr int dim = Grid::dimension; + static constexpr int dow = Grid::dimensionworld; + + using LocalCoordinate = FieldVector<ctype, dim>; + using GlobalCoordinate = FieldVector<ctype, dow>; + + using ProjectionBase = VirtualFunction<LocalCoordinate, GlobalCoordinate>; + + public: + + /// Default constructor + GridFactory() + {} + + /// Insert a vertex into the coarse grid + virtual void insertVertex(GlobalCoordinate const& pos) override + { + GlobalCoordinate new_pos(pos); + Projection::project(new_pos); + + coordinates.push_back(new_pos); + factory.insertVertex(new_pos); + } + + /// Insert an element into the coarse grid. + virtual void insertElement(GeometryType const& type, + std::vector<unsigned int> const& vertices) override + { + assert( type.isSimplex() ); + + std::vector<GlobalCoordinate> corners(vertices.size()); + for (std::size_t i = 0; i < vertices.size(); ++i) + corners[i] = coordinates[vertices[i]]; + + factory.insertElement(type, vertices, std::shared_ptr<ProjectionBase>(new Projection(std::move(corners))) ); + } + + virtual void insertBoundarySegment(std::vector<unsigned int> const& /*vertices*/) override + { + AMDiS::warning("Not implemented!"); + } + + /// Finalize grid creation and hand over the grid. + Grid* create() + { + return factory.createGrid(); + } + + virtual GridWrapper* createGrid() override + { + AMDiS::warning("Should not be created. Use non-virtual method `create()` instead, to create the underlying grid!"); + return nullptr; + } + + private: + + // buffers for the mesh data + std::vector<GlobalCoordinate> coordinates; + + GridFactory<Grid> factory; + }; + +} // end namespace Dune diff --git a/examples/SphereMapping.hpp b/examples/SphereMapping.hpp new file mode 100644 index 0000000000000000000000000000000000000000..3e7b0ec421f47d83972444080b021def63ad696d --- /dev/null +++ b/examples/SphereMapping.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include <algorithm> +#include <array> + +#include <dune/common/function.hh> +#include <dune/common/fvector.hh> + +namespace AMDiS +{ + /** + * \brief Mapping class mapping from a triangle with points on the unit sphere onto the sphere + * with theta in [0, pi] and phi in [0, 2*pi) + */ + template <int dim, int dow, class Radius, class ct = double> + class SphereMapping + : public Dune::VirtualFunction<Dune::FieldVector<ct, dim>, Dune::FieldVector<ct, dow> > + { + using LocalCoordinate = Dune::FieldVector<ct, dim>; + using GlobalCoordinate = Dune::FieldVector<ct, dow>; + + public: + template <class VertexContainer> + SphereMapping(VertexContainer&& vertices) + { + std::copy_n(vertices.begin(), dim+1, vertices_.begin()); + } + + void evaluate(LocalCoordinate const& x, GlobalCoordinate& y) const + { + // calculate global coordinate + auto shapeFunctions = evaluateShapeFunctions(x); + y = 0.0; + for(size_t i = 0; i < dim+1; i++) + for(size_t j = 0; j < dow; j++) + y[j] += vertices_[i][j]*shapeFunctions[i]; + project(y); + } + + GlobalCoordinate evaluateShapeFunctions(LocalCoordinate const& x) const + { + GlobalCoordinate out; + out[0] = 1.0; + for (size_t i=0; i<2; i++) + { + out[0] -= x[i]; + out[i+1] = x[i]; + } + return out; + } + + static void project(GlobalCoordinate& y) + { + // project it on the unit sphere + y /= y.two_norm() / Radius::eval(y); + } + + + private: + std::array<GlobalCoordinate, dim+1> vertices_; + }; + +} // end namespace AMDiS diff --git a/examples/init/surface.dat.2d b/examples/init/surface.dat.2d new file mode 100644 index 0000000000000000000000000000000000000000..aa1dadaec3ac1e9728e8a17ac5644e0bddfae473 --- /dev/null +++ b/examples/init/surface.dat.2d @@ -0,0 +1,14 @@ +surfaceMesh->macro file name: ./macro/sphere_flat.3d +surfaceMesh->global refinements: 1 + +surface->mesh: surfaceMesh + +surface->solver->name: bicgstab +surface->solver->max iteration: 10000 +surface->solver->relative tolerance: 1.e-6 +surface->solver->info: 10 +surface->solver->left precon: ilu + +surface->output[0]->filename: surface.2d +surface->output[0]->name: u +surface->output[0]->output directory: ./output diff --git a/examples/surface.cc b/examples/surface.cc new file mode 100644 index 0000000000000000000000000000000000000000..42df82861d6f0882037327ed5a16c9cd9e530abe --- /dev/null +++ b/examples/surface.cc @@ -0,0 +1,60 @@ +#include <amdis/AMDiS.hpp> +#include <amdis/Integrate.hpp> +#include <amdis/LocalOperators.hpp> +#include <amdis/ProblemStat.hpp> +#include <amdis/common/Literals.hpp> + +#include <dune/grid/albertagrid/albertareader.hh> +#include <dune/foamgrid/foamgrid.hh> + +#include "FactoryParametrization.hpp" +#include "SphereMapping.hpp" + +using namespace AMDiS; +using namespace Dune::Indices; + +using Grid = Dune::FoamGrid<2,3>; +using Basis = LagrangeBasis<typename Grid::LeafGridView, 1>; + +struct UnitRadius +{ + template <class T> + static double eval(T&&) { return 1.0; } +}; + + +// solve the equation -laplace(u) - k^2 u = f on the sphere +int main(int argc, char** argv) +{ + AMDiS::init(argc, argv); + + std::string gridName = Parameters::get<std::string>("surface->mesh").value(); + std::string gridFileName = Parameters::get<std::string>(gridName + "->macro file name").value(); + + using Param = Parametrization<Grid, SphereMapping<2, 3, UnitRadius>>; + Dune::GridFactory<Param> gridFactory; + Dune::AlbertaReader<Param>().readGrid(gridFileName, gridFactory); + std::unique_ptr<Grid> grid(gridFactory.create()); + + ProblemStat<Basis> prob("surface", *grid); + prob.initialize(INIT_ALL); + + auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0); + prob.addMatrixOperator(opL, 0, 0); + + double kappa = Parameters::get<double>("surface->kappa").value_or(1.0); + auto opM = makeOperator(tag::test_trial{}, -kappa*kappa); + prob.addMatrixOperator(opM, 0, 0); + + auto opForce = makeOperator(tag::test{}, X(0) + X(1) + X(2)); + prob.addVectorOperator(opForce, 0); + + AdaptInfo adaptInfo("adapt"); + prob.assemble(adaptInfo); + prob.solve(adaptInfo); + + prob.writeFiles(adaptInfo); + + AMDiS::finalize(); + return 0; +}