Merge branch 'example/surface_grid' into 'master'

Example/surface grid

See merge request !30

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

examples/FactoryParametrization.hpp

+#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

examples/SphereMapping.hpp

+#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

examples/init/surface.dat.2d

+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

examples/surface.cc

+#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;
+}