PolyFEM

The highlevel functions of PolyFEM are in polyfem::State
, a class containing the application state and depending on all the PolyFEM submodules.
polyfem::assembler
assembles the system matrix and right hand side (rhs). The assembler relies on finite element basis (polyfem::basis
), quadrature rules (polyfem::quadrature
), and for transient problems on polyfem::time_integrator
.
polyfem::basis
contains the implementation of finite element basis functions.
The module supports the following 2 dimensional elements:
LagrangeBasis2d.hpp
)LagrangeBasis2d.hpp
)SplineBasis2d.hpp
)PolygonalBasis2d.hpp
)The module supports the following 3 dimensional elements:
LagrangeBasis3d.hpp
)LagrangeBasis3d.hpp
)SplineBasis3d.hpp
)PolygonalBasis3d.hpp
)polyfem::mesh
contains data structures and algorithms for simplicial and polygonal meshes.
The interface is generic (Mesh2D.hpp
, Mesh3D.hpp
). There are two specializations: a conforming implementation that supports conforming polygonal meshes (CMesh2D.hpp
, CMesh3D.hpp
) without hanging nodes, and a nonconforming version (CMesh2D.hpp
, CMesh3D.hpp
) specialized for simplicial meshes only but supporting hanging nodes, which are used to implement adaptive hrefinement.
Many input formats are supported (see the JSON documentation for the complete list), and it outputs meshes in VTU (the format used by paraview/VTK).
In PolyFEM the mesh is used to define a collection of dofs (in case of linear elements, they correspond to the nodes of the mesh, but this is not the case for other elements). The mesh is then discarded and not used by other modules.
polyfem::quadrature
contains a set of quadrature rules for all elements implemented in polyfem::basis
. The quadrature library is mostly independent from the rest of the codebase and it is based on the opensource library quadpy.
For transient problems, polyfem::time_integrator
implements firstorder (Implicit Euler ImplicitEuler.hpp
) and higherorder time integrators (Implicit Newmak ImplicitNewmark.hpp
, BDF BDF.hpp
).
Polyfem contains right hand side, boundary conditions, exact solutions, and initial conditions for many common PDEs. We refer to the documentation in the polyfem::problem
section for more details. The major physics supported by polyfem are:
From these, 1, 3–5, and 7 supports both static and transient problems, and 3–5 model contact and friction forces using the Incremental Potential Contact formulation (https://github.com/ipcsim/ipctoolkit).
polyfem::autogen
contains python scripts to generate C++ code for:
polyfem::basis
polyfem::quadrature
polyfem::assembler
The code autogenerated by the python scripts is directly committed in the code repository, as the autogeneration is time consuming.
Smaller utilities are grouped in the polyfem::utils
module, including mesh IO, hashing functions, logging, matrix utilities, parallelization, rasterization, function interpolation (RBFInterpolation.hpp
), and timers (Timer.hpp
).
In this section, a brief, highlevel code walkthrough for a simple example problem is given.
For this example, we solve the Laplacian on a (triangular) platehole mesh with Dirichlet boundary conditions and a constant RHS. We choose P1 elements implicitly (as this is the default). See the JSON specification for more details.
Internally, PolyFEM keeps a polyfem::State
object that contains the entire state of the solve. The mesh, basis, boundary conditions, stiffness/mass matrices, RHS, and userselected settings are all stored in polyfem::State
. The polyfem::State
class is also responsible for dispatch, meaning it contains member functions for the basic steps of a solve, including building a basis, assembling matrices, and actually performing the linear or nonlinear solve.
After creating a polyfem::State
object and loading the input mesh, the next step in the code is to build the basis. In PolyFEM, there are two separate notions of basis: the finite element (FE) basis and the geometric basis. (They are the same in the isoparametric case.) The former is used to construct the solution whereas the latter is used to represent the geometric mapping, which is especially useful to implement curved elements.
Both the FE and geometric basis are stored as vectors of basis::ElementBases
. Each basis::ElementBases
object stores the basis functions for a given element as a vector of basis::Basis
. Each basis::Basis
stores functions representing a local basis function and its gradient.
The highlevel function that constructs the FE and geometric bases is State::build_basis
. In our example, it uses basis::LagrangeBasis2d
to create the vector of basis::ElementBases
.
In order to evaluate local basis functions and construct mass/stiffness matrices, PolyFEM makes use of a series of data storage classes. The lowest level of these classes is assembler::AssemblyValues
, which stores the evaluation and gradient of a single local basis function on its element's quadrature points. A vector of basis::Local2Global
objects is also stored. This is used to constrain nodes at the interface between different order elements (allowing for refinement).
The assembler::ElementAssemblyValues
class stores per element information, including quadrature points in real space, the geometric mapping's determinant, and a vector of assembler::AssemblyValues
objects (one for each local basis function).
Finally, assembler::AssemblyValsCache
stores a vector of assembler::ElementAssemblyValues
, boosting performance by caching the per element information. Its data can be accessed through its compute
function.
These three classes are the core data storage objects used during matrix assembly. These classes use the previously constructed basis::ElementBases
objects to evaluate basis functions and their gradients.
After building the basis, the next step is to assemble the RHS. This is done through the State::assemble_rhs
function, which itself creates a assembler::RhsAssembler
object. The assembler::RhsAssembler::assemble
function is then called in order to construct the RHS vector of the weak form using assembler::AssemblyValsCache
and the userprovided RHS function. In our simple example, the latter is just $10$. Note that the RHS function and boundary conditions are stored in assembler::GenericProblem
.
The virtual class assembler::Assembler
provides the interface used to actually construct matrices through its assemble
function. For our simple example, the relevant subclass is assembler::LinearAssembler
, which is itself a virtual class that defines how to assemble global matrices. However, the function that defines how to evaluate local matrix entries according to the chosen bilinear form is a pure virtual function. In our example, the stiffness matrix is built using assembler::Laplacian
, which is a subclass of assembler::LinearAssembler
that defines the local matrix entries (ie grad u dot grad v). The assembler::Mass
class is similar except for building the mass matrix.
The relevant highlevel functions for matrix assembly are State::assemble_mass_mat
and State::build_stiffness_mat
. Note that the stiffness matrix is built later on (in the State::solve_linear
function).
After assembling the necessary matrices/vectors, the State::solve_problem
function dispatches State::init_linear_solve
then State::solve_problem
, which solves the linear system.
A diagram showing the relationships discussed above is included below. Each block represents its own class. polyfem::State
is shown in orange and basis::LagrangeBasis2d
in red. The basis storage classes are blue, the value storage classes purple, and the assemblyrelated classes green. Relationships between classes are shown with arrows. Bold arrows represent calling a member function and dashed ones represent inheritance.
In order to boost performance, PolyFEM caches sparsity structure when assembling matrices, taking advantage of the fact that the structure of stiffness/mass matrices remains unchanged assuming the basis is unchanged. The details can be found in MatrixCache.hpp
. A brief overview of polyfem::utils::SparseMatrixCache
is included here for reference.
In PolyFEM, sparse matrices are stored in Compressed Sparse Column (CSC) form. Thus, the matrix is represented as three arrays: two for storing nonzero indices and one for storing values.
Consider the first time a polyfem::utils::SparseMatrixCache
object is filled with values. Each time SparseMatrixCache::add_value
is called, the values and their global indices are placed in a temporary buffer. It is assumed that a given polyfem::utils::SparseMatrixCache
always receives each entry from a given element in the same order.
The first time the matrix is retrieved using SparseMatrixCache::get_matrix
, SparseMatrixCache::prune
is called, which writes the buffered values to the actual matrix. Then, the stored index information is used to create a mapping from the ordering of values received from a given element to the corresponding CSCstyle index. In future iterations, instead of recomputing the matrix structure, polyfem::utils::SparseMatrixCache
uses this cached index to save the given value directly to the CSCstyle values array at the correct position.
All the C++ dependencies required to build the code are included and are compiled statically in a single executable. PolyFEM is compiled using automatic builds on Windows, macOS, and Linux using:
After compilation, unit and system tests (tests/main.cpp
) can be run with:
Polyfem can be added to an existing cmake
project with
and linked with
in your cmake script. We do not currently support other building systems. PolyFEM will download the dependencies that it needs in the build folder.