BVPy: A FEniCS-based Python package to ease the expression and study of boundary value problems in Biology

BVPy is a python library to easily implement and study numerically Boundary Value Problems (BVPs) and Initial Boundary Value Problems (IBVPs) through the Finite Element Method (FEM). BVPy proposes an intuitive Application Programming Interface (API) to harness and combine the core functionalities of three powerful libraries: FEniCS (Logg et al., 2012) provides the core data structures and solving algorithms; Gmsh (Geuzaine & Remacle, 2009) defines the domains and their meshing; and Meshio (Schlömer et al., 2020) handles data reading and writing. Initially built in the context of developmental biology and morphomechanics, its purpose is to enable all users, even with little to none experience in FEM, to quickly and efficiently estimate the behavior of a wide variety of fields (scalars, vectors, tensors) on biologically relevant structures, inspired by biophysical and biochemical processes (morphogene patterning, active matter mechanics, active transports...). Despite this biological motivation, the BVPy library has been implemented in an agnostic manner that makes it suitable for many other scientific context.

• Provide a high-level API to soften the learning curve of FEMs. So users with little to none knowledge can parametrize, run and monitor simulations based on built-in templates. • Enable users experienced in FEM-based modeling to develop and fully customize de novo templates that could, in turn, be used by non-specialists.

State of the field
Although an exhaustive inventory of the existing FEM solutions is far beyond the reach of this paper, it is worth mentioning a few existing alternatives; especially in the context of computational morphomechanics, in which BVPy has been developed.
Some FEM-frameworks provide GUI such as Sofa (Faure et al., 2012), FEBio (Maas et al., 2012) or MorphoMechanX (see (Sapala et al., 2018) for example of use) the yet-to-be published add-on to the MorphoGraphX software (Reuille et al., 2015). While such userfriendly solutions are accessible to a wide range of users; their "monolithic" construction reduces versatility and prevents more experienced users to tune them at will. Moreover, their code sources are not always freely available.
Besides these GUI-based solutions, some "lighter" and more open frameworks, such as FREEFEM (Hecht, 2012) or FEniCS (Logg et al., 2012) are also available. Such solutions appear more versatile and opensource. But their use is usually restricted to people already familiar with the theory of Finite Elements.
A gap currently exists between "closed," userfriendly, softwares on one hand and "open," technical frameworks on the other. BVPy attempts to fill this gap. Our strategy was to harness FEniCS richness into an intuitive and evolutive API. To that end, the library architecture maps the conceptual mathematical components of BVPs to abstract classes built on FEniCS and Gmsh core components (see the library general description below). From these abstract classes, experienced users can implement concrete classes, dedicated to their very specific need. Once implemented, these concrete classes can easily be instantiated and combined in a manner that do not require specific FEM-related skills (see the Example of use below). Although the library comes with a few of these concrete classes, covering classic equations (e.g. Poisson's and Helmholtz's equations, linear elasticity and hyperelasticity, reaction-diffusion equations…); the full potential of the library resides in its ability to integrate de novo classes addressing genuine equations and problems.

Library general description
By definition, a BVP corresponds to a (set) differential equation(s) together with a set of constraints defined on the boundaries of the integration domain, see Equation 1.
The above formalization can be extended as follow in the case of first-order IBVPs (i.e. where time derivatives are limited to first order): the integration domain and its boundaries can be split into one-dimension time-line and orthogonal spatial hyperplanes: Assuming that the spatial hyperplane Ω t0 does not evolve with time, such IBVPs can be implemented by coupling the resolution of time-dependent BVPs with timeintegration schemes.

Scope and range
In terms of equations: In its current version (1.0.0), BVPy can handle non-homogeneous second order Partial Differential Equations (PDEs) with potentially non-linear first order terms but linear higher order ones. The unknowns within these PDEs can be scalar fields, vector fields or second-order tensor fields. Systems of coupled PDEs can be implemented. We extended the initial notion of BVP to encompass IBVPs featuring first-order time derivatives, as described by Equation 2.
In terms of domains: BVPy provides geometrical primitives to generate simple 2D and 3D domains (rectangles, cubes, ellipoids, torus, …) as well as basic Constructive Solid Geometry (CSG) functionalities (addition, substraction and intersection) to combine these primitives into more complex geometries. The library can also handle triangulated, as well as piecewisepolygonal, meshes generated elsewhere. In the current version, .txt , .ply and .obj files are accepted as inputs. This versatility comes however with a drawback, for the domain generation procedure is not yet compatible with parallelization processes.
In terms of boundary conditions: Classic Dirichlet and Neumann boundary conditions can be implemented as well as combinations of these along the domain boundaries. Periodic boundary conditions can also be defined.

Organization
BVPy has been developed and organized as close as possible to the BVP and IBVP mathematical formulations given in Equation 1 and Equation 2. The library is built around the following key components, see Besides these main components, the bvpy.utils module gathers some useful "housekeeping" functions regrouped thematically into sub-modules, e.g. visu, io… See online documentation.

Example of use
The following example is loosely inspired by (Zhao et al., 2020). The idea is to estimate the mechanical stress distribution within a 2D cross section of a pressurized plant tissue with heterogeneous rigidity. The goal here is to illustrate how parsimonious, intuitive and yet insightful, BVPy simulations can be. The recorded results can be processed by third-party software, e.g. Figure 2 shows the visualization of the corresponding stress field within the Paraview software. More examples are available in the tutorial section of the online documentation.