Dune-CurvedGrid -- A Dune module for surface parametrization
DDune-CurvedGrid – A Dune module for surfaceparametrization
Simon Praetorius and Florian Stenger Technische Universit¨at Dresden, Institut f¨ur Wissenschaftliches Rechnen,Zellescher Weg 12–14, 01069 Dresden, Germany( [email protected] ). Technische Universit¨at Dresden, Institut f¨ur Wissenschaftliches Rechnen,Zellescher Weg 12–14, 01069 Dresden, Germany( [email protected] ).October 8, 2020
Abstract
In this paper we introduce and describe an implementation of curved (surface) geometrieswithin the
Dune framework for grid-based discretizations. Therefore, we employ the abstrac-tion of geometries as local-functions bound to a grid element, and the abstraction of a grid asconnectivity of elements together with a grid-function that can be localized to the elementsto provide element local parametrizations of the curved surface.
Numerical computations on curved surfaces are an important tool for studying physical phenom-ena in thin structures, on curved boundaries of domains, and on interfaces between two bulkregions. Those problem arise, e.g., in fluid dynamics when considering free-surface flows, in inter-facial transport problems on biological membranes or fluid droplets, in mathematical geosciences,and in the physics of thin films with vanishing film thickness. All these applications require thediscretization of partial differential equations (PDEs) on the embedded surface, involving quan-tities like curvature, surface measures, normal vectors, and covariant derivatives that need to beavailable by a numerical method. An overview about some of these applications can be found inBothe and Reusken (2017); Nestler et al. (2017); Jankuhn et al. (2018); Freeden and Schreiner(2009).In discretization methods based on finite elements or finite volumes one is often faced with theproblem of representing the curved geometry by numerical grids and discretizing a (geometric)partial differential equation on this approximation of the actual smooth surface. One distinguishesbetween implicit surface representations using cut cells, level-sets or diffuse interfaces, see R¨atzand Voigt (2006); Burman et al. (2018); Olshanskii and Reusken (2017), and explicit surfacerepresentations by triangulations and surface finite element or surface finite volume schemes.For the latter, the representation with piece-wise flat elements is a widely used lowest orderapproximation that is easy to construct and to implement with standard software. Unfortunately,it has some drawbacks. The local flat approximation leads to vanishing curvature inside theelements and this can sometimes lead to non-converging numerical schemes, as shown in Heine(2004) for the discrete mean-curvature vector and Weingarten map, in Fritz (2013) for a finite-element approximation of the Ricci curvature, and in Hansbo et al. (2019) for the discretization1 a r X i v : . [ c s . M S ] O c t f a surface vector Laplacian . The numerical error involved in a discretization of PDEs oncurved surfaces depends on two properties of the discretization, the representation of the objectivefunction and the representation of the geometry. Thus, a higher-order scheme is only possible withalso a higher-order description of the surface approximation.Thus, better approximations of the geometry than pice-wise flat are necessary. Those canbe found in piece-wise polynomial approximations or even exact representations of the surfaceif an analytical description of the surface’s geometry is available. Numerical libraries for theimplementation of partial-differential equation solvers do not always provide a usable interface forsuch non-linear geometry transformation or representation of such curved grids and often do notprovide the utilities necessary to implement such simulations efficiently.The present work tries to fill this gap for the numerical library Dune , Bastian et al. (2020),a framework for the discretization of grid-based numerical problems. This modular library iscentered around the abstraction of a grid interface and the coupling to various different gridimplementations for structured or unstructured grids, with or without adaptivity, sequential orparallel traversal, for volume and surface grids and additionally allowing to wrap any grid imple-mentation to extend or modify its functionality. This concept of wrapped grids, called meta-grids in Dune , is the basis of our implementation. We provide a wrapper around any grid implementa-tion in
Dune , transforming a (piece-wise) flat reference grid into a curved grid using a (non-linear)geometry transformation.All implementations discussed in this paper can be found in individual
Dune modules inpublicly available source code repositories: Praetorius and Stenger (2020a); Praetorius (2020,2019); Praetorius and Stenger (2020c). The code examples and results in this manuscript aresummarized and made available in the repository Praetorius and Stenger (2020b).
For a motivation of the introduced functionality in
Dune-CurvedGrid , we consider the polyno-mial approximation of a spherical surface by local polynomials of some given order . C++ code // 1. Construct a reference grid std :: unique_ptr refGrid = Gmsh4Reader < FoamGrid <2,3> >:: createGridFromFile (" sphere . msh "); // 2. Define the geometry mapping auto sphere = []( const auto & x) { return x / x. two_norm (); }; auto sphereGridFct = analyticDiscreteFunction ( sphere , * refGrid , order ); // 3. Wrap the reference grid to build a curved grid CurvedGrid grid {* refGrid , sphereGridFct };
The curved grid is build on top of a piece-wise flat reference grid provided by a
Gmsh file andrepresented by a Dune-FoamGrid , Sander et al. (2017). Although this reference grid is createdusing the
Gmsh4Reader which is introduced in section 7, the grid could also be created by any othermethod . The polynomial representation is given by interpolating an analytic coordinate projec-tion, sphere , i.e., a simple normalization of global coordinates, locally into a polynomial space ofgiven polynomial order . This is achieved by the function wrapper AnalyticDiscreteFunction ,that provides evaluation of values and derivatives of the given function in terms of its local discreteapproximation, see section 5.The library provides grid-functions that can be used to represent geometry mappings. Both,analytical and discrete representations, are implemented and thus even evolving parametrizationsas solutions of PDEs or other external descriptions are possible. In these three examples a stabilization technique or additional higher-order geometric knowledge allows toovercome these flat element limitations The
Gmsh format is described in Geuzaine and Remacle (2009) For example
Dune-Grid provides it’s own
GmshReader which supports the older file format version 2.
CurvedGrid fulfills the
Dune grid requirements and can be used instead of any other regular
Dune grid. It provides the geometric mappings for grid elements and for element intersections,allowing for continuous and discontinuous discretization schemes build on top of this grid.
In section 2 the mathematical foundation is laid, describing smooth surfaces and their (polyno-mial) approximation, the grid and its geometry transforms as well as grid-functions. This sectionintroduces the notation and defines the objects and mappings that are implemented in the library.The subsequent section 3 introduces the interface for the geometry classes representing the localmapping of coordinates. This is followed by a section about the class interface for the
CurvedGrid in section 4 that implements the grid wrapper providing locally defined curved geometries. Thedescription of the actual surface or its approximation requires grid-functions, introduced in sec-tion 5, and projection mappings that are shown by examples in section 6. File-readers and writersfor curved geometries are introduced in section 7. This concludes the implementation aspects ofthe library. In section 8 numerical validation of the implementation is given by analyzing knownerror bounds for geometric quantities like distance between surfaces, normal vectors, and meancurvature. This is followed by a selection of numerical examples of finite-element problems oncurved surfaces.
Let Γ ⊂ R m +1 be an oriented, connected, and smooth m –dimensional manifold. Γ can be describedin multiple ways, e.g., by a parametrization over a reference domain, by an implicit representationas level-set of a function, or by closest-point projection of coordinates on another manifold ina close neighborhood of Γ. All these descriptions have advantages and disadvantages and aresummarized in Dziuk and Elliott (2013). While continuous descriptions allow to extract geometricmeasures and characteristics of the surface, like its metric or curvature, they are complicatedto use in numerical computations. Hence, an approximation, or piece-wise representation of thesurface for local evaluation of quantities and data is desirable. Such a representation might be given by a piece-wise flat surface Γ h , topologically equivalent to thesmooth surface Γ. This reference surface is composed of finitely many regular and quasi-uniform(flat) m –dimensional elements of diameter h . The collection of these patches, typically simplicesor hyper-cubes, is denoted by G h and is called the grid representation of Γ h , withΓ h = (cid:91) e ∈G h e (1)where e denotes an element of the grid. We assume that the patches do not overlap, i.e., for e , e ∈ G h we have that int( e ) ∩ int( e ) = ∅ and if e ∩ e = I (cid:54) = ∅ and dim( I ) = m −
1, it iscalled an intersection of e and e and is assumed to be a subset of an ( m − e and e , respectively.Each element of the grid e ∈ G h is parametrized over a reference element ˆ e ⊂ R m by aninvertible and differentiable mapping µ e : ˆ e → e , called the geometry mapping of e . Additionally,we assume that there exists a bijective mapping X : Γ h → Γ, such that the smooth surface can berepresented by the union of (non-overlapping) mapped patches, i.e.,Γ = (cid:91) e ∈G h X ( e ) = (cid:91) e ∈G h X ( µ e (ˆ e )) =: (cid:91) e ∈G h X e (ˆ e ) . (2)With this property, we call Γ h the reference manifold or reference domain of Γ and the family { X e } e ∈G h its reference parametrization. 3 .2 (Higher-order) Approximations of the Manifold The reference manifold Γ h from the last section is not used directly for an approximate discretiza-tion of functions on Γ, since it does not necessarily approximate the smooth manifold well enough.It just provides a reference domain for the parametrization X . For numerical computations anddiscretizations, we need another manifold in the proximate neighborhood of Γ. For a pice-wisepolynomial surface approximation, we follow the general notation of Demlow (2009).Therefore, let X k := I kh X ∈ P k ( e ) be a k th-order polynomial (Lagrange) interpolation of themapping X on the element e of the reference manifold, with Lagrange nodes sitting on the smoothsurface Γ. P k ( e ) denotes the space of polynomials on e of degree at most k and I the (component-wise) Lagrange interpolation operator. The interpolation can be expressed in terms of local basisfunctions in the reference element ˆ e , by X k ( µ e (ˆ x )) := X ke (ˆ x ) = n k (cid:88) j =1 ξ j · φ j (ˆ x ) , for ˆ x ∈ ˆ e , (3)with ξ j := X e (ˆ x j ) ∈ Γ for { ˆ x j } j =1 ...n k local Lagrange nodes, the corresponding local Lagrangebasis functions { φ j } j =1 ...n k , and n k the number of local basis functions of P k (ˆ e ). The nodes andbasis functions fulfill the nodal interpolation property φ i (ˆ x j ) = δ ij . Thus, the mapped Lagrangenodes ξ j sit on the smooth surface Γ, see Figure 1 for an illustration.Then, the k th-order approximation Γ kh of the smooth surface Γ can be obtained by the unionof (non-overlapping) mapped elements, mapped by X k :Γ kh = (cid:91) e ∈G h X k ( e ) = (cid:91) e ∈G h X k ( µ e (ˆ e )) = (cid:91) e ∈G h X ke (ˆ e ) . (4)In case of k = 1 we speak of a piece-wise flat or polyhedral surface grid Γ h . Since the mappings X ke are locally smooth, we obtain a piece-wise differentiable manifold.ˆ e X e (ˆ e )ˆ x ˆ x ˆ x ˆ x ˆ x ˆ x ξ ξ ξ ξ ξ ξ Figure 1: Lagrange parametrization of order k = 2 with Lagrange nodes on vertices and edges. The reference manifold Γ h is composed of elements e that build the grid G h , or vice versa, thegrid defines the manifold. The (higher-order) mapped elements and the elements mapped to thesmooth manifold also form grids, namely G kh := { X k ( e ) } e ∈G h and G := { X ( e ) } e ∈G h . (5)Since, we assume that all manifolds Γ, Γ h , and Γ kh have the same dimension, also the grids arecomposed of elements e of that same dimension m . We speak of a conforming grid if the surfaceis continuous and the non-empty intersection of each two elements e and e is an l –dimensional4acet of both elements, called sub-entity s (cid:22) e and s (cid:22) e . We say s has co-dimension c = m − l .In conforming grids the sub-entities of co-dimension one are the intersections I of two elements.Corresponding to the reference element ˆ e there is a reference element ˆ s ⊂ R l of the sub-entity s of e . The relation between the geometries of ˆ s and ˆ e is given by the invertible anddifferentiable mapping η s,e : ˆ s → ˆ e , called the local-geometry mapping between sub-entity andelement. With this, the parametrization of the real sub-entity s is given by the chain of η and µ ,i.e., s = µ e ( η s,e (ˆ s )), see Figure 2 for an illustration. If this parametrization is equivalent to thedirect mapping from the reference element ˆ s to s , i.e., if it holds µ e ◦ η s,e = µ s for all sub-entities s of e , we call the grid twist-free , see Dedner and Nolte (2012). This property is assumed tohold at least for intersections, in the following. Examples of twist-free grids are Dune::OneDGrid , Dune::YaspGrid , and
Dune::ALUGrid .ˆ s ˆ e e X ( e ) s X ( s )Γ h Γ µ e XX e η s,e µ s Figure 2: Coordinate mappings µ and η between reference element, flat element, and curvedelement. An additional mapping X s,e can be defined by chaining of η, µ and X , i.e. X s,e = X ◦ µ e ◦ η s,e . When discussing the mapping X or its polynomial variants X k , we often define it by its localrepresentation X e and X ke , respectively, with e the element it is defined on. Those functions,defined via their local element variant, are called grid-functions in the following and are directlyconnected to their local variant, called local-functions . The evaluation of the grid-function in globalcoordinates x ∈ R n might be an expensive operation, whereas the evaluation in the correspondinglocal coordinate ˆ x can be easily defined. An example is the evaluation of a discrete function bylinear combination of evaluations of local basis functions inside the grid elements.In general, we denote by f e : ˆ e → R , f e = f ◦ µ e the local-function bound to the element e . Itis associated to the grid-function f : e → R with range R . If f is smooth or at least differentiableinside the element e , we denote by Df : e → L ( e, R ) its derivative as linear mapping, oftenrepresented as a matrix R dim( R ) × n . The corresponding localized derivative ( Df ) e is then given by( Df ) e : ˆ e → L ( e, R ) , ( Df ) e = Df ◦ µ e , (6)with the same linear mapping in the range as the global Df . This notation follows Engwer et al.(2017). If instead just a local jacobian D ( f e ) : ˆ e → L (ˆ e, R ) is available or requested, we have therelation D ( f e ) = ( Df ) e · D ( µ e ).Note that the geometry mapping X is a grid-function with local-function X e . Also, theparametrized geometry mappings X k and X ke are grid-function and local-function, respectively.5hese mappings are differentiable, by differentiating their local basis functions φ j , i.e., D ( X ke ) = n k (cid:88) j ξ j ⊗ D ( φ j ) . (7)In the next section we want to introduce the implementation of the element geometry mappings X e and X ke and then in the subsequent section a wrapper to transform the reference grid G h into G or G kh using the geometry mappings and its global variants X and X k from above. A Geometry is a mapping from local coordinates in R l to global coordinates in R n , where the localcoordinates are in the coordinate system of an entity e which this geometry belongs to. The entitycan be an element of the grid, or a sub-entity of any co-dimension. This flexibility requires thegeometry parametrization to be evaluable in different coordinate systems and also its derivativesto be available for the corresponding coordinate transformations.Depending on how the parametrization of the geometry is given, different implementations areprovided. The ParametrizedGeometry expects the mapping X e and a local finite-element andconstructs the local interpolation X ke internally, whereas the LocalFunctionGeometry wraps agiven X e or X ke directly. The first implementation,
ParametrizedGeometry , expects only a callable function that mapsentity local coordinates to global coordinates. This mapping is internally interpolated into a localfinite-element space, e.g., local Lagrange functions, that allows to evaluate values and Jacobians ofthe parametrization from linear combinations of evaluated local basis functions and its derivatives.
C++ code // < dune / curvedgeometry / parametrizedgeometry .hh > template < class LocalFiniteElement , int coorddim , class Traits = (...) > class ParametrizedGeometry ; The template parameters are defined by
LocalFiniteElement
A class representing a local finite-element in the sense of
Dune-LocalFunctions . coorddim Dimension n ≥ l of the global coordinates this geometry maps into. Traits (optional) Parameters for internal optimization of operations.This geometry mapping directly corresponds to the description of the discrete geometry X k where the input function represents the mapping X combined with the local-to-global mapping µ e ◦ η s,e . Thus, the actual input is X s,e . In the geometry, the local interpolation X ks,e is representedby the interpolation coefficients { X s,e (ˆ x j ) } j , i.e., Lagrange nodes on the surface Γ, and the setof local basis functions { φ j } associated to these nodes, see (3). The jacobian of the geometrymapping can thus be provided by evaluating the gradients of the local basis functions and itslinear combination with the stored coefficients, see (7).The LocalFiniteElement parameter here is the crucial input characterizing which type oflocal basis functions and local interpolation to use for calculating and representing the (Lagrange)nodes.
Dune-LocalFunctions provides various implementations of local finite-elements, likeLagrange functions on all supported geometry types. The corresponding local finite-element canbe obtained either by an explicit instantiation if the geometry type is known and identical for allelements, or by using a local finite-element cache. The latter provides the local finite-element ofone kind for all geometry types by type-erasure or variadic visitors, see the example below.6 .2 Geometry with Differentiable Parametrization
The second implementation,
LocalFunctionGeometry , expects a mapping for coordinates andadditionally the jacobian of that mapping so that the geometry Jacobian can be represented di-rectly by the function. We expect this mapping function to be compatible with a
LocalFunction interface, to be defined below. A local-function f e typically can be evaluated only in elementlocal coordinates of the element e , denoted by R m , but not in codim > R l . In order to allow the geometry to be defined also for these entities, or even elementintersections, this geometry implementation is parametrized additionally with a LocalGeometry coordinate transform η . This coordinate transform maps the entity-local coordinates to the grid-element local coordinates where the local-function can be evaluated in. Thus, the geometry map-ping is a chaining R l → R m → R n with l ≤ m ≤ n , that is, X = f e ◦ η = f ◦ µ e ◦ η , where f is theglobal grid-function associated to the local-function f e that is bound to an element e . C++ code // < dune / curvedgeometry / localfunctiongeometry .hh > template < class LocalFunction , class LocalGeometry , class Traits = (...) > class LocalFunctionGeometry ; template < class LocalFunction , class ctype , int dim > using ElementLocalFunctionGeometry = LocalFunctionGeometry < LocalFunction , DefaultLocalGeometry < ctype ,dim , dim >, LocalFunctionGeometryTraits < ctype > >;
This geometry is parametrized with the types
LocalFunction , LocalGeometry , and
Traits that fulfill the following requirements:
LocalGeometry represents a geometric mapping from an entity of codim c to the element with codim
0. The geometry is bound to the domain element whereas the
LocalFunction can bebound to the range element of this geometric mapping. Thus, it is a differentiable function η : R l → R m with l = m − c that fulfills a reduced Dune::Geometry concept, i.e., thereis a local-to-global mapping from the coordinate system of the codim- c entity to the elementgeometry. Let η be of type LocalGeometry and ˆ x of type LocalCoordinate , then the expression η (ˆ x ) ∈ R m results in a type that is the domain type of the Localfunction .The derivative of this parametrization must be accessible by evaluating the geometry method η .jacobianTransposed( ˆ x ) that returns the transposed of the jacobian of η , that is, D (cid:62) η with Dη : R l → L ( R n , R l ) ∼ = R n × l . LocalFunction represents a differentiable mapping f e : R m → R n with given derivative ( Df ) e : R m → L ( R n , R n ) ∼ = R n × n . It is required that f e is a model of the concept Callable , i.e., let x be oftype LocalGeometry::GlobalCoordinate , then the expression f e ( x ) must result in a valid typedenoted by the GlobalCoordinate of the
LocalFunctionGeometry . f e must be differentiable ,i.e., there must exist a function derivative( f ) whose return type is another model of the LocalFunction concept. It returns the global derivative D ( f ) e of the grid-function f associatedto f e . Traits (optional) is a class holding parameters for the implementation of the geometry, like thetolerance and iteration limit for a Newton solver implementing the global-to-local function. Addi-tionally it allows to specify element properties that cannot be deduced from the
LocalFunction or LocalGeometry directly, like the
GeometryType of the element, if there is only one.The definition of the
LocalFunction follows the definition of localized functions in Engweret al. (2017). Especially, the definition of the range of the derivative of the local-function asderivative w.r.t. global coordinates is taken from there. The final jacobianTransposed of thegeometry mapping X is thus given by the chaining D (cid:62) X = D (cid:62) η · ( D (cid:62) µ e ◦ η ) · ( Df ◦ µ e ◦ η ) (cid:62) . Notethat the mappings X , η , and µ e are geometry mappings and thus provide only the transposed oftheir Jacobians, whereas the localized function f e provides the non-transposed Jacobian derivative.In order to provide the transposed of the final geometry Jacobian D X , it needs to be a linear map7hat is transposible , i.e. it can be applied in a transposed fashion to a vector, by implementing themethod mtv() or it must be representable as a matrix. The two geometry implementations can now be used directly to parametrize a surface while travers-ing a flat reference grid G h . The following examples show the wrapping of a flat geometry into a LocalFunctionGeometry and
ParametrizedGeometry , respectively.At first we have to define a reference grid that provides the actual elements e and the mapping X for the curved geometry, here constructed using the VtkReader introduced in section 7.
C++ code // Construct a reference grid auto refGrid = VtkReader < FoamGrid <2,3> >:: createGridFromFile (" sphere . vtu "); // Define the geometry mapping auto sphere = SphereProjection <3, double >{1.0}; auto sphereGridFct = analyticGridFunction < FoamGrid <2,3>>( sphere ); The mapping sphere , see section 6, is a differentiable function that can be transformed intoa differentiable grid-function sphereGridFct using the wrapper
AnalyticGridFunction , see sec-tion 5, that is provided by the library.For the construction of a
LocalFunctionGeometry , we have to provide a local-function of thatgrid-function and a
LocalGeometry . In case of wrapping the grid element geometry, i.e., codim is zero, this local-geometry mapping is the identity. An efficient implementation is given by theclass
DefaultLocalGeometry . C++ code // Define a local - function from the grid - function auto sphereLocalFct = localFunction ( sphereGridFct ); auto localGeometry = DefaultLocalGeometry < double ,2,2> {}; // traverse the reference grid for ( const auto & e : elements ( refGrid -> leafGridView ())) { // bind the local - function to the grid element sphereLocalFct . bind (e); // construct the LocalFunctionGeometry from LocalFunction and LocalGeometry LocalFunctionGeometry localFctGeometry (e. type () , sphereLocalFct ,localGeometry ); // ( optionally ) unbind from the element , i.e., free memory and unset variables sphereLocalFct . unbind (); } A local-function must be bound to an element (and optionally unbound at the end of usage).The type supports class-template argument deduction and if the geometry is constructed on thegrid element, the
LocalGeometry argument can even be omitted, defaulting to
DefaultLocalGeometry in this case:
C++ code LocalFunctionGeometry localFctGeometry ( referenceElement (e), sphereLocalFct );
Note, it is necessary to pass the element type as
Dune::ReferenceElement , in order to allow thededuction of the element dimension.Similarly, we can construct a Lagrange
ParametrizedGeometry by using the sphere functionfrom above. Therefore, we have to either use the local-function wrapper or have to construct thelocal-to-global mapping from reference element coordinates to element coordinates directly.8 local finite-element can be provided by using a local finite-element cache, or by explicitinstantiation. Both variants are shown in the example below.
C++ code // < dune / localfunctions / lagrange / lagrangelfecache .hh > LagrangeLocalFiniteElementCache < double , double , 2, order > lfeCache ; // traverse the reference grid for ( const auto & e : elements ( refGrid -> leafGridView ())) { // projection from local coordinates auto X_e = [& sphere , geo =e. geometry ()]( const auto & local ) { return sphere ( geo . global ( local )); }; // construct the ParametrizedGeometry from lfe cache ParametrizedGeometry curvedGeometry (e. type () , lfeCache . get (e. type ()), X_e ); // construct the ParametrizedGeometry from local finite - element auto lfe = LagrangeSimplexLocalFiniteElement < double , double , 2, order > {}; ParametrizedGeometry curvedGeometry2 (e. type () , lfe , X_e ); } Instead of wrapping the geometries manually while traversing the flat reference grid, the whole gridcan be wrapped. This allows to return the wrapped geometry directly in a call to e.geometry() instead of the flat element geometry. The library provides such a grid wrapper with the class
CurvedGrid which is an implementation of G kh or G , depending on the element parametrizationprovided.The class signature is given by C++ code // < dune / curvedgrid / grid .hh > template < class GridFunction , bool useInterpolation = false > class CurvedGrid ; // constructor using Lagrange geometry interpolation template < class GridFunction > template < class RefGrid > CurvedGrid < GridFunction , true > :: CurvedGrid ( const RefGrid &, const GridFunction &, int oder ); // constructor using LocalFunctionGeometry template < class GridFunction > template < class RefGrid > CurvedGrid < GridFunction , false > :: CurvedGrid ( const RefGrid &, const GridFunction &); with template parameters GridFunction the type of a grid-function associated with a reference grid
RefGrid
The reference grid the curved grid is based on useInterpolation (optional) if true, uses Lagrange
ParametrizedGeometry , otherwise constructa geometry of type
LocalFunctionGeometry
This class allows to locally construct both the
LocalFunctionGeometry and the
ParametrizedGeometry ,depending on the properties of the grid-function passed to the grid wrapper and the order pa-rameter given to the grid. If the given order is positive, it is assumed that a local interpolation9hould be constructed and thus a
ParametrizedGeometry with Lagrange local finite-element isused as local geometry parametrization. Otherwise, if order ≤ and the grid-function is locallydifferentiable, a LocalFunctionGeometry is used.The corresponding grid geometry is then similar to
C++ code // < dune / curvedgrid / geometry .hh > template < class ct , int dim > using LFE_t = typename LagrangeLFECache
In the following examples we construct both a wrapper using the
ParametrizedGeometry and the
LocalFunctionGeometry .At first, we construct the grid by wrapping a callable representing the geometry mapping X . C++ code // Construct a reference grid auto refGrid = Gmsh4Reader < FoamGrid <2,3> >:: createGridFromFile (" sphere . msh "); // Define the geometry mapping auto sphere = []( const auto & x) { return x / x. two_norm (); }; // Wrap the reference grid to build a curved grid CurvedGrid grid {* refGrid , sphere , order };
In this example, a grid-function is automatically constructed from the callable sphere as aninstance of
AnalyticGridFunction . If class template-argument deduction from
C++17 cannot beused, a generator function curvedSurfaceGrid(*refGrid, sphere) is provided.In the second example we construct a grid-function first which either uses a local interpolationinternally, or is given as a differentiable function as above.
C++ code // Define a discrete grid - function on the reference grid auto gridFct = discreteGridViewFunction <3>( refGrid -> leafGridView () , order ); // Interpolate the parametrization into the grid - function Functions :: interpolate ( gridFct . basis () , gridFct . coefficients () , sphere ); // Wrap the reference grid to build a curved grid CurvedGrid grid {* refGrid , gridFct };
Here, in the example we use a generic discrete function with range type
FieldVector
Dune-Functions global basis.
C++ code power <3>( lagrange ( order ), blockedInterleaved ()) The basis is stored together with a coefficient vector inside the
DiscreteGridViewFunction , seesection 5. Note that this grid-function is associated to a
GridView and not the whole grid, sincethe global basis is bound to a
GridView . Each time the (reference) grid changes, e.g., by localrefinement or parallel load balancing, the grid-function must be updated as well, using10 ++ code gridFct . update ( refGrid -> leafGridView ()); In order to construct the
ParametrizedGeometry or LocalFunctionGeometry and thus the curvedgrid, parametrizations in form of mappings X e or local-functions must be provided. Various grid-functions are implemented in Dune-CurvedGrid to simplify the construction and to act asreference implementations:
AnalyticGridFunction
Implementation of a grid-function that can be bound to any entity inthe grid given by a
Callable , mapping global coordinates of Γ h to a range type. This range typedefines the global coordinates in the curved geometry. If the callable is differentiable so is thegrid-function. It can thus be used in the LocalFunctionGeometry .The
AnalyticGridFunction can be constructed by
C++ code // < dune / curvedgrid / gridfunctions / analyticgridfunction .hh > template < class Grid , class Functor > class AnalyticGridFunction ; template < class Grid , class Functor > auto analyticGridFunction ( Functor && functor ) -> AnalyticGridFunction
Similarly to
AnalyticGridFunction this grid-function is constructedfrom a
Callable , mapping global coordinates to global coordinates, but the mapping is interpo-lated into a Lagrange basis first. Thus, this grid-function does not represent an exact geometrybut an approximation. Moreover, it provides derivatives by differentiating the local basis func-tions, see (7). It can be used to parametrize
LocalFunctionGeometry .The
AnalyticDiscreteFunction can be constructed by
C++ code // < dune / curvedgrid / gridfunctions / analyticdiscretefunction .hh > template < class Grid , class Functor , int order = -1> class AnalyticDiscreteFunction ; template < class Functor , class Grid > auto analyticDiscreteFunction ( Functor && functor , const Grid &, int order ) -> AnalyticDiscreteFunction
Dune-Functions as module dependency.
DiscreteGridViewFunction
This grid-function is restricted to a specific
GridView of the gridand is build by a set of global basis functions and a coefficient vector, both stored inside thisgrid-function. It can be used as
LocalFunctionGeometry , since the basis functions provide aderivative and thus the grid-function is differentiable. Additionally, the coefficient vector, i.e.,the vector of Lagrange nodes on the smooth surface Γ, can be modified and thus evolving gridscan be parametrized easily.The
DiscreteGridViewFunction can be constructed by11 ++ code // < dune / curvedgrid / gridfunctions / discretegridviewfunction .hh > template < class GridView , int components = GridView :: dimensionworld , int ORDER = -1, class T = double > class DiscreteGridViewFunction ; template < int components , int ORDER = -1, class T = double , class GridView > auto discreteGridViewFunction ( const GridView & gridView , int order = ORDER ) -> DiscreteGridViewFunction < GridView , components , ORDER , T>; This grid-function is not as general as the
DiscreteGlobalBasisFunction of Dune-Functions ,i.e., the range type is fixed to
FieldVector
Dune-Functions as moduledependency.
In order to test PDE discretizations or in benchmarks geometry parametrizations for simple shapesmust be available. A common example is the sphere parametrization used in all the examplesabove. But additionally, shapes with less symmetry might be of interest for benchmarks andnumerical validation. We have implemented the geometries of the sphere, ellipsoid and torus assimple parametrizable shapes. Those three geometries can be obtained by
C++ code // < dune / curvedgrid / geometries / sphere .hh > template < int dim , class ctype = double > class SphereProjection ; // sphere radius r template < class Grid , class T> auto sphereGridFunction (T r) { auto sphere = SphereProjection < Grid :: dimensionworld ,T>{r}; return analyticGridFunction < Grid >( sphere ); } for the sphere parametrization. C++ code // < dune / curvedgrid / geometries / ellipsoid .hh > template < class ctype = double > class EllipsoidProjection ; // major axis a , b , and c template < class Grid , class T> auto ellipsoidGridFunction (T a, T b, T c) { auto ellipsoid = EllipsoidProjection
The explicit surface representation is based on a flat grid that approximates the smooth surfacewith higher-resolution than the target grid we want to construct from a reference grid. Thiswould allow to run simulations on a coarse grid, while the surface is only given by a very fine grid.Additionally, it allows to construct higher-order geometries for a surface that is given only withpiece-wise flat geometries.
C++ code // < dune / curvedgrid / geometries / explicitsurface .hh > template < class ctype = double > class ExplicitSurfaceProjection ; // Constructor with grid and an option to activate caching template < class ctype > template < class Grid > ExplicitSurfaceProjection < ctype > :: ExplicitSurfaceProjection ( const Grid & grid , bool cached = true ); The input to construct the
ExplicitSurfaceProjection is a (surface) grid representing thehigh-resolution (piece-wise) flat surface. Internally, the vertices of the grid are stored in a fastsearch tree, a KDTree implementation based on Blanco and Rai (2014) which supports nearestneighbor search. Each time the projection is evaluated for a global coordinate x the closest vertexin the high-resolution surface grid is searched. Afterwards the adjacent elements are considered.For each of them the closest point to x is determined by orthogonal projection. The closest foundpoint is then returned.This approach only works well if the high-resolution surface grid is of sufficient quality, i.e.,no overly acute-angled elements occur. Otherwise the adjacent elements of the closest vertexneed not contain the actual closest point to x . When choosing a method for creating the grid aDelaunay-triangulation, for instance, is a good candidate. Grids not fulfilling the element-qualitycondition can typically be adapted without loosing their important features by using the meshingtool meshconv , Stenger (2020).The decision to only consider the elements adjacent to the unique closest vertex is a compromisemade for performance-reasons. By extending the nearest neighbor search in the KDTree to severalclosest vertices and considering adjacent elements to all of them one could allow for lower qualitysurface grids at the cost of sacrificing performance.13 h X −→ Γ fine h Figure 3: Coarse grid (left) used as reference domain Γ h for parametrization with closest-pointprojection to fine grid Γ fine h (right) of Stanford-bunny geometry. The Coarse grid is obtained byfeature-preserving coarsening of the fine-grid, see Stenger (2020); Valette and Chassery (2004). If the surface is given implicitly as the zero-level set of a higher-order function, the closest-pointprojection must be calculated iteratively using either a Newton-method or a fixed-point iteration.The essential property that is used in these algorithms, is that the normal vector of the surface ina proximate neighborhood is given by the normalized gradient of the implicit function.Thus, the input to the
ImplicitSurfaceProjection is a differentiable function ψ providingthe surface as Γ = { x ∈ R n : ψ ( x ) = 0 } with n Γ ( x ) = ∇ ψ ( x ) / (cid:107)∇ ψ ( x ) (cid:107) at x ∈ Γ.Given an initial guess for the projected point x , the authors of Persson (2004); Nitschke(2014) describe a scheme to iteratively compute better guesses of the projection of x to Γ byapproximating the closest-point property X ( x ) = x − d ( x ) n Γ ( X ( x )) with a representation of theapproximate distance d ( x ) ≈ ψ ( x ) / (cid:107)∇ ψ ( x ) (cid:107) and the normal vector representation from above: x i +1 = x i − ∇ ψ ( x i ) ψ ( x i ) (cid:107)∇ ψ ( x i ) (cid:107) . (8)This scheme applies this relation iteratively, eventually converging to a point on Γ near the closest-point X ( x ). It is implemented in the class C++ code // < dune / curvedgrid / geometries / implicitsurface .hh > template < class Functor > class SimpleImplicitSurfaceProjection ; template < class Functor > SimpleImplicitSurfaceProjection < Functor > :: SimpleImplicitSurfaceProjection ( const Functor & psi , int maxIter = 10); where the maximal number of iterations in the iterative scheme is given by maxIter .An improved version of that scheme that takes the iterate from the simple scheme as initialapproximation of the closest-point projection and uses this point to get a better estimate for theactual distance of x to Γ is proposed in Demlow and Dziuk (2007):˜ x i +1 = x i − ∇ ψ ( x i ) ψ ( x i ) (cid:107)∇ ψ ( x i ) (cid:107) , dist = sign( ψ ( x )) (cid:107) ˜ x i +1 − x (cid:107) x i +1 = x − ∇ ψ (˜ x i +1 ) dist (cid:107)∇ ψ (˜ x i +1 ) (cid:107) (9)The computational demand is higher than in the simple scheme, but it converges to the actualclosest point on Γ. This scheme is implemented in the class14 ++ code // < dune / curvedgrid / geometries / implicitsurface .hh > template < class Functor > class ImplicitSurfaceProjection ; template < class Functor > ImplicitSurfaceProjection < Functor > :: ImplicitSurfaceProjection ( const Functor & psi , int maxIter = 10); We consider a surface with genus two, given by the function ψ ( x, y, z ) = 2 y ( y − x )(1 − z ) + ( x + y ) − (9 z − − z ) . A reference surface can be obtained by extracting the zero-level set contour, e.g., by usingthe tool ParaView, see Ahrens et al. (2005), by a surface Delaunay triangulation combined witha surface projection, see Persson and Strang (2004); Persson (2004), or by reconstructing theimplicitly defined surface using some fast marching algorithm, see Engwer and N¨ußing (2017). Wefollowed the first approach, combined with a coarsening of the obtained surface grid using Stenger(2020), see also Valette and Chassery (2004) for a similar approach.Applying the implicit projection methods from above to such a coarse reference grid, resultsin very fast convergence of both schemes to a small multiple of the machine epsilon tolerance.iter. error(8) 3 2 . · − (9) 5 2 . · − (a) Iteration counts (b) Grid view (c) Levelset view Figure 4: Surface extracted from the implicit description as zero-level set of ψ , using an implicitprojection method for the higher-order surface approximation. The colored plane illustrates acut through the function ψ . Left-hand side shows a table with the number of iterations and thecorresponding error in the approximate distance, i.e., error = | ψ ( x i ) / (cid:107)∇ ψ ( x i ) (cid:107)| . In all the examples above, a reference grid G h is provided by reading a surface mesh from file. Dune provides a multitude of grid file readers, but lacks support for a reader that can read curvedgeometries directly. This would allow to not only start from a reference grid with an analyticalprojection, but to provide a discrete representation of the curved surface from the beginning. Manymeshing tools allow to directly construct such curved meshes and provide a file format that is ableto represent the additional nodes for a parametrization. We have implemented two readers, the
VtkReader for the VTK file format and a
Gmsh4Reader for the
Gmsh file format. Both associatedmeshing and visualization tools, ParaView and
Gmsh , allow to design curved geometries and toexport these in the mentioned file formats.In addition to file readers, the result of a numerical computation must be exported to allowvisualization and postprocessing. The tool of choice is ParaView, Ahrens et al. (2005), supporting15he VTK file format also for curved geometries. We have implemented a grid and data writer forthis file format.
The module dune-vtk provides grid readers and writers with flexible input and output policies.Writing curved grids requires additional points to be written in the cells of the grid. Those pointscan be associated with Lagrange basis functions allowing a parametrization of the grid also forthe visualization.
VTK supports higher order cell types including Lagrange parametrizations of cells since version9. This allows to directly write the curved geometries to files. A corresponding output policy,called data-collector in dune-vtk , is added to support these cell types. This data-collector isresponsible for transforming a grid-view into a list of point coordinates and a connectivity table.Additionally, it collects values associated to the point coordinates, if data should be written tothe VTK file.For writing higher-order Lagrange parametrized grids, in addition to the corner vertices ofgrid elements internal Lagrange nodes are written to the file. The connectivity table collectsthese nodes in a specific order so that they can be associated to Lagrange basis functions. Thecorresponding data-collector is called LagrangeDataCollector : C++ code // < dune / vtk / datacollectors / lagrangedatacollector .hh > namespace Vtk { template < class GridView , int ORDER = -1> class LagrangeDataCollector ; } The template parameter
GridView represents the type of the grid-view that shall be writ-ten and the second (optional) parameter
ORDER represents the Lagrange polynomial order of thecell parametrization. This second parameter is optional since also runtime polynomial order issupported, by passing the polynomial order parameter to the constructor instead.
C++ code LagrangeDataCollector ( const GridView & gridView , int order = ORDER )
In case no constructor parameter for the order is given, the template
ORDER parameter is usedas default value. Note that either in the template parameter or in the constructor argument apositive value for order must be given.The corresponding writer object can be instantiated by either passing a data-collector objector by letting the writer construct it with the given grid-view object.
C++ code using DataCollector = Vtk :: LagrangeDataCollector < GridView ,4>; using Writer = VtkUnstructuredGridWriter < GridView , DataCollector >; // a) default construct the data - collector with the passed gridView Writer vtkWriter1 { gridView }; // b) construct the data - collector before and pass it to the writer DataCollector dataCollector { gridView }; Writer vtkWriter2 ( dataCollector );
Note, we are using an unstructured-grid writer to generate a .vtu file that represents the grid.16igure 5: Three different approximations of the sphere visualized using the VTK writer withParaView. Shown are the element edges and the Lagrange nodes. Left: reference grid, Center:Lagrange parametrization with polynomial order k = 4, Right: Lagrange parametrization withpolynomial order k = 1 and two grid refinements. When reading a higher-order grid from file, we need to construct both the reference grid andthe parametrization. The reference grid could be obtained by evaluating the higher-order gridrepresentation in the element’s corner vertices, whereas the parametrization must be extractedfrom the additional Lagrange nodes stored in the file.In order to read these nodes and to construct an element connectivity, an input policy called grid-creator — the analog to the data-collector as output policy in the
VtkWriter — must beprovided. It creates local vertex coordinates and element indices from the fields read from file andpasses those to a
GridFactory to create the actual grid.The grid itself does not contain the additional Lagrange nodes, but a parametrization orcoordinate mapping describing the higher-order geometries. Thus, one needs to associate thesenodes to a local Lagrange basis. Since the connection of nodes to corresponding correctly numberedbasis functions might be difficult to make manually by a user of the reader, we provide insteada grid-function representation of the higher-order geometries parametrized over the extractedreference grid. This grid-function is represented by the grid-creator itself.The grid construction, i.e., extracting the reference grid from a potentially higher-order La-grange VTK file, works by just defining the grid type to be constructed. A
GridFactory parametrizedby the grid type is responsible for the actual construction.
C++ code // < dune / vtk / vtkreader .hh > // < dune / vtk / gridcreators / lagrangegridcreator .hh > using Grid = FoamGrid <2,3>; using Creator = Vtk :: LagrangeGridCreator < Grid >; auto grid = VtkReader
C++ code using Grid = FoamGrid <2,3>; GridFactory < Grid > factory ; Vtk :: LagrangeGridCreator creator { factory }; VtkReader reader { creator }; reader . read (" filename . vtu "); // construct the reference grid auto grid = factory . createGrid (); // the creator itself is a grid - function auto localFct = localFunction ( creator ); // traverse the constructed reference - grid for ( auto const & e : elements ( grid -> leafGridView ())) { localFct . bind (e); auto global = localFct ( referenceElement (e). position (0,0)); } Thus, the grid-function can be used to fill any other storage to represent the geometry, e.g. byinterpolating into a
DiscreteGridViewFunction , or can be used directly for the parametrizationof the
CurvedGrid : C++ code CurvedGrid curvedGrid {*grid , creator , creator . order () };
The MSH 4 format implemented in
Gmsh , see Geuzaine and Remacle (2009),supports a multitudeof different cell types. Among them are higher order triangles and quadrilaterals usable for oursurface parametrizations with orders up to 11 (as of version 4.1 of the format). The dune-gmsh4 module provides a file reader that handles most of the features the format provides. Currentlynot supported is the partitioned version of the file format as well as 2d polygonal elements withmore than four vertices.As with the VTK reader, we need to construct both the reference grid and a suitable parametriza-tion when reading a higher-order grid. Therefore, dune-gmsh4 provides a similar
LagrangeGrid-Creator as dune-vtk which represents a grid-function that can be used as parametrization whenconstructing the curved grid. C++ code // < dune / gmsh4 / gmsh4reader .hh > // < dune / gmsh4 / gridcreators / lagrangegridcreator .hh > using Grid = FoamGrid <2,3>; GridFactory < Grid > factory ; Gmsh4 :: LagrangeGridCreator creator { factory }; Gmsh4Reader reader { creator }; reader . read (" filename . msh "); auto grid = factory . createGrid (); CurvedGrid curvedGrid {*grid , creator , creator . order () };
Reading a grid with an implicitly generated grid creator and thus ignoring the parametrization,works analogously to the
VtkReader example. dune-gmsh4 does not implement a file writer for the MSH format. The VTK writer is rec-ommended for output of higher order geometries.
In order to verify the implementation and to test different geometric representations, we have firstanalyzed the difference between the discrete surface and the continuous surface. This numericalverification considers the difference between the smooth surface quantities, the closest-point pro-jection X , the surface normal n , and the mean curvature H , in the L ∞ (Γ h ) norm. In Demlow(2009); Hansbo et al. (2019), upon many others, the following estimates are shown:18 roposition 1 For h small enough, we have the estimates (cid:107) X − X k (cid:107) L ∞ (Γ h ) ≤ Ch k +1 , (cid:107) n ◦ X − n kh (cid:107) L ∞ (Γ kh ) ≤ Ch k , (cid:107) H ◦ X − H kh (cid:107) L ∞ ( e ) ≤ Ch k − (10) for e ∈ G kh , with C a generic constant independent of the mesh parameter h . The corresponding orders are shown for three smooth geometries, the unit sphere, an ellipsoidwith major axis (1 , . , . , k and finally, evaluated the three quantities by iterating over the reference surface. Theerror norms are shown in Figure 6.grid width h grid width h grid width h s ph e r ee lli p s o i d t o r u s k = 1 k = 2 k = 3 h h h h (cid:107) X − X k (cid:107) L ∞ (cid:107) X − X k (cid:107) L (cid:107) n ◦ X − n kh (cid:107) L (cid:107) H ◦ X − H kh (cid:107) L Figure 6: Geometric error norms, normalized by error of largest grid-size, evaluated for threedifferent geometries, the unit sphere, an ellipsoid, and a torus. In various dashed lines, the idealconvergence lines h p are shown. We consider the vector Helmholtz equation for a test of the surface parametrization. The cor-responding intrinsic surface formulation reads: Find the tangential-vector field u ∈ H (Γ , T Γ)such that (cid:0) ∇ Γ u, ∇ Γ v (cid:1) Γ + (cid:0) u, v (cid:1) Γ = (cid:0) f, v (cid:1) Γ , ∀ v ∈ H (Γ , T Γ) , (11)with ∇ Γ the covariant derivative of the vector fields, ( · , · ) Γ the generic L inner product on Γ, and f a given tangential vector field.For the discretization of this equation we follow the ideas of Nestler et al. (2017, 2019); Hansboet al. (2019); Gross et al. (2018); Jankuhn and Reusken (2019) and represent vector field in anembedding space – in this case the Euclidean space T Γ ∼ = R – and transform the corresponding19ovariant derivatives into the embedding space. By allowing the vector field to also have non-tangential components, the equation can be written as a coupled system of scalar-valued equationswith a penalization term to enforce tangentiality: Find the vector field u ∈ [ H (Γ , R )] such that (cid:0) ∇ Γ P u, ∇ Γ P v (cid:1) Γ + (cid:0) P u, P v (cid:1) Γ + ω (cid:0) n · u, n · v (cid:1) Γ = (cid:0) f, P v (cid:1) Γ , ∀ v ∈ [ H (Γ , R )] , (12)with ω (cid:29) P = Id − n ⊗ n the tangential projection operator w.r.t.the surface normal vector n . For extended vector fields u , the surface covariant derivative can beexpressed in terms of the Euclidean derivative ∇ in the ambient space, by ∇ Γ P u = P ∇ ( P u ) P = P ( u ⊗ ∇ Γ ) − ( n ⊗ ∇ Γ )( n · u ). The expression u ⊗ ∇ Γ means the component-wise surface gradient.In order to discretize this equation, we approximate Γ by Γ kh and H (Γ , R ) by V rh,k , the Lagrangefinite-element space of order r , given by V rh,k := (cid:8) v ∈ C (Γ kh ) : v ◦ X ke ∈ P r , ∀ e ∈ G h (cid:9) Following the analysis of Hansbo et al. (2019) and Jankuhn and Reusken (2019); Gross et al.(2018) the normal vector involved in the covariant derivative and mass-matrix term can have thesame approximation order as the geometry, but the normals involved in the penalty term shouldbe at least one order better. We denote this “better” normal by ˜ n . Additionally, the scaling for thepenalization factor should be of order h − , thus we take ω = β/h with β = 10 in the numericalexperiments below.A detailed numerical analysis of this equation for even tensor-valued functions can be foundin Hardering and Praetorius (2020), giving more insight into the penalization term and the errorbehavior for tangential parts of the solution.The resulting discrete variational formulation reads: Find the vector field u h ∈ [ V rh,k ] suchthat (cid:0) ∇ Γ kh P h u h , ∇ Γ kh P h v h (cid:1) Γ kh + (cid:0) P h u h , P h v h (cid:1) Γ kh + βh − (cid:0) ˜ n h · u h , ˜ n h · v h (cid:1) Γ kh = (cid:0) f ◦ X k , P h v h (cid:1) Γ kh , ∀ v ∈ [ V rh,k ] . (13)where the inner product and derivatives have to be understood element-wise and locally, i.e.,( A, B ) Γ kh = (cid:80) e ∈G kh (cid:82) e A : B dΓ. The challenge hereby is to integrate over the parametrized ge-ometries and to provide normal vectors of differing approximation order. Since we have (10), thehigher-order normals can be obtained by constructing a local geometry of parametrization order k + 1. C++ code template < int order > using LFE_t = LagrangeSimplexLocalFiniteElement < double , double ,2, order >; // traverse the reference grid for ( const auto & e : elements ( refGrid -> leafGridView ())) { // projection from local coordinates auto X_e = [& sphere , geo =e. geometry ()]( const auto & local ) { return sphere ( geo . global ( local )); }; // construct the CurvedGeometries from local parametrization ParametrizedGeometry geometry (e. type () , LFE_t
CurvedGrid , one could iterate over the reference grid G h insteadand locally construct the curved geometries of order k and k + 1 to obtain the surface element andthe surface normal vectors. Following the example in Nestler et al. (2019) we define an exact solution u ∗ := rot n ( xyz ) of (11)and construct the corresponding load-vector function f := − div Γ ∇ Γ u ∗ + u ∗ . A coarse grid ofthe sphere is explicitly provided using a Gmsh mesh with nearly equal sized elements.In the following numerical test the geometry parametrization and function parametrizationtake the same polynomial order, i.e., r = k .level grid width h error ( k = 1) eoc error ( k = 2) eoc error ( k = 3) eoc0 7 . · − . · − — 2 . · − — 1 . · − —1 3 . · − . · − .
809 3 . · − .
713 1 . · − . . · − . · − .
000 4 . · − .
121 7 . · − . . · − . · − .
998 5 . · − .
099 4 . · − . . · − . · − .
000 6 . · − .
064 2 . · − . . · − . · − .
000 8 . · − .
036 1 . · − . L -error of linear ( k = 1), quadratic ( k = 2), and cubic ( k = 3) iso-parametric finiteelements for the vector Helmholtz equation with experimental order of convergence (eoc) in thegrid width h that tends towards 2, 3, and 4, respectively. (a) Sphere (b) Ellipsoid Figure 7: Solution of the vector Helmholtz equation on the sphere and on the ellipsoid.
Let Γ h ⊂ R be a smooth closed and stationary reference surface. A map X : Γ h × [0 , T ] → R then defines a parametrization of a family of surfaces Γ( t ) ⊂ R over this reference manifold:Γ( t ) = { X ( x, t ) : x ∈ Γ h } . (14)The evolution of this family of surfaces is characterized by its velocity v ( X, t ) ∈ R at each point X = X ( x, t ) ∈ Γ( t ), ∂ t X ( x, t ) = v ( X ( x, t ) , t ) . (15) The corresponding symbolic computations are done using sympy .
21e consider the surface evolution driven by its mean curvature, the geometric mean curvatureflow flow, see Deckelnick et al. (2005); Kov´acs et al. (2019); Dziuk (1990). Therefore, we introduce H := tr( κ ) the mean curvature of the surface with extended Weingarten map κ = − n ⊗ ∇ Γ , andthe surface evolution v = − Hn Utilizing the geometric identity ∆ Γ X = − Hn , see, e.g., Dziuk (1990); Dziuk and Elliott (2013),a weak formulation of the evolution law can be written: For all t ∈ [0 , T ], find X ( · , t ) ∈ [ H (Γ h )] such that (cid:90) Γ( t ) ∂ t x ( t ) · y dΓ( t ) = − (cid:90) Γ( t ) (cid:88) i (cid:79) Γ( t ) x i ( t ) · (cid:79) Γ( t ) y i dΓ( t ) , ∀ y ∈ [ H (Γ( t ))] . (16)with the surface identity x ( t ) = X ( · , t ). Note that on the right-hand side of that equation, we findthe component-wise surface gradient of the parametrization.In order to discretize this equation, we introduce a splitting of the time-interval [0 , T ] intodiscrete time steps 0 < t < t < . . . < t N = T with generic time step width τ = t s − t s − and denote by X s ∼ = X ( · , t s ) the parametrization at time step t s . Correspondingly, we denote byΓ s = X s (Γ h ) the surface at that time step. Since the grid-function X s is parametrized over thereference surface Γ h we replace the integration over Γ s by an integration over the reference surfaceusing a transformation of the surface elements dΓ h → dΓ s .For the discretization in space, we introduce the finite-element space V rh of Lagrange finite-elements on Γ h , V rh := (cid:8) v ∈ C (Γ h ) : v ◦ X e ∈ P r , ∀ e ∈ G h (cid:9) . Then we get the discrete variational formulation by simple Euler discretization in time: Let X be a given initial parametrization. For all s = 1 , . . . , N , find X s ∈ [ V rh ] such that (cid:90) Γ h ( X s − X s − ) · Y dΓ s − = − (cid:90) Γ h τ (cid:88) i (cid:79) Γ s − X it · (cid:79) Γ s − Y i dΓ s − , ∀ Y ∈ [ V rh ] . (17)So, while traversing the reference grid, we need the geometry of the curved grid from the lasttime step Γ s − = X s − (Γ h ). This is given by the grid-function X s − :Initially we construct a DiscreteGridViewFunction : C++ code auto X = discreteGridViewFunction ( refGrid -> leafGridView () , order ); auto X_e = localFunction (X); // interpolate the initial surface parametrization auto perturbedSphere = []( auto const & x) { return ...; }; Functions :: interpolate (X. basis () , X. coefficients () , perturbedSphere );
This grid-function additionally provides a global basis that can be localized to an element:
C++ code auto localView = X. basis (). localView (); // traverse the reference grid for ( const auto & e : elements (X. basis (). gridView ())) { // bind the local function to the element X_e . bind (e); // bind the localized basis to the element localView . bind (e); // the localized basis provides a local finite - element auto const & localFE = node . child (0). finiteElement (); auto const & localBasis = localFE . localBasis (); // ... }
22n each element of the reference grid, we can construct a curved geometry. This can be used toobtain the integration element and the transform of the local gradients of the local basis-functionsto the actual domain of the curved element.
C++ code LocalFunctionGeometry geometry ( referenceElement (e), X_e ); const auto & quadRule = QuadratureRules < double ,2>:: rule (e. type () , quad_order ); for ( const auto & qp : quadRule ) { // integration element dG_ {s -1} double dS = geometry . integrationElement (qp. position ()) * qp. weight (); // the inverse of the transposed geometry Jacobian auto Jtinv = geometry . jacobianInverseTransposed (qp. position ()); // evaluate the local basis Jacobians in the quadrature point std :: vector < FieldMatrix < double ,1,2>> shapeGradients ; localBasis . evaluateJacobian (qp. position () , shapeGradients ); // transform the local basis Jacobians to the real element std :: vector < FieldVector < double ,3>> gradients ( shapeGradients . size ()); for ( std :: size_t i = 0; i < gradients . size (); ++i) Jtinv .mv( shapeGradients [i][0], gradients [i]); } The evolution of a perturbed initial sphere can be found in Figure 8. It starts the evolution bysmoothing high curvature regions while continuously shrinking the surface. Eventually the surfacegets sphere-like with a radius that tends to zero.Figure 8: Mean-curvature flow of a perturbed spherical surface with parametrization of polynomialorder 2 at four different time-steps in the evolution.
We have implemented parametrized and curved geometries for the
Dune framework by wrappinggrid-functions or differentiable functions into the
Geometry interface defined by the
Dune-Grid module. Additionally, we have implemented a wrapper for flat grids providing a curved geometryon traversal.It is shown in several examples how these wrappers provide high flexibility while preservingsimple usage patterns. In a numerical study we have verified the implementation by showingclassical geometry error bounds to be achieved.The
Dune modules
Dune-CurvedGeometry , Dune-CurvedGrid and additionally
Dune-Vtk and
Dune-Gmsh4 provide not only the geometry and grid wrappers but also utility functionsto simplify the work with curved geometries. These utilities include some reference geometries andsurface projections, as well as grid-functions for various requirements. The file readers and writersfor the VTK format and the
Gmsh format allow for visualization of curved grids and additionally23o construct curved geometries in these tools that can be imported into
Dune grids and thus beused for the internal representation of these surfaces.The implementation of the
Dune modules is neither restricted to only surface parametriza-tions nor to grids without boundaries. Those topics are subject to further investigation and mayrequire additional tools and interface extensions. An example is a non-complete Lagrange ele-ment parametrization with different number of nodes on sub-entities of the same co-dimension.This may allow for boundary adapted parametrizations while the inner elements are quasi flat.Another topic for further studies is the differentiability of the local-functions used in the gridparametrization. Automatic or numeric differentiation of the projections X or X e , as well asimplicit differentiation of the levelset function ψ , could be a way to allow for exact geometryparametrizations of more surfaces. Acknowledgments
This work was supported by German Research Foundation (DFG), Research Unit
Vector- andTensor-Valued Surface PDEs (FOR 3013)
References
J. Ahrens, B. Geveci, and C. Law.
ParaView: An End-User Tool for Large Data Visualization .Visualization Handbook. Elsevier, 2005. ISBN 978-0123875822.P. Bastian, M. Blatt, A. Dedner, N.-A. Dreier, C. Engwer, R. Fritze, C. Gr¨aser, D. Kempf,R. Kl¨ofkorn, M. Ohlberger, and O. Sander. The dune framework: Basic concepts and recentdevelopments. arXiv, 2020. URL https://dune-project.org , 1909.13672.J. L. Blanco and P. K. Rai. nanoflann: a C++ header-only fork of FLANN, a library for nearestneighbor (NN) with kd-trees, 2014. URL https://github.com/jlblancoc/nanoflann .D. Bothe and A. Reusken, editors.
Transport Processes at Fluidic Interfaces . Springer InternationalPublishing, 2017. doi: 10.1007/978-3-319-56602-3.E. Burman, P. Hansbo, M. G. Larson, and A. Massing. Cut finite element methods for partialdifferential equations on embedded manifolds of arbitrary codimensions.
ESAIM: MathematicalModelling and Numerical Analysis , 52(6):2247–2282, 2018. doi: 10.1051/m2an/2018038.K. Deckelnick, G. Dziuk, and C. M. Elliott. Computation of geometric partial differential equationsand mean curvature flow.
Acta Numerica , 14:139–232, 2005. doi: 10.1017/S0962492904000224.A. Dedner and M. Nolte. Construction of local finite element spaces using the generic referenceelements. In A. Dedner, B. Flemisch, and R. Kl¨ofkorn, editors,
Advances in DUNE , pages 3–16,Berlin, Heidelberg, 2012. Springer Berlin Heidelberg. ISBN 978-3-642-28589-9. doi: 10.1007/978-3-642-28589-9 1.A. Demlow. Higher-order finite element methods and pointwise error estimates for ellipticproblems on surfaces.
SIAM Journal on Numerical Analysis , 47(2):805–827, 2009. doi:10.1137/070708135.A. Demlow and G. Dziuk. An adaptive finite element method for the laplace–beltrami operatoron implicitly defined surfaces.
SIAM Journal on Numerical Analysis , 45(1):421–442, 2007. doi:10.1137/050642873.G. Dziuk. An algorithm for evolutionary surfaces.
Numerische Mathematik , 58:603–611, 1990.doi: 10.1007/BF01385643. 24. Dziuk and C. M. Elliott. Finite element methods for surface pdes.
Acta Numerica , 22:289–396,2013. doi: 10.1017/S0962492913000056.C. Engwer and A. N¨ußing. Geometric reconstruction of implicitly defined surfaces and domainswith topological guarantees.
ACM Transactions on Mathematical Software , 44(2), 2017. doi:10.1145/3104989.C. Engwer, C. Gr¨aser, S. M¨uthing, and O. Sander. The interface for functions in the dune-functionsmodule.
Archive of Numerical Software , 5(1):95–110, 2017. doi: 10.11588/ans.2017.1.27683.C. Engwer, C. Gr¨aser, S. M¨uthing, and O. Sander. Function space bases in the dune-functionsmodule. arXiv, 2018, 1806.09545.W. Freeden and M. Schreiner.
Spherical Functions of Mathematical Geosciences . Springer BerlinHeidelberg, 2009. doi: 10.1007/978-3-540-85112-7.H. Fritz. Isoparametric finite element approximation of ricci curvature.
IMA Journal of NumericalAnalysis , 33(4):1265–1290, 2013. doi: 10.1093/imanum/drs037.C. Geuzaine and J.-F. Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities.
International Journal for Numerical Methods in Engineering , 79(11):1309–1331, 2009. doi: 10.1002/nme.2579. URL https://gmsh.info .S. Gross, T. Jankuhn, M. A. Olshanskii, and A. Reusken. A trace finite element method forvector-laplacians on surfaces.
SIAM Journal on Numerical Analysis , 56(4):2406–2429, 2018.doi: 10.1137/17M1146038.P. Hansbo, M. G. Larson, and K. Larsson. Analysis of finite element methods for vector laplacianson surfaces.
IMA Journal of Numerical Analysis , 04 2019. doi: 10.1093/imanum/drz018.H. Hardering and S. Praetorius. Higher-order surface finite-elements for the tensor-laplacian, 2020.in preparation.C. J. Heine. Isoparametric finite element approximation of curvature on hypersurfaces. Technicalreport, Abteilung f¨ur Angewandte Mathematik, Universit¨at Freiburg, Hermann-Herder-Straße10, 79104 Freiburg i.Br., Germany, 2004.T. Jankuhn and A. Reusken. Trace finite element methods for surface vector-laplace equations.arXiv, 2019, 1904.12494.T. Jankuhn, M. Olshanskii, and A. Reusken. Incompressible fluid problems on embedded surfaces:modeling and variational formulations.
Interfaces and Free Boundaries , 20(3):353–377, 2018.doi: 10.4171/IFB/405.B. Kov´acs, B. Li, and C. Lubich. A convergent evolving finite element algorithm for mean cur-vature flow of closed surfaces.
Numerische Mathematik , 143:797–853, 2019. doi: 10.1007/s00211-019-01074-2.M. Nestler, I. Nitschke, S. Praetorius, and A. Voigt. Orientational Order on Surfaces: The Couplingof Topology, Geometry, and Dynamics.
Journal of Nonlinear Science , 28(1):147–191, 2017. doi:10.1007/s00332-017-9405-2.M. Nestler, I. Nitschke, and A. Voigt. A finite element approach for vector- and tensor-valuedsurface pdes.
Journal of Computational Physics , 389:48–61, 2019. doi: 10.1016/j.jcp.2019.03.006.I. Nitschke. Diskretes ¨Außeres Kalk¨ul (DEC) auf Oberfl¨achen ohne Rand. Master’s thesis, Tech-nische Universit¨at Dresden, Dresden, 2014. 25. A. Olshanskii and A. Reusken. Trace Finite Element Methods for PDEs on Surfaces. In
Lecture Notes in Computational Science and Engineering , pages 211–258. Springer InternationalPublishing, 2017. doi: 10.1007/978-3-319-71431-8 7.P.-O. Persson.
Mesh Generation for Implicit Geometries . PhD thesis, Department of Mathematics,MIT, Dec 2004.P.-O. Persson and G. Strang. A simple mesh generator in MATLAB.
SIAM Review , 46(2):329–345,2004. doi: 10.1137/S0036144503429121.S. Praetorius. Dune-Vtk – grid reader and writer for the vtk file format, 2019. URL https://gitlab.mn.tu-dresden.de/iwr/dune-vtk .S. Praetorius. Dune-CurvedGeometry – parametrizations of curved geometries, 2020. URL https://gitlab.mn.tu-dresden.de/iwr/dune-curvedgeometry .S. Praetorius and F. Stenger. Dune-CurvedGrid – meta grid for wrapping element geometries into acurved geometries, 2020a. URL https://gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid .S. Praetorius and F. Stenger. Example codes of this manuscript collected in a dune module, 2020b.URL https://gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid-examples .S. Praetorius and F. Stenger. Dune-Gmsh4 – grid reader and writer for the gmsh-4 file format,2020c. URL https://gitlab.mn.tu-dresden.de/iwr/dune-gmsh4 .A. R¨atz and A. Voigt. PDE’s on surfaces — a diffuse interface approach.
Communicationsin Mathematical Sciences , 4(3):575–590, 09 2006. URL https://projecteuclid.org:443/euclid.cms/1175797557 .O. Sander, T. Koch, N. Schr¨oder, and B. Flemisch. The
Dune-FoamGrid implementation forsurface and network grids.
Archive of Numerical Software , 5(1):217–244, 2017. doi: 10.11588/ans.2017.1.28490.F. Stenger. meshconv: a tool for various mesh-conversions and mesh-transformations., 2020. URL https://gitlab.mn.tu-dresden.de/iwr/meshconv . v3.20.S. Valette and J.-M. Chassery. Approximated centroidal voronoi diagrams for uniform polygonalmesh coarsening.