Geometry models

Geometry models#

The geometry model is responsible for describing the domain in which we want to solve the equations. A domain is described in deal.II by a coarse mesh and, if necessary, an object that characterizes the boundary. Together, these two suffice to reconstruct any domain by adaptively refining the coarse mesh and placing new nodes generated by refining cells onto the surface described by the boundary object. The geometry model is also responsible for marking different parts of the boundary with different boundary indicators for which one can then, in the input file, select whether these boundaries should be Dirichlet-type (fixed temperature) or Neumann-type (no heat flux) boundaries for the temperature, and what kind of velocity conditions should hold there. In deal.II, a boundary indicator is a number of type types::boundary_id, but since boundaries are hard to remember and get right in input files, geometry models also have a function that provide a map from symbolic names that can be used to describe pieces of the boundary to the corresponding boundary indicators. For example, the simple box geometry model in 2d provides the map {"left"\rightarrow0, "right"\rightarrow1, "bottom"\rightarrow2,"top"\rightarrow3}, and we have consistently used these symbolic names in the input files used in this manual.

To implement a new geometry model, you need to overload the aspect::GeometryModel::Interface class and use the ASPECT_REGISTER_GEOMETRY_MODEL macro to register your new class. The implementation of the new class should be in namespace aspect::GeometryModel.

Specifically, your new class needs to implement the following basic interface:

template <int dim>
    class aspect::GeometryModel::Interface
    {
      public:
        virtual
        void
        create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const = 0;

        virtual
        double
        length_scale () const = 0;

        virtual
        double depth(const Point<dim> &position) const = 0;

        virtual
        Point<dim> representative_point(const double depth) const = 0;

        virtual
        double maximal_depth() const = 0;

        virtual
        std::set<types::boundary_id_t>
        get_used_boundary_indicators () const = 0;

        virtual
        std::map<std::string,types::boundary_id>
        get_symbolic_boundary_names_map () const;

        static
        void
        declare_parameters (ParameterHandler &prm);

        virtual
        void
        parse_parameters (ParameterHandler &prm);
    };

The kind of information these functions need to provide is extensively discussed in the documentation of this interface class at aspect::GeometryModel::Interface. The purpose of the last two functions has been discussed in the general overview of plugins above.

The create_coarse_mesh function does not only create the actual mesh (i.e., the locations of the vertices of the coarse mesh and how they connect to cells) but it must also set the boundary indicators for all parts of the boundary of the mesh. The deal.II glossary describes the purpose of boundary indicators as follows:

In a Triangulation object, every part of the boundary is associated with a unique number (of type types::boundary_id) that is used to identify which boundary geometry object is responsible to generate new points when the mesh is refined. By convention, this boundary indicator is also often used to determine what kinds of boundary conditions are to be applied to a particular part of a boundary. The boundary is composed of the faces of the cells and, in 3d, the edges of these faces.

By default, all boundary indicators of a mesh are zero, unless you are reading from a mesh file that specifically sets them to something different, or unless you use one of the mesh generation functions in namespace GridGenerator that have a ‘colorize’ option. A typical piece of code that sets the boundary indicator on part of the boundary to something else would look like this, here setting the boundary indicator to 42 for all faces located at \(x=-1\):

for (typename Triangulation<dim>::active_cell_iterator
         cell = triangulation.begin_active();
       cell != triangulation.end();
       ++cell)
    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
      if (cell->face(f)->at_boundary())
        if (cell->face(f)->center(true)[0] == -1)
          cell->face(f)->set_boundary_indicator (42);

This calls functions TriaAccessor::set_boundary_indicator. In 3d, it may also be appropriate to call TriaAccessor::set_all_boundary_indicators instead on each of the selected faces. To query the boundary indicator of a particular face or edge, use TriaAccessor::boundary_indicator.

The code above only sets the boundary indicators of a particular part of the boundary, but it does not by itself change the way the Triangulation class treats this boundary for the purposes of mesh refinement. For this, you need to call Triangulation::set_boundary to associate a boundary object with a particular boundary indicator. This allows the Triangulation object to use a different method of finding new points on faces and edges to be refined; the default is to use a StraightBoundary object for all faces and edges. The results section of step-49 has a worked example that shows all of this in action.

