Multimesh FEM in FEniCS

 

August Johansson and Anders Logg

 

FEniCS'14 in Paris 2014-06-17

 

Simula Research Laboratory
Chalmers University of Technology

Motivation

Methodology

Cut FEM / Nitsche

  • Couple solutions on different meshes
  • Provably stable discretizations
  • Optimal order convergence
  • Based on Nitsche's method / DG
\[ \sum_i (\nabla u, \nabla v)_{\Omega_i} - \sum_j \left \{ (\langle \nabla u \cdot n \rangle, [v])_{\Gamma_j} + (\langle \nabla v \cdot n \rangle, [u])_{\Gamma_j} \right\} \\ + \sum_j \frac{\alpha}{h}([u], [v])_{\Gamma_j} + \sum_k ([\nabla u], [\nabla v])_{O_k} \]

Implementation

Design considerations

  • Independece of CGAL / GTS / external libraries
  • Minimal / orthogonal changes to FEniCS
  • Simple implementation
  • Robust implementation
  • Efficient implementation

New in FEniCS

  • UFL: new measure dc = Measure("custom")
  • UFC: new integral type custom_integral
  • FFC: code generation for custom_integral
  • DOLFIN:
    • Assembly for custom_integral
    • New family of classes MultiMeshFoo
    • Collision detection in 1D, 2D, 3D

UFL


# Custom measures (FIXME: prettify)
dc0 = dc(0, metadata={"num_cells": 1})
dc1 = dc(1, metadata={"num_cells": 2})
dc2 = dc(2, metadata={"num_cells": 2})

# Prettify code
dx = dx + dc0
di = dc1
do = dc2

# Parameters
alpha = 4.0

# Bilinear form
a = dot(grad(u), grad(v))*dx \
  - dot(avg(grad(u)), jump(v, n))*di \
  - dot(avg(grad(v)), jump(u, n))*di \
  + alpha/h*jump(u)*jump(v)*di \
  + dot(jump(grad(u)), jump(grad(v)))*do

# Linear form
L = f*v*dx

UFC



virtual std::size_t num_cells() const = 0;

virtual void tabulate_tensor(double* A,
                             const double * const * w,
                             const double* vertex_coordinates,
                             std::size_t num_quadrature_points,
                             const double* quadrature_points,
                             const double* quadrature_weights,
                             const double* facet_normals,
                             int cell_orientation) const = 0;

 

  • Computes cell tensor A on num_cells() cells
  • One integral class - multiple uses
  • Geometry decided by quadrature_points / weights

FFC


void tabulate_tensor(...)
{
  [...]

  // Set quadrature weights
  const double* W = quadrature_weights;

  [...]

  // Evaluate basis functions on cell 0
  static double FE0_values_0[12];
  for (unsigned int ip = 0; ip < num_quadrature_points; ip++)
  {
    // Get current quadrature point and compute values of basis functions
    const double* x = quadrature_points + ip*2;
    const double* v = vertex_coordinates + 0;
    customvectorintegral_finite_element_1::_evaluate_basis_all(FE0_values_0, x, v, cell_orientation);

    // Copy values to table FE0_C1
    for (std::size_t i = 0; i < 6; i++)
      FE0_C1[ip][0 + i] = FE0_values_0[2*i + 1];

    // Copy values to table FE0_C0
    for (std::size_t i = 0; i < 6; i++)
      FE0_C0[ip][0 + i] = FE0_values_0[2*i + 0];
  }

[...]

  // Evaluate basis function derivatives on cell 0
  static double FE0_dvalues_1_0[24];
  for (unsigned int ip = 0; ip < num_quadrature_points; ip++)
  {
    // Get current quadrature point and compute values of basis function derivatives
    const double* x = quadrature_points + ip*2;
    const double* v = vertex_coordinates + 0;
    customvectorintegral_finite_element_1::_evaluate_basis_derivatives_all(1, FE0_dvalues_1_0, x, v, cell_orientation);

    // Copy values to table FE0_C1_D10
    for (std::size_t i = 0; i < 6; i++)
      FE0_C1_D10[ip][0 + i] = FE0_dvalues_1_0[4*i + 2];

    // Copy values to table FE0_C0_D10
    for (std::size_t i = 0; i < 6; i++)
      FE0_C0_D10[ip][0 + i] = FE0_dvalues_1_0[4*i + 0];

    // Copy values to table FE0_C1_D01
    for (std::size_t i = 0; i < 6; i++)
      FE0_C1_D01[ip][0 + i] = FE0_dvalues_1_0[4*i + 3];

    // Copy values to table FE0_C0_D01
    for (std::size_t i = 0; i < 6; i++)
      FE0_C0_D01[ip][0 + i] = FE0_dvalues_1_0[4*i + 1];
  }

  // Continue with standard quadrature code
  • Call evaluate_basis_[derivatives_]all cells
  • Fill in tables FE0_C0_D10 etc
  • Reuse standard quadrature code

DOLFIN

  • Multimesh assembly
    • uncut_cells, cut_cells, interface, overlap
  • New classes
    • MultiMesh
    • MultiMesh{Function, FunctionSpace, SubSpace}
    • MultiMesh{Assembler, DofMap, Form, DirichletBC}
  • Collision detection
    • BoundingBoxTree
    • IntersectionTriangulation
    • SimplexQuadrature
assemble_cut_cells

void MultiMeshAssembler::assemble_cut_cells(GenericTensor& A,
                                            const MultiMeshForm& a)
{
  // Iterate over parts
  for (std::size_t part = 0; part < a.num_parts(); part++)
  {
    // Get integral for cut cells
    custom_integral = ufc_part.get_custom_integral(0);

    // Get cut cells and quadrature rules
    const auto& cut_cells = multimesh->cut_cells(part);
    const auto& quadrature_rules = multimesh->quadrature_rule_cut_cells(part);

    // Iterate over cut cells
    for (auto it = cut_cells.begin(); it != cut_cells.end(); ++it)
    {
      // Get quadrature rule for cut cell
      const auto& qr = quadrature_rules.at(*it);

      // Tabulate cell tensor
      custom_integral->tabulate_tensor(ufc_part.A.data(),
                                       ufc_part.w(),
                                       vertex_coordinates.data(),
                                       num_quadrature_points,
                                       qr.first.data(),
                                       qr.second.data(),
                                       0,
                                       ufc_cell.orientation);

      // Add entries to global tensor
      A.add(ufc_part.A.data(), dofs);
    }
  }
}

Computational geometry

Essential algorithms

  • Collision detection mesh - point
  • Collision detection mesh - facet
  • Collision detection mesh - cell
  • Collision detection mesh - mesh
  • Quadrature on cut cells
  • Quadrature on interfaces
  • Quadrature on overlaps

Collision detection

  • Axis-aligned bounding box (AABB) trees
  • Implemented by the class hierarchy BoundingBoxTree
    • BoundingBoxTree
    • GenericBoundingBoxTree
    • BoundingBoxTree1D
    • BoundingBoxTree2D
    • BoundingBoxTree3D
  • Inlining of specialized functions

Build tree


GenericBoundingBoxTree::_build()
{
  dolfin_assert(begin < end);

  // Create empty bounding box data
  BBox bbox;

  // Reached leaf
  if (end - begin == 1)
  {
    // Get bounding box coordinates for leaf
    const unsigned int entity_index = *begin;
    const double* b = leaf_bboxes.data() + 2*gdim*entity_index;

    // Store bounding box data
    bbox.child_0 = num_bboxes(); // child_0 == node denotes a leaf
    bbox.child_1 = entity_index; // index of entity contained in leaf
    return add_bbox(bbox, b, gdim);
  }

  // Compute bounding box of all bounding boxes
  double b[MAX_DIM];
  std::size_t axis;
  compute_bbox_of_bboxes(b, axis, leaf_bboxes, begin, end);

  // Sort bounding boxes along longest axis
  auto middle = begin + (end - begin) / 2;
  sort_bboxes(axis, leaf_bboxes, begin, middle, end);

  // Split bounding boxes into two groups and call recursively
  bbox.child_0 = _build(leaf_bboxes, begin, middle, gdim);
  bbox.child_1 = _build(leaf_bboxes, middle, end, gdim);

  // Store bounding box data. Note that root box will be added last.
  return add_bbox(bbox, b, gdim);
}

Search tree


void GenericBoundingBoxTree::_compute_collisions()
{
  // Get bounding box for current node
  const BBox& bbox = tree.get_bbox(node);

  // If point is not in bounding box, then don't search further
  if (!tree.point_in_bbox(point.coordinates(), node))
    return;

  // If box is a leaf (which we know contains the point), then add it
  else if (tree.is_leaf(bbox, node))
  {
    // child_1 denotes entity for leaves
    const unsigned int entity_index = bbox.child_1;

    // If we have a mesh, check that the candidate is really a collision
    if (mesh)
    {
      // Get cell
      Cell cell(*mesh, entity_index);
      if (cell.collides(point))
        entities.push_back(entity_index);
    }

    // Otherwise, add the candidate
    else
      entities.push_back(entity_index);
  }

  // Check both children
  else
  {
    _compute_collisions(tree, point, bbox.child_0, entities, mesh);
    _compute_collisions(tree, point, bbox.child_1, entities, mesh);
  }
}

Collision detection

Quadrature

  • Uncut cells
    trivial (reuse)
  • Interface
    easy
  • Cut cells
    relatively easy
     

Quadrature on cut cells


\[\int_{A \setminus \cup_i B_i} = \int_A - \sum_i \int_{A \cap B_i}\]
  • \(A \setminus \cup_i B_i\) is often nonconvex (meshing NP hard)
  • \(A \cap B_i\) is convex and easy to mesh
  • Even better: \(A \cap B_i\) is a tet-tet intersection
  • Triangulate \(A \cap B_i\) and compute quadrature points

 

Reduce all computational geometry to pairwise
intersection and triangulation of tetrahedra (triangles).

Collision detection

Examples

 

(Preliminary results)

Poisson on 2 meshes

Poisson on 3 meshes

Set background mesh

Add a mesh

Add another mesh

Solution on mesh 0

Solution on mesh 0 + 1

Solution on mesh 0 + 1 + 2

Stokes on 3 meshes

Set background mesh

Add a mesh

Add another mesh

Solution on mesh 0

Solution on mesh 0 + 1

Solution on mesh 0 + 1 + 2

Conclusions

  • Works now:
    • Code generation toolchain (UFL/UFC/FFC)
    • Assembly on cut cells, interfaces, overlaps
    • Scalar and mixed systems

    •  

  • In progress / future work:
    • Further testing, especially for multiple meshes
    • Holes and external interfaces
    • Syntactic sugar, interface improvements
    • Parallel collision detection