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

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.

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.

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.

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.

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.