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.
• Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/363
• Rebased and Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/371
• Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/378
• Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/381
• Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/383
• Merged - https://bitbucket.org/fenics-project/dolfin/pull-requests/385
• Closed - https://bitbucket.org/fenics-project/fiat/pull-requests/38
• Merged - https://bitbucket.org/fenics-project/fiat/pull-requests/39
• Closed - https://bitbucket.org/fenics-project/ffc/pull-requests/73
• Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/77
• Closed - https://bitbucket.org/fenics-project/ffc/pull-requests/82
• Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/83
• Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/84
• Merged - https://bitbucket.org/fenics-project/ffc/pull-requests/85
• Open - https://bitbucket.org/fenics-project/ffc/pull-requests/87
• Merged - https://bitbucket.org/fenics-project/ufl/pull-requests/77
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
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:
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
hex_element clearly has the shape of quadrilateral and hexahedron.
# 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))
# 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.
# 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)))
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
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.
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,
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
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.
• 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
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
• 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.
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.