Section author: Andreas Dedner <a.s.dedner@warwick.ac.uk>, Robert Klöfkorn <robert.kloefkorn@iris.no>, Martin Nolte <nolte.mrtn@gmail.com>

Grid Views: Adaptivity and Moving Domains

When constructing a grid, the object returned to Python is always the so called LeafGridView. Without any refinement this is simply a view on all the elements of the grid. As soon as the grid is refined the leaf grid view changes so that it always contains the leaf elements the grid, i.e., the elements on the finest level. Since it is a read only view refinement is carried out using the underlying hierarchical grid, i.e.,

grid.hierarchicalGrid.globalRefine(1)

For a given hierarchical grid one can use different views, i.e., a view on all the elements of a given level:

from dune.grid import cartesianDomain
from dune.alugrid import aluConformGrid, aluSimplexGrid
domain = cartesianDomain([0,0],[1,1],[1,1])
# first construct a grid using quartering as refinement strategy
aluView = aluSimplexGrid(domain)
hGrid = aluView.hierarchicalGrid
hGrid.globalRefine(5)
for level in range(hGrid.maxLevel):
    print("level:",level, "number of elements:",hGrid.levelView(level).size(0))
Created parallel ALUGrid<2,2,simplex,nonconforming> from input stream.
WARNING (ignored): Could not open file 'alugrid.cfg', using default values 0 < [balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.

You are using DUNE-ALUGrid, please don't forget to cite the paper:
Alkaemper, Dedner, Kloefkorn, Nolte. The DUNE-ALUGrid Module, 2016.
level: 0 number of elements: 2
level: 1 number of elements: 8
level: 2 number of elements: 32
level: 3 number of elements: 128
level: 4 number of elements: 512
# now construct a grid using bisection as refinement strategy
aluView = aluConformGrid(domain)
hGrid = aluView.hierarchicalGrid
hGrid.globalRefine(10)
for level in range(hGrid.maxLevel):
    print("level:",level, "number of elements:",hGrid.levelView(level).size(0))
GridParameterBlock: Parameter 'refinementedge' not specified, defaulting to 'ARBITRARY'.
WARNING (ignored): Could not open file 'alugrid.cfg', using default values 0 < [balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.

You are using DUNE-ALUGrid, please don't forget to cite the paper:
Alkaemper, Dedner, Kloefkorn, Nolte. The DUNE-ALUGrid Module, 2016.
Created parallel ALUGrid<2,2,simplex,conforming> from input stream.

level: 0 number of elements: 2
level: 1 number of elements: 4
level: 2 number of elements: 8
level: 3 number of elements: 16
level: 4 number of elements: 32
level: 5 number of elements: 64
level: 6 number of elements: 128
level: 7 number of elements: 256
level: 8 number of elements: 512
level: 9 number of elements: 1024

DUNE-FEM provides a number of additional views which will be discussed further in this chapter:

  • dune.fem.view.adaptiveLeafGridView: this view should be used when the grid is supposed to be locally adapted. The view is still on the leaf elements of the grid but way data is attached to the entities of the grid is optimized for frequent grid changes as caused by local adaptivity. Its usage is shown in the following example.

  • dune.fem.view.geometryGridView: this is an example of a meta grid view which is constructed from an existing grid view and replaces some aspect - in this case the geometry of each element using a given grid function. This concept makes it easy to perform simulations for example on complex domains or on moving grids as shown in Evolving Domains.

  • dune.fem.view.filteredGridView: allows the user to select a subset of a given grid view by providing a filter on the elements. In its simplest version only the iterator is replaced with an iterator over the elements in the filter but in addition it is also possible to obtain a new index set with indices restricted to the elements in the filter.

FilteredGridView examples

In this example we use filteredGridView to filter out parts of a given grid view.

from ufl import SpatialCoordinate, dot
from dune.grid import cartesianDomain
from dune.alugrid import aluConformGrid as leafGridView
from dune.fem.view import filteredGridView
from dune.fem.space import lagrange

Create a host grid view of the underlying grid as usual.

gridView = leafGridView( cartesianDomain([0,0],[1,1],[16,16]) )
Created parallel ALUGrid<2,2,simplex,conforming> from input stream.
GridParameterBlock: Parameter 'refinementedge' not specified, defaulting to 'ARBITRARY'.
WARNING (ignored): Could not open file 'alugrid.cfg', using default values 0 < [balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.

You are using DUNE-ALUGrid, please don't forget to cite the paper:
Alkaemper, Dedner, Kloefkorn, Nolte. The DUNE-ALUGrid Module, 2016.

Now create a filteredGridView with a simple callable specifying which element belongs to the new domain.

filteredView = filteredGridView(gridView, lambda e: e.geometry.center.two_norm > 0.5, domainId=1)
space = lagrange(filteredView, order=2)
x = SpatialCoordinate(space)
solution = space.interpolate(dot(x,x),name="solution")
solution.plot()
print("number of dofs:", solution.size,\
      "integral over filtered domain",solution.integrate())
filteredgridview_nb_files/filteredgridview_nb_5_0.png
number of dofs: 1089 integral over filtered domain 0.6413014729817705

Revert the filter and create a filteredGridView with also provides an overloaded index set with a consecutive index for the entities belonging to the new domain.

filteredView = filteredGridView(gridView, lambda e: e.geometry.center.two_norm < 0.5, domainId=1,
                                useFilteredIndexSet=True)
space = lagrange(filteredView, order=2)
x = SpatialCoordinate(space)
solution = space.interpolate(dot(x,x),name="solution")
solution.plot()
print("number of dofs:", solution.size,\
      "integral over filtered domain",solution.integrate())
filteredgridview_nb_files/filteredgridview_nb_7_0.png
number of dofs: 239 integral over filtered domain 0.025365193684895825

Compare with the original grid view.

space = lagrange(gridView, order=2)
solution = space.interpolate(dot(x,x),name="solution")
solution.plot()
print("number of dofs:", solution.size,\
      "integral over filtered domain",solution.integrate())
filteredgridview_nb_files/filteredgridview_nb_9_0.png
number of dofs: 1089 integral over filtered domain 0.6666666666666667

Dynamic Local Grid Refinement and Coarsening

For refining and coarsening a grid locally the dune.fem module provides a function adapt. The storage of all discrete functions will be automatically resized to accommodate the changes in the grid but the resulting dof vector will not be initialized. To prolong and restrict data from the old to the new grid, the corresponding discrete functions have to be passed to the dune.fem.adapt method:

fem.adapt(u1,u2,...,uN)

The module dune.fem also provides a globalRefine(level,*dfs) method, where a negative level globally coarsens the grid. If discrete functions are passed in they will be prolong (restricted), the dof vectors of all other dof vectors will be resized.

Note

if the underlying storage of a discrete function is stored on the Python side as a numpy array, i.e., vec = uh.as_numpy was called, then access to vec will be undefined after a grid modification since the underlying buffer change will not have been registered.

The module dune.fem provides a function for marking elements for refinement/coarsening:

def mark(indicator, refineTolerance, coarsenTolerance=0,
    minLevel=0, maxLevel=None):

where indicator is a grid function. An element \(T\) is marked for refinement if the value of indicator on \(T\) is greater then refineTolerance and coarsened if the value is less then coarsenTolerance. The element \(T\) is not refined if its level is already at maxLevel and not coarsened if its level it at minLevel. This method can for example be used to refine the grid according to an equal distribution strategy by invoking

dune.fem.mark(indicator, theta/grid.size(0))

where theta is a given tolerance.

A layered Doerfler strategy is also available

def doerflerMark(indicator, theta, maxLevel=None, layered=0.05):

The following two examples showcase adaptivity: the first one using a residual a-posteriori estimator for an elliptic problem, the second one shows adaptivity for a time dependent phase field model for crystal growth. At the end of this section a dual weighted residual approach is used to optimize the grid with respect to the error at a given point. While the first two examples can be implemented completely using the available Python bindings the final example requires using a small C++ snippet which is easy to integrate into the Python code.

Evolving Domains

As mentioned above DUNE-FEM provides a grid view that makes it easy to exchange the geometry of each entity in the grid. To setup such a grid view one first needs to construct a standard grid view, i.e., a leafGridView and define a grid function over this view using for example a discrete function, a UFL function, or one of the concepts described in the section Grid Function. Note that the topology of the grid does not change, i.e., how entities are connected with each other. The following shows an example of how to change a grid of the unit square into a grid of a diamond shape:

# <markdowncell>
import matplotlib
matplotlib.rc( 'image', cmap='jet' )
from matplotlib import pyplot
from ufl import sqrt, SpatialCoordinate, triangle, as_vector
from dune.grid import structuredGrid
from dune.fem.view import geometryGridView
from dune.fem.function import uflFunction
square = structuredGrid([0,0],[1,1],[10,10])
x = SpatialCoordinate(triangle)
transform = as_vector([ (x[0]+x[1])/sqrt(2), (-x[0]+x[1])*sqrt(2) ])
gridFunction = uflFunction(square,ufl=transform,order=1,name="diamond")
diamond = geometryGridView(gridFunction)
# <markdowncell>
diamond.plot()
_images/geoview_nb_2_0.png

By using a discrete function to construct a geometry grid view, it becomes possible to simulate problems on evolving domains where the evolution is itself the solution of the partial differential equation. We demonstrate this approach based on the example of surface mean curvature flow first in its simplest setting and then with the evolution of the surface depending on values of a computed surface quantity satisfying a heat equation on the surface:

Todo

add a coupled surface diffusion/evolution problem

https://zenodo.org/badge/DOI/10.5281/zenodo.3706994.svg