The second use of boundary indicators is to describe not only which geometry object to use on a particular boundary but to select a part of the boundary for particular boundary conditions. […]

Note: Boundary indicators are inherited from mother faces and edges to their children upon mesh refinement. Some more information about boundary indicators is also presented in a section of the documentation of the Triangulation class.

Two comments are in order here. First, if a coarse triangulation’s faces already accurately represent where you want to pose which boundary condition (for example to set temperature values or determine which are no-flow and which are tangential flow boundary conditions), then it is sufficient to set these boundary indicators only once at the beginning of the program since they will be inherited upon mesh refinement to the child faces. Here, at the beginning of the program is equivalent to inside the create_coarse_mesh() function of the geometry module shown above that generates the coarse mesh.

Secondly, however, if you can only accurately determine which boundary indicator should hold where on a refined mesh – for example because the coarse mesh is the cube \([0,L]^3\) and you want to have a fixed velocity boundary describing an extending slab only for those faces for which \(z>L-L_{\text{slab}}\) – then you need a way to set the boundary indicator for all boundary faces either to the value representing the slab or the fluid underneath after every mesh refinement step. By doing so, child faces can obtain boundary indicators different from that of their parents. deal.II triangulations support this kind of operations using a so-called post-refinement signal. In essence, what this means is that you can provide a function that will be called by the triangulation immediately after every mesh refinement step.

The way to do this is by writing a function that sets boundary indicators and that will be called by the Triangulation class. The triangulation does not provide a pointer to itself to the function being called, nor any other information, so the trick is to get this information into the function. C++ provides a nice mechanism for this that is best explained using an example:

#include <deal.II/base/std_cxx1x/bind.h>

    template <int dim>
    void set_boundary_indicators (parallel::distributed::Triangulation<dim> &triangulation)
    {
      ... set boundary indicators on the triangulation object ...
    }

    template <int dim>
    void
    MyGeometry<dim>::
    create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const
    {
      ... create the coarse mesh ...

      coarse_grid.signals.post_refinement.connect
        (std_cxx1x::bind (&set_boundary_indicators<dim>,
                          std_cxx1x::ref(coarse_grid)));

    }

What the call to std_cxx1x::bind does is to produce an object that can be called like a function with no arguments. It does so by taking the address of a function that does, in fact, take an argument but permanently fix this one argument to a reference to the coarse grid triangulation. After each refinement step, the triangulation will then call the object so created which will in turn call set_boundary_indicators<dim> with the reference to the coarse grid as argument.

This approach can be generalized. In the example above, we have used a global function that will be called. However, sometimes it is necessary that this function is in fact a member function of the class that generates the mesh, for example because it needs to access run-time parameters. This can be achieved as follows: assuming the set_boundary_indicators() function has been declared as a (non-static, but possibly private) member function of the MyGeometry class, then the following will work:

#include <deal.II/base/std_cxx1x/bind.h>

    template <int dim>
    void
    MyGeometry<dim>::
    set_boundary_indicators (parallel::distributed::Triangulation<dim> &triangulation) const
    {
      ... set boundary indicators on the triangulation object ...
    }

    template <int dim>
    void
    MyGeometry<dim>::
    create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const
    {
      ... create the coarse mesh ...

      coarse_grid.signals.post_refinement.connect
        (std_cxx1x::bind (&MyGeometry<dim>::set_boundary_indicators,
                          std_cxx1x::cref(*this),
                          std_cxx1x::ref(coarse_grid)));
    }

Here, like any other member function, set_boundary_indicators implicitly takes a pointer or reference to the object it belongs to as first argument. std::bind again creates an object that can be called like a global function with no arguments, and this object in turn calls set_boundary_indicators with a pointer to the current object and a reference to the triangulation to work on. Note that because the create_coarse_mesh function is declared as const, it is necessary that the set_boundary_indicators function is also declared const.

Note

For reasons that have to do with the way the parallel::distributed::Triangulation is implemented, functions that have been attached to the post-refinement signal of the triangulation are called more than once, sometimes several times, every time the triangulation is actually refined.