PaperlayerView exact revision on arXiv

Source-linked AI summary

High order unfitted finite element methods on level set domains using isoparametric mappings

Christoph Lehrenfeld

arXiv:1509.02762v1math.NAcs.CE

TL;DR

High-order integration is difficult for unfitted finite element methods on implicitly defined level-set domains because standard quadrature is inaccurate on cut elements. The paper combines a piecewise planar unfitted method with a parametric mapping of the underlying mesh, yielding an isoparametric approach whose numerical examples show accurate results for resolved interfaces and robust second-order behavior when interfaces are underresolved.

  • Problem

    Standard quadrature is unsuitable for implicitly defined cut domains, while high-order integration remains a major challenge for unfitted finite element methods.

  • Method

    The method uses an explicit mesh mapping Ψh to improve a piecewise planar level-set interface reconstruction, producing an isoparametric unfitted finite element formulation.

  • Results

    The numerical examples show high accuracy for smooth, well-resolved interfaces and robust standard second-order accuracy when interfaces are not well resolved.

  • Takeaways & Limitations

    The geometry-based approach supports high-order unfitted computations while remaining robust and fairly simple to implement across interface, boundary, and surface problems.

  • Takeaways & Limitations

    A rigorous error analysis is missing, and stable discretizations and suitable linear solvers for high-order unfitted methods remain difficult to develop and analyze.

Abstract

from arXiv · show

We introduce a new class of unfitted finite element methods with high order accurate numerical integration over curved surfaces and volumes which are only implicitly defined by level set functions. An unfitted finite element method which is suitable for the case of piecewise planar interfaces is combined with a parametric mapping of the underlying mesh resulting in an isoparametric unfitted finite element method. The parametric mapping is constructed in a way such that the quality of the piecewise planar interface reconstruction is significantly improved allowing for high order accurate computations of (unfitted) domain and surface integrals. This approach is new. We present the method, discuss implementational aspects and present numerical examples which demonstrate the quality and potential of this method.

1. Introduction

The paper introduces a geometry-based isoparametric unfitted finite element approach for high-order integration on implicitly defined level-set domains. It improves piecewise planar interface reconstructions through a mesh mapping while preserving explicit representations and practical implementation features.

  • Standard quadrature on implicitly defined cut domains is inadequate because cut-element integrands lack the smoothness required for accurate integration.
  • The proposed method combines a piecewise planar unfitted discretization with a parametric mapping of the underlying mesh, rather than the sub-triangulation or interface.The resulting finite element formulation is naturally isoparametric.
  • The mapping Ψh transforms the piecewise planar interface and its associated subdomains into explicit high-order approximations of the desired level-set geometry.Quadrature points and weights from the planar configuration can therefore be reused after transformation.
  • The approach provides an explicit high-order geometry approximation, positive interface and subdomain quadrature weights, unchanged cut topology, and straightforward integration into existing codes.
  • The paper presents the approach and implementation aspects but does not aim at a thorough error analysis, which remains ongoing research.
  • The construction of Ψh is designed for implementation and is developed explicitly for simplicial meshes with an approximate signed distance level-set function.

2. Construction of the mesh transformation Ψh

The transformation Ψh is designed to map a piecewise planar interface toward the target interface while preserving element quality and locality. Its construction uses an explicit mapping, finite-element approximation, and robustness controls.

  • Transformation requirements: Ψh is sought in a continuous vector-valued piecewise-polynomial space and must approximate the target interface with distance O(h^{k+1}).The admissible transformation also satisfies shape-regularity and locality constraints.
  • Transformation requirements: The mapped elements should retain shape regularity comparable to the original mesh under the composition of affine and parametric transformations.A constant C>1 bounds the transformed-element regularity relative to the original element.
  • Transformation requirements: The deformation is restricted to the vicinity of the interface, while the transformation remains the identity away from the interface region.Locality is desirable for efficiency but is not essential to the construction.
  • Construction: The deformation is essentially a local high-order correction because the piecewise linear interface is already a second-order approximation.The requirements do not uniquely determine Ψh, so different choices are possible.
  • Construction: For simplicial meshes and an approximate signed distance level set, the paper constructs Ψh explicitly from an ideal mapping, a modification, and a finite-element approximation.The construction targets the conditions for interface approximation, identity away from the interface, and shape regularity.
  • Construction: The algorithm for determining Ψh and its principal computational properties is summarized in section 2.7.

2.1. Notation and assumptions

The analysis assumes a smooth interface represented implicitly by a signed distance function and approximated on a shape-regular simplicial mesh. The exact level set is projected into a continuous piecewise-polynomial space, while its linear interpolant defines the planar interface.

  • Assumptions: The interface is assumed C^m-smooth with m≥2, and the mesh size satisfies h<δ within a tubular neighborhood of the interface.These assumptions provide the neighborhood in which the exact level-set function is sufficiently regular.
  • Level-set approximation: The exact level-set function is approximated by φ∈V_h^k, with 1≤k≤m−1, through a suitable projection operator.In practice, only φ may be known, so approximation-error assumptions are imposed on the projection.
  • Level-set approximation: The projection error is assumed to satisfy L2 and uniform bounds of order h^{k+1} in the interface neighborhood.The constants depend only on k and Γ.
  • Planar reconstruction: The zero level Γ1={Ihφ=0} of the nodal linear interpolant is piecewise planar and provides an explicit representation of the interface and subdomains.Γ1 is a second-order accurate approximation to Γ.
  • Cut-element regions: Cut elements form TΓ and their region ΩΓ, with neighboring elements extending the computational region to TΓ,+ and ΩΓ,+.

2.2. Ideal transformation Ψ

The ideal mapping Ψ transports piecewise-linear level sets to corresponding level sets of the exact level-set function. In the non-ideal case, discontinuous derivatives can make Ψ nonsmooth, motivating a smoothing or localization modification.

  • Ideal construction: The ideal construction assumes φ=φex, so the level-set gradient is continuous in the cut-element region and Ψ is defined there.The ideal case also identifies the target interface Γ with the zero level of φ.
  • Level-set mapping: Ψ maps every level set {Ihφ=c} of the piecewise-linear approximation onto the corresponding level set {φ=c}.Consequently, the planar interface Γ1 is mapped onto the target interface Γ.
  • Pointwise definition: For each point x, the mapped point x*=x+r(x)s is selected on the matching level set using the shortest-distance search direction.The signed-distance structure makes the corresponding point unique under the stated neighborhood and smoothness assumptions.
  • Ideal regularity: When φ=φex, continuity of ∇φex yields a continuous ideal transformation, and mesh vertices remain fixed because their interpolated and exact level-set values coincide.
  • Non-ideal case: For φ≠φex, discontinuities of derivatives across element interfaces can make Ψ only Lipschitz continuous within elements and introduce shifted kinks.The search direction may be projected for continuity, or discontinuity may be accepted before the final finite-element projection.
  • Evaluation: Evaluating the pointwise mapping by Newton iteration can require neighboring-element evaluations and encounter level-set kinks, complicating nonlinear solution and parallel implementation.A subsequent modification addresses both issues.

2.3. Localized (piecewise smooth) transformation

The localized transformation ΨT replaces the level-set function on each cut element by a polynomial extension, producing elementwise smooth mappings and avoiding neighbor-element evaluations. Projection then yields a continuous finite-element deformation Ψh.

  • Localized mapping: The localized construction defines ΨT elementwise on cut elements using the polynomial extension ET(φ) of the level-set function.This avoids evaluating finite-element functions from neighboring elements.
  • Localized mapping: Using ET(φ) in place of φ makes ΨT smooth within each element, as illustrated by the localized transformation figure.The construction maps planar level-set points to corresponding points on locally extended level sets.
  • Approximation: For an interface-resolved mesh, the extension error in a neighborhood ω(T) is bounded by O(h^{k+1}).The bound combines the approximation error of φ with a Bramble-Hilbert estimate for the extension.
  • Finite-element projection: Although ΨT is discontinuous across element boundaries, its jumps are expected to be O(h^{k+1}).A projection ΠΨ is used to obtain a continuous finite-element transformation Ψh.
  • Finite-element projection: The deformation is represented by a finite-element function Ψh∈[V_h^k]^d whose construction can use nodal interpolation, an L2 projection, or another projection.The paper states that rigorous justification of the associated assumption remains ongoing research.
  • Error estimate: The mapped interface Γh:=Ψh(Γ1) is expected to satisfy the same error estimate, with the localized case differing by an additional interpolation term.For ΨT, the term involving φ∘ΨT−Ihφ is not exactly zero but is estimated separately.

2.5. An Oswald-type projection

The Oswald-type projection avoids global linear systems by projecting into a discontinuous polynomial space element by element, then averaging into the continuous finite element space.

  • The construction applies to scalar or vector-valued projections on volume or interface regions, including Πφ, Πs, and ΠΨ.
  • The projection consists of a discontinuous finite-element projection followed by averaging into V^k_h.
  • The discontinuous space contains piecewise polynomials that may be discontinuous across element boundaries.
  • Because continuity constraints are absent, the L2 projection can be computed element by element and implemented efficiently.
  • The averaging operator is a high-order version of the Oswald interpolation operator and averages coefficients over elements supporting each basis function.

2.6. Shape regularity

Shape regularity is straightforward for smooth, well-resolved interfaces but requires separate treatment when the interface is non-smooth or under-resolved.

  • The analysis distinguishes well-resolved smooth interfaces from non-smooth or under-resolved interfaces when assessing shape regularity.
  • The method seeks conditions that ensure shape regularity even when the interface is not sufficiently resolved by the mesh.

Simplex transformations:.

The transformed element is constructed through affine and nonlinear mappings, and its shape regularity is controlled through the Jacobian of the combined transformation.

  • The affine map Φ_h transforms the reference simplex into a planar simplex, while Ψ_h produces the final curved element.
  • The final element satisfies T = Θ_h(ˆT) with Θ_h = Ψ_h ◦ Φ_h, combining the curving and affine transformations.
  • The mappings fix the simplex vertices, so ˆΨ_h(ˆx_V) = ˆx_V and Ψ_h(x_V) = x_V.
  • If the reference curving map and affine map are well behaved, the combined transformation is also well behaved.
  • Shape regularity is quantified using the relative spectral condition number of the Jacobian, with the initial mesh assumed shape regular and the nonlinear contribution controlled through κ(∇ˆΨ_h).
  • For sufficiently fine interface resolution, the deformation becomes arbitrarily small, ensuring shape regularity of the transformed mesh.

The resolved case:.

The barrier step prevents poor element deformations in under-resolved regions, but using it limits the geometry approximation to essentially second order.

  • Under-resolved interfaces can cause self-intersections, arbitrarily small angles, condition-number blow-up, or singular Jacobians.
  • A barrier step limits the deformation to guarantee the quality of the resulting mesh.
  • For sufficiently small γ, the limitation step ensures shape regularity independently of interface resolution.
  • Using the limitation step reduces the geometry approximation to essentially second-order accuracy when the limited deformation differs from the original deformation.
  • When the interface is resolved, the limitation is unnecessary; in the worst under-resolved case, the method essentially recovers piecewise-linear approximation quality while retaining shape-regular deformed elements.

2.7. Summary of computational aspects

The computational procedure constructs the deformation Ψh elementwise, limits excessive displacements, and averages local contributions into global coefficients. The resulting mesh changes only near the interface, can achieve high-order interface resolution when sufficiently resolved, and remains shape regular.

  • Algorithm: The algorithm computes element contributions to Ψh at integration points, including Newton evaluation and displacement limitation, then assembles and averages global coefficient vectors.The search direction is computed separately using an analogous projection procedure.
  • Deformation properties: Only elements near Γ are deformed; the mesh remains unchanged in Ω\ΩΓ,+.
  • Deformation properties: Where the interface is sufficiently resolved, the deformation represents it with high-order accuracy.
  • Deformation properties: Shape regularity is guaranteed in all cases, including unresolved-interface configurations.In the worst case, the method recovers the quality of the piecewise linear approximation while preserving shape regularity.

3. Numerical examples I: Geometry approximation

Geometry experiments show high-order convergence once interfaces are sufficiently resolved, while deformation limiting preserves robustness on coarse meshes. The three-dimensional gyroid exhibits a longer pre-asymptotic regime but ultimately reaches O(h^(k+1)) behavior.

  • Setup: The implementation uses γ = 0.1, a 10^-14 nonlinear tolerance, a continuous piecewise linear search field, and the projected transformation Ψh = Πo(ΨT).The nonlinear solve required only 2–3 iterations per point in the reported examples.
  • Error measure: The geometrical error is the maximum distance between the discrete interface Γh and exact interface Γ, combining level-set and interface-reconstruction errors.The reported quantity is dist(Γh, Γ).
  • Flower-shape example: For the flower example, sufficiently fine interface resolution yields the optimal geometrical-error order k + 1, whereas coarse levels 0–2 lack high-order convergence because deformation is limited.
  • Gyroid example: For the gyroid, the asymptotic geometrical-error rate is O(h^(k+1)), but the pre-asymptotic regime extends through refinement level 4 with only O(h^2) decay.Once the piecewise linear interface resolves the geometry, deformation limitation is no longer necessary; increasing k by one reduces the finest-level error by two orders of magnitude.
  • Overall findings: The experiments conclude that the method is highly accurate for smooth, well-resolved interfaces and retains standard second-order accuracy when interfaces are not well resolved.

4. Numerical examples II: A high order unfitted isoparametric finite element method for an interface problem

The paper combines an isoparametric transformation with an extended finite element and Nitsche discretization for an unfitted interface problem. The method achieves optimal convergence, while high-order unfitted solvers remain numerically challenging and the method lacks a rigorous error analysis.

  • Discretization: The interface discretization combines an extended finite element space, an isoparametric mapping, and a Nitsche formulation for the unfitted interface problem.
  • Stabilization: Heaviside averaging weights based on the undeformed cut configuration yield stability for arbitrary polynomial degrees, independently of interface cut position.
  • Limitations: Efficient robust linear solvers for high-order unfitted finite element systems remain a challenging open problem.The paper notes that rigorous robust solver analysis currently covers only the k = 1 case cited there.
  • Numerical results: The geometry error is expected to satisfy O(h^(k+1)), while the volume errors show predicted optimal convergence rates and the interface-condition error converges with O(h^(k+1)).The interface-condition rate is reported as half an order better than predicted.
  • Numerical results: The method is combined with a standard Nitsche-XFEM discretization and exhibits optimal-order convergence for an unfitted interface problem.
  • Limitations: A rigorous error analysis is missing, and extensions to non-simplex meshes and moving interfaces remain future work.
Loading 1509.02762v1…