Notebook

Develop assembly of finite element forms on quadrilateral and hexahedral meshes in FEniCS¶

by Ivan Yashchuk¶

This notebook contains the overview of the work done for the Google Summer of Code 2017.

This project involved implementing Lagrange finite elements for quadrilaterals and hexahedrons. The project idea aimed at being able to assemble and solve the simplest PDE, a Poisson's equation, in 2D (quadrilateral mesh) and 3D (hexahedral mesh) in FEniCS. FEniCS Project consists of several repositories, main changes were added to FIAT and FFC. ___

List of the contributions for the GSoC project:¶

Description¶

Many parts to assemble and solve on quad/hex meshes were already in FEniCS, but there were missing links to get the pieces working as a whole.

It was decided to construct Lagrange finite elements on quadrilateral and hexahedron as a tensor product of 1D Lagrange elements on interval, which is mathematically correct.

FIAT has already had a class for constructing TensorProductElements which creates a finite element that is the tensor product of two existing finite elements.

It was possible to create quad and hex element in FIAT:

In [1]:

import FIAT

# Create interval reference cell
interval_ref_el = FIAT.reference_element.UFCInterval()
# Define polynomial degree for the element
degree = 1
# Create 1D Lagrange element
interval_element = FIAT.lagrange.Lagrange(interval_ref_el, degree)
# Create quadrilateral element
quad_element = FIAT.tensor_product.TensorProductElement(interval_element, interval_element)
# Create hexahedral element
hex_element = FIAT.tensor_product.TensorProductElement(quad_element, interval_element)

Reference cell of the quad_element and hex_element clearly has the shape of quadrilateral and hexahedron.

In [2]:

# Print coordinates of reference cell vertices
print(quad_element.ref_el.get_vertices())

((0.0, 0.0), (0.0, 1.0), (1.0, 0.0), (1.0, 1.0))

In [3]:

# Print coordinates of reference cell vertices
print(hex_element.ref_el.get_vertices())

((0.0, 0.0, 0.0), (0.0, 0.0, 1.0), (0.0, 1.0, 0.0), (0.0, 1.0, 1.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (1.0, 1.0, 0.0), (1.0, 1.0, 1.0))

However, FFC automatic code generation did not work for these elements.

The first step was to build an appropriate interface between FFC and FIAT so that C++ code could be generated by FFC to evaluate the tensor product finite element basis functions on the quad/hex cell geometry.

First rough implementation that took care of getting the data from FIAT and computing intermediate representation for the code generation for quadhex case in FFC was submitted in FFC PR #73 and FIAT PR #38. Both declined as the changes were moved to other places in the later pull requests.

These changes made possible assembling the simplest form on quadhex mesh.

In [4]:

# Importing FEniCS components
from fenics import *
# Creating unit square domain devided into 4 x 6 cells
mesh = UnitQuadMesh.create(4, 6)
# Assembling the expression. It should return 1 for unit square domain
print(assemble(1.0 * dx(mesh)))

1.0

Assembling functions defined as Expression and interpolation was fixed with FFC PR #77.

Topological dimension of TensorProductElement's reference cell is stored as the tuple of dimensions in different “directions”, it comes from the tensor product operation, FFC can work only if topological dimension is an integer. This issue is described in my blogpost here. A new wrapper class that reconstructs a FIAT element defined on a TensorProductCell to be defined on Quadrilateral or Hexahedron cells was added with FIAT PR #39. Short description of the class is in my blogpost.

Changes in FFC to support new FIAT class were added in FFC PR #83.

The issue with DirichletBC was fixed with DOLFIN PR #371 and FFC PR #84. The problem was that definition of vertices and facets for quadrilateral and hexahedron cells was different in different parts of FEniCS.

With the above introduced changes it is possible to solve the Poisson's equation from the demo just by changing the mesh to one consisting of quadrilaterals or hexahedrons.

Projection of Constant function onto the Lagrange FunctionSpace for quads and hexes was fixed in FFC PR #85 allowing to solve Stokes or Hyperelasticity type of problems correctly, where it is common to have constant loading.

As the result of this GSoC project, Q and dQ elements (see Periodic Table of the Finite Elements) were implemented in FEniCS. All needed changes are in master branch already. One can try it with any of demos which uses CG or DG FunctionSpace, by changing the mesh to one consisting of quadrilaterals instead of triangles or hexahedrons instead of tetrahedrons.

quadrilateral_mesh = UnitQuadMesh.create(n, m)

hexahedral_mesh = UnitHexMesh.create(n, m , k)

Assembling and solving the PDE on quad and hex mesh also work in parallel!

Weak scaling perfomance with ~640000 DOFs per process was tested for hex mesh and the results are here. ___

Further improvements:¶

• Function evaluation at a point does not work

from fenics import *
mesh = UnitQuadMesh.create(3, 3)
V = FunctionSpace(mesh, 'P', 1)
u = Function(V)
u(0.1, 0.1) # Returns "not implemented" error

It requires evaluate_basis method of the finite element class, which is not implemented. CiarletElement interface shall be implemented for elements on quadrilateral and hexahedron cells. Changes are needed in FIAT/expansions.py, FIAT/polynomial_set.py, ffc/uflacs/backends/ufc/jacobian.py.

• Discontinuous Galerkin Poisson demo does not work as computation of Circumradius is not implemented for quad and hex cells. Therefore h = CellSize(mesh) does not work.

• Many parts of DOLFIN that interact with quadrilateral and hexahedral mesh is not implemented yet. Collision detection, mesh refinement.

• Changing TensorProductElement class to accept list of FiniteElements will simplify the code to generate tensor product element of more than two elements.

To conclude, the original goals of the project are met and exceeded, so I consider GSoC as successfully completed. I hope the community will appreciate my work and contributions :) ___ Some additional information about my journey through the GSoC 2017 can be found in my blog. ___

Demonstration¶

Cahn–Hilliard equation solved on a quadrilateral mesh:

alt text

Hyperelasticity problem solved on a hexahedral mesh: alt text

Lid driven cavity Stokes flow solved on a quadrilateral mesh: alt text

And Poisson's equation solved on a quad mesh: alt text

In [ ]: