Log In Sign Up

Veamy: an extensible object-oriented C++ library for the virtual element method

This paper summarizes the development of Veamy, an object-oriented C++ library for the virtual element method (VEM) on general polygonal meshes, whose modular design is focused on its extensibility. The two-dimensional linear elastostatic problem has been chosen as the starting stage for the development of this library. The theory of the VEM in which Veamy is based upon is presented using a notation and a terminology that resemble the language of the finite element method in engineering analysis. Several examples are provided to demonstrate the usage of Veamy, and in particular, one of them features the interaction between Veamy and the polygonal mesh generator PolyMesher. Veamy is free and open source software.


page 37

page 38


A Velocity-based Moving Mesh Virtual Element Method

We present a velocity-based moving mesh virtual element method for the n...

Efficient multi-level hp-finite elements in arbitrary dimensions

We present an efficient algorithmic framework for constructing multi-lev...

SplineLib: A Modern Multi-Purpose C++ Spline Library

This paper provides the description of a novel, multi-purpose spline lib...

The deal.II finite element library: design, features, and insights

deal.II is a state-of-the-art finite element library focused on generali...

The deal.II finite element library

deal.II is a state-of-the-art finite element library focused on generali...

Medusa: A C++ Library for solving PDEs using Strong Form Mesh-Free methods

Medusa, a novel library for implementation of strong form mesh-free meth...

1 Introduction

When the Galerkin weak formulation of a boundary-value problem such as the linear elastostatic problem is solved numerically, the trial and test displacements are replaced by their discrete representations using basis functions. Herein, we consider basis functions that span the space of functions of degree 1 (i.e., affine functions). Due to the nature of some basis functions, the discrete trial and test displacement fields may represent linear fields plus some additional functions that are non-polynomials or high-order monomials. Such additional terms cause inhomogeneous deformations, and when present, integration errors appear in the numerical integration of the stiffness matrix leading to stability issues that affect the convergence of the approximation method. This is the case of polygonal and polyhedral finite element methods Talischi-Paulino:2014 ; talischi:2015 ; francis:LSP:2016 , and meshfree Galerkin methods dolbow:1999:NIO ; chen:2001:ASC ; babuska:2008:QMM ; babuska:2009:NIM ; ortiz:2010:MEM ; ortiz:2011:MEI ; duan:2012:SOI ; duan:2014:CEF ; duan:2014:FPI ; ortiz:VAN:2015 ; ortiz:IRN:2015 .

The virtual element method BeiraoDaVeiga-Brezzi-Cangiani-Manzini-Marini-Russo:2013 (VEM) has been presented to deal with these integration issues. In short, the method consists in the construction of an algebraic (exact) representation of the stiffness matrix without the explicit evaluation of basis functions (basis functions are virtual). In the VEM, the stiffness matrix is decomposed into two parts: a consistent matrix that guarantees the exact reproduction of a linear displacement field and a correction matrix that provides stability. Such a decomposition is formulated in the spirit of the Lax equivalence theorem (consistency stability convergence) for finite-difference schemes and is sufficient for the method to pass the patch test cangiani:2015 . Recently, the virtual element framework has been used to correct integration errors in polygonal finite element methods Gain-Talischi-Paulino:2014 ; Manzini-Russo-Sukumar:2014 ; BeiraodaVeiga-Lovadina-Mora:2015 and in meshfree Galerkin methods ortiz:CSMGMVEM:2017 .

Some of the advantages that the VEM exhibits over the standard finite element method (FEM) are:

  • Ability to perform simulations using meshes formed by elements with arbitrary number of edges, not necessarily convex, having coplanar edges and collapsing nodes, while retaining the same approximation properties of the FEM.

  • Possibility of formulating high-order approximations with arbitrary order of global regularity daVeiga:VEMAR:2013 .

  • Adaptive mesh refinement techniques are greatly facilitated since hanging nodes become automatically included as elements with coplanar edges are accepted cangiani:PEE:2017 .

In this paper, object-oriented programming concepts are adopted to develop a C++ library, named Veamy, that implements the VEM on general polygonal meshes. The current status of this library has a focus on the linear elastostatic and Poisson problems in two dimensions, but its design is geared towards its extensibility. Veamy uses Eigen library eigenweb for linear algebra, and Triangle shewchuk96b and Clipper clipperweb are used for the implementation of its polygonal mesh generator, Delynoi delynoiweb , which is based on the constrained Voronoi diagram. Despite this built-in polygonal mesh generator, Veamy is capable of interacting straightforwardly with PolyMesher Talischi:POLYM:2012 , a polygonal mesh generator that is widely used in the VEM and polygonal finite elements communities.

In presenting the theory of the VEM, upon which Veamy is built, we adopt a notation and a terminology that resemble the language of the FEM in engineering analysis. The work of Gain et al. Gain-Talischi-Paulino:2014 is in line with this aim and has inspired most of the notation and terminology used in this paper.

In Veamy

’s programming philosophy entities commonly found in the VEM and FEM literature such as mesh, degree of freedom, element, element stiffness matrix and element force vector, are represented by objects. In contrast to some of the well-established free and open source object-oriented FEM codes such as FreeFEM++ 

hecht:FFEM:2012 , FEniCS alnaes:FENI:2015 and Feel++ prud:FEEL:2012 , Veamy does not generate code from the variational form of a particular problem, since that kind of software design tends to hide the implementation details that are fundamental to understand the method. On the contrary, since Veamy’s scope is research and teaching, in its design we wanted a direct and balanced correspondence between theory and implementation. In this sense, Veamy is very similar in its spirit to the 50-line MATLAB implementation of the VEM Sutton:VEM:2017 . However, compared to this MATLAB implementation, Veamy is an improvement in the following aspects:

  • Its core VEM numerical implementation is entirely built on free and open source libraries.

  • It offers the possibility of using a built-in polygonal mesh generator, whose implementation is also entirely built on free and open source libraries. In addition, it allows a straightforward interaction with PolyMesher Talischi:POLYM:2012 , a popular and widely used MATLAB-based polygonal mesh generator.

  • It is designed using the object-oriented paradigm, which allows a safer and better code design, facilitates code reuse and recycling, code maintenance, and therefore code extension.

  • Its initial release implements both the two-dimensional linear elastostatic problem and the two-dimensional Poisson problem.

We are also aware of the MATLAB Reservoir Simulation Toolbox lie:MRS:2016 , which provides a module for first- and second-order virtual element methods for Poisson-type flow equations that was developed as part of a master thesis klemetsdal:VEM:2016 . The toolbox also implements a module dedicated to the VEM in linear elasticity for geomechanics simulations.

Veamy is free and open source software, and to the best of our knowledge is the first object-oriented C++ implementation of the VEM.

The remainder of this paper is structured as follows. The model problem for two-dimensional linear elastostatics is presented in Section 2. Section 3 summarizes the theoretical framework of the VEM for the two-dimensional linear elastostatic problem. Also in this section, the VEM element stiffness matrix for the two-dimensional Poisson problem is given. The object-oriented implementation of Veamy is described and explained in Section 4. In Section 5, some guidelines for the usage of Veamy’s built-in polygonal mesh generator are given. Several examples that demonstrate the usage of Veamy and a performance comparison between VEM and FEM are presented in Section 6. The paper ends with some concluding remarks in Section 7.

2 Model problem

The Galerkin weak formulation for the linear elastostatic problem is considered for presenting the main ingredients of the VEM. Consider an elastic body that occupies the open domain and is bounded by the one-dimensional surface whose unit outward normal is . The boundary is assumed to admit decompositions and , where is the essential (Dirichlet) boundary and is the natural (Neumann) boundary. The closure of the domain is . Let be the displacement field at a point of the elastic body when the body is subjected to external tractions and body forces . The imposed essential (Dirichlet) boundary conditions are . The Galerkin weak formulation, with being the arbitrary test function, gives the following expression for the bilinear form:



is the Cauchy stress tensor and

is the gradient operator. The gradient of the displacement field can be decomposed into its symmetric (

) and skew-symmetric (

) parts, as follows:




is the strain tensor, and


is the skew-symmetric gradient tensor that represents rotations. The Cauchy stress tensor is related to the strain tensor by


where is a fourth-order constant tensor that depends on the material of the elastic body.

Substituting (2) into (1) and noting that because of the symmetry of the stress tensor, results in the following simplification of the bilinear form:


which leads to the standard form of presenting the weak formulation: find such that


where and are the displacement trial and test spaces defined as follows:

where the space includes linear displacement fields.

In the Galerkin approximation, the domain is partitioned into disjoint non overlapping elements. This partition is known as a mesh. We denote by an element having an area of and a boundary that is formed by edges of length . The partition formed by these elements is denoted by , where is the maximum diameter of any element in the partition. The set formed by the union of all the element edges in this partition is denoted by , and the set formed by all the element edges lying on is denoted by . On this partition, the trial and test displacement fields are approximated using basis functions, and hence and are replaced by the approximations and , respectively. The bilinear and linear forms are then obtained by summation of the contributions from the elements in the mesh, as follows:

In general, the weak form integrals are not available in closed-form expressions since functions in , and in particular its basis, are not necessarily polynomial functions. Therefore, these integrals are evaluated using quadrature with the potential of introducing quadrature errors making them mesh-dependent. If that is the case, the convergence of the numerical solution will be affected. To reflect this, a superscript is added to the symbols that represent the bilinear and linear forms. Thus, the Galerkin solution is sought as the solution of the global system that results from the weak formulation described by the following discrete bilinear and linear forms:

respectively, with the corresponding discrete global trial and test spaces defined respectively as follows:

In the preceding discussion, we have implied that is inexact due to its evaluation using numerical quadrature — in this case, is said to be uncomputable. The situation is completely different in the VEM approach: is not evaluated using numerical quadrature. Instead, the displacement field is computed through projection operators that are tailored to achieve an algebraic (exact) evaluation of — in this case, is said to be computable.

3 The virtual element method

In standard two-dimensional finite element methods, the partition is usually formed by triangles and quadrilaterals. In the VEM, the partition is formed by elements with arbitrary number of edges, where triangles and quadrilaterals are particular instances. We refer to these more general elements as polygonal elements.

3.1 The polygonal element

Let the domain be partitioned into disjoint non overlapping polygonal elements with straight edges. The number of edges and nodes of a polygonal element are denoted by . The unit outward normal to the element boundary in the Cartesian coordinate system is denoted by . Fig. 1 presents a schematic representation of a polygonal element for , where the edge of length and the edge of length are the element edges incident to node , and and are the unit outward normals to these edges, respectively.

Figure 1: Schematic representation of a polygonal element of edges

3.2 Projection operators

As in finite elements, for the numerical solution to converge monotonically it is required that the displacement approximation in the polygonal element can represent rigid body modes and constant strain states. This demands that the displacement approximation in the element is at least a linear polynomial strang:2008:AAO . In the VEM, projection operators are devised to extract the rigid body modes, the constant strain states and the linear polynomial part of the motion at the element level. The spaces where these components of the motion reside are given next.

The space of linear displacements over is defined as


where is defined through the mean value of a function over the element nodes given by


where is the number of nodes of coordinates that define the polygonal element111Eq. (9) in fact defines any ‘barred’ term that appears in this paper.; is a second-order tensor and thus can be uniquely expressed as the sum of a symmetric and a skew-symmetric tensor. Let the symmetric and skew-symmetric tensors be denoted by and , respectively. The spaces of rigid body modes and constant strain states over are defined, respectively, as follows:


Note that the space of linear displacements is the direct sum of the spaces given in (10) and (11), that is, .

The extraction of the components of the displacement field in the three aforementioned spaces is achieved through the following projection operators:


for extracting the rigid body modes,


for extracting the constant strain states, and


for extracting the linear polynomial part. And since , the projection operators satisfy the relation


We know by definition that the space includes linear displacements. This means that . Thus, any can be decomposed into three terms, as follows:


that is, into their rigid body modes, their constant strain states and their additional non-polynomial or high-order functions, respectively.

The explicit forms of the projection operators that are defined through (12)-(14) are given in Ref. Gain-Talischi-Paulino:2014 and are summarized as follows: let the cell-average of the strain tensor be defined as


where the divergence theorem has been used to transform the volume integral into a surface integral. Similarly, the cell-average of the skew-symmetric gradient tensor is defined as


Note that and are constant tensors in the element.

On using the preceding definitions, the projection of onto the space of rigid body modes is written as


where and are the rotation and translation modes of , respectively. And the projection of onto the space of constant strain states is given by


Hence, by (15) the projection of onto the space of linear displacements is written as


The projection operator satisfies some important energy-orthogonality conditions that are invoked when constructing the VEM bilinear form. The proofs can be found in Ref. Gain-Talischi-Paulino:2014 . The energy-orthogonality conditions are given next.

The projection satisfies:


The condition (22a) means that is energetically orthogonal to . The condition (22b) emanates from condition (22a) after replacing and using the fact that rigid body modes have zero strain, that is .

3.3 The VEM bilinear form

Substituting the VEM decomposition (16) into the bilinear form (6) leads to the following splitting of the bilinear form at element level:


where the symmetry of the bilinear form, the fact that and do not contribute in the bilinear form (both have zero strain as they belong to the space of rigid body modes), and the energy-orthogonality condition (22b) have been used.

The first term on the right-hand side of (3.3) is the bilinear form associated with the constant strain states that provides consistency (it leads to the consistency stiffness) and the second term is the bilinear form associated with the additional non-polynomial or high-order functions that provides stability (it leads to the stability stiffness). We come back to these concepts later in this section.

3.4 Projection matrices

The projection matrices are constructed by discretizing the projection operators. We begin by writing the projections and in terms of their space basis. To this end, consider the two-dimensional Cartesian space and the skew-symmetry of 333Note that and .. The projection (19) can be written as follows:


where the basis for the space of rigid body modes is:


Similarly, on considering the symmetry of , the projection (20) can be written as


where the basis for the space of constant strain states is:


On each polygonal element of edges with nodal coordinates denoted by , the trial and test displacements are locally approximated as


where and are assumed to be the canonical basis functions having the Kronecker delta property (i.e., Lagrange-type functions), and and are nodal displacements. The canonical basis functions are also used to locally approximate the components of the basis for the space of rigid body modes:


and the components of the basis for the space of constant strain states:


The discrete version of the projection to extract the rigid body modes is obtained by substituting (28) and (29) into (24), which yields










In (36), appeared because of the discretization of (see (18)) and is given by


Similarly, substituting (28) and (30) into (26) leads to the following discrete version of the projection to extract the constant strain states:








In (41), is also given by (37) but in this case it stems from the discretization of (see (17)).

The matrix form of the projection to extract the polynomial part of the displacement field is then .

For the development of the element consistency stiffness matrix, it will be useful to have the following alternative expression for the discrete projection to extract the constant strain states:


3.5 VEM element stiffness matrix

The decomposition given in (3.3) is used to construct the approximate mesh-dependent bilinear form in a way that is computable at the element level. To this end, we approximate the quantity , which is uncomputable, with a computable one given by and define


where its right-hand side, as it will be revealed in the sequel, is computed algebraically. The decomposition (43) has been proved to be endowed with the following crucial properties for establishing convergence BeiraoDaVeiga-Brezzi-Cangiani-Manzini-Marini-Russo:2013 ; BeiraodaVeiga-Brezzi-Marini:2013 :

For all and for all in

  • Consistency: and

  • Stability: two constants and , independent of and of , such that


The discrete version of the VEM element bilinear form (43) is constructed as follows. Substitute (3.4) into the first term of the right-hand side of (43) (note that when is used instead of , is replaced by the column vector of nodal displacements , which has the same structure of ); use (31) and (38) to obtain , where . Also, note that . Then, substitute the expressions for and into the second term of the right-hand side of (43). This yields


where is the identity () matrix and . Using Voigt notation and observing that (the identity (33) matrix), in (3.5) we have used that , where is the constitutive matrix for an isotropic linear elastic material given by


for plane strain condition, and


for plane stress condition, where is the Young’s modulus and is the Poisson’s ratio.

The first term on the right-hand side of (3.5) is the consistency part of the discrete VEM element bilinear form that provides patch test satisfaction when the solution is a linear displacement field (condition (44) is satisfied). The second term on the right-hand side of (3.5) is the stability part of the discrete VEM element bilinear form and is dependent on the matrix . This matrix must be chosen such that condition (45) holds without putting at risk condition (44) already taken care of by the consistency part. There are quite a few possibilities for this matrix (see for instance BeiraoDaVeiga-Brezzi-Cangiani-Manzini-Marini-Russo:2013 ; BeiraodaVeiga-Brezzi-Marini:2013 ; Gain-Talischi-Paulino:2014 ). Herein, we adopt given by Gain-Talischi-Paulino:2014


where is the scaling parameter and is typically set to 1.

From (3.5), the final expression for the VEM element stiffness matrix is obtained respectively as the summation of the element consistency and stability stiffness matrices, as follows:


where we recall that . Note that and , which are given in (35) and (40), respectively, are easily computed using the nodal coordinates of the element. However, in order to compute and (see their expressions in (36) and (41), respectively), we need some knowledge of the basis functions so that and can be determined. Observe that is computed using (9), which requires the knowledge of the basis functions at the element nodes. And is computed using (37), which requires the knowledge of the basis functions on the element edges. Hence, everything we need to know about the basis functions is their behavior on the element boundary.

We have already mentioned that the basis functions in the VEM are assumed to be Lagrange-type functions. This provides everything we need to know about them on the boundary of an element: basis functions are piecewise linear (edge by edge) and continuous on the element edges, and have the Kronecker delta property. Therefore, can be computed simply as


and can be computed exactly using a trapezoidal rule, which gives


where are the components of and is the length of the edge incident to node as defined in Fig. 1.

The adoption of (51) and (52) in the VEM, results in an algebraic evaluation of the element stiffness matrix. This also means that the basis functions are not evaluated explicitly — in fact, they are never computed. Thus, basis functions are said to be virtual. In addition, the knowledge of the basis functions in the interior of the element is not required, although the linear approximation of the displacement field everywhere in the element is computable through the projection (21). Therefore, a more specific discrete global trial space than the one already given in Section 2 can be built by assembling element by element the local space defined as BeiraoDaVeiga-Brezzi-Cangiani-Manzini-Marini-Russo:2013 ; BeiraodaVeiga-Lovadina-Mora:2015


3.6 VEM element body and traction force vectors

For linear displacements, the body force can be approximated by a piecewise constant. Typically, this piecewise constant approximation is defined as the cell-average . Thus, the body force part of the discrete VEM element linear form can be computed as follows BeiraoDaVeiga-Brezzi-Cangiani-Manzini-Marini-Russo:2013 ; BeiraodaVeiga-Brezzi-Marini:2013 ; artioli:AO2DVEM:2017 :




Hence, the VEM element body force vector is given by


The traction force part of the VEM element linear form is similar to the integral expression given in (54), but the integral is one dimension lower. Therefore, on considering the element edge as a two-node one-dimensional element, the VEM element traction force vector can be computed on an element edge lying on the natural (Neumann) boundary, as follows:




and .

3.7 -norm and -seminorm of the error

To assess the accuracy and convergence of the VEM, two global error measures are used. The relative -norm of the displacement error defined as


and the relative -seminorm of the displacement error given by


where the strain appears in Voigt notation and (see (3) and (20)).

3.8 VEM element stiffness matrix for the Poisson problem

The VEM formulation for the Poisson problem is derived similarly to the VEM formulation for the linear elastostatic problem. However, herein we give the VEM stiffness matrix for the Poisson problem by reducing the solution dimension in the two-dimensional linear elastostatic VEM formulation. The following reductions are used: the displacement field reduces to the scalar field , the strain is simplified to , the rotations become , and the constitutive matrix is replaced by the identity () matrix. Hence, the VEM projections for the Poisson problem become and . The matrices that result from the discretization of the projection operators are simplified to


On using the preceding matrices, the projection matrix is and the final expression for the VEM element stiffness matrix is written as


where is the identity () matrix and has been used in the stability stiffness as this represents a suitable choice for in the Poisson problem BeiraoDaVeiga-Brezzi-Cangiani-Manzini-Marini-Russo:2013 .

4 Object-oriented implementation of VEM in C++

In this section, we introduce Veamy, a library that implements the VEM for the linear elastostatic and Poisson problems in two dimensions using object-oriented programming in C++. For the purpose of comparison with the VEM, a module implementing the standard FEM is available within Veamy for the solution of the two-dimensional linear elastostatic problem using three-node triangular finite elements. In Veamy, entities such as element, degree of freedom, constraint, among others, are represented by C++ classes.

Veamy uses the following external libraries:

  • Triangle shewchuk96b , a two-dimensional quality mesh generator and Delaunay triangulator.

  • Clipper clipperweb , an open source freeware library for clipping and offsetting lines and polygons.

  • Eigen eigenweb , a C++ template library for linear algebra.

Triangle and Clipper are used in the implementation of Delynoi delynoiweb , a polygonal mesh generator that is based on the constrained Voronoi diagram. The usage of our polygonal mesh generator is covered in Section 5.

Veamy is free and open source software and is available from Netlib ( as the na51 package. In addition, a website ( is available, where the software is maintained. After downloading and uncompressing the software, the main directory “Veamy-2.1/” is created. This directory is organized as follows. The source code that implements the VEM is provided in the folder “veamy/” and the subfolders therein. External libraries that are used by Veamy are provided in the folder “lib/.” The folder “matplots/” contains MATLAB functions that are useful for plotting meshes and the VEM solution, and for writing a PolyMesher Talischi:POLYM:2012 mesh and boundary conditions to a text file that is readable by Veamy. A detailed software documentation with graphical content can be found in the tutorial manual that is provided in the folder “docs/.” Several tests are located in the folder “test/.” Some of these tests are covered in the tutorial manual and in Section 6 of this paper. Veamy supports Linux and Mac OS machines only and compiles with g++, the GNU C++ compiler (GCC 7.3 or newer versions should be used). The installation procedure and the content that comprises the software are described in detail in the README.txt file (and also in the tutorial manual), which can be found in the main directory.

The core design of Veamy is presented in three UML diagrams that are intended to explain the numerical methods implemented (Fig. 2), the problem conditions inherent to the linear elastostatic and Poisson problems (Fig. 3), and the computation of the -norm and -seminorm of the errors (Fig. 4).

4.1 Numerical methods

The Veamy library is divided into two modules, one that implements the VEM and another one that implements the FEM. Fig. 2 summarizes the implementation of these methods. Two abstract classes are central to the Veamy library, Calculator2D and Element. Calculator2D is designed in the spirit of the controller design pattern. It receives the ProblemDiscretization subclasses with all their associated problem conditions, creates the required structures, applies the boundary conditions and runs the simulation. Calculator2D, as an abstract class, has a number of methods that all inherited classes must implement. The two most important are the one in charge of creating the elements, and the one in charge of computing the element stiffness matrix and the element (body and traction) force vector. We implement two concrete Calculator2D classes, called Veamer and Feamer, with the former representing the controller for the VEM and the latter for the FEM.

Figure 2: UML diagram for the Veamy library. VEM and FEM modules

On the other hand, Element is the class that encapsulates the behavior of each element in the domain. It is in charge of keeping the degrees of freedom of the element and its associated stiffness matrix and force vector. Element contains methods to create and assign degrees of freedom, assemble the element stiffness matrix and the element force vector into the global ones. An Element has the information of its defining polygon (the three-node triangle is the lowest-order polygon) along with its degrees of freedom. Element has two inherited classes, VeamyElement and FeamyElement, which represent elements of the VEM and FEM, respectively. They are in charge of the computation of the element stiffness matrix and the element force vector. Algorithm 1 summarizes the implementation of the linear elastostatic VEM element stiffness matrix in the VeamyElement class using the notation presented in Sections 3.4 and 3.5.

, , , for each node in the polygonal element do
      Get incident edges Compute the unit outward normal vector to each incident edge Compute and , and insert them into and , respectively Compute and , and insert them into and , respectively
end for
Compute , , , , Compute Output:
Algorithm 1 Implementation of the VEM element stiffness matrix for the linear elastostatic problem in the VeamyElement class

The element force vector is computed with the aid of the abstract classes BodyForceVector and TractionVector. Each of them has two concrete subclasses named VeamyBodyForceVector and FeamyBodyForceVector, and VeamyTractionVector and FeamyTractionVector, respectively.

Even though we have implemented the three-node triangular finite element only as a means to comparison with the VEM, we decided to define FeamyElement as an abstract class so that more advanced elements can be implemented if desired. Finally, each FeamyElement concrete implementation has a ShapeFunction

concrete subclass, representing the shape functions that are used to interpolate the solution inside the element. For the three-node triangular finite element, we include the

Tri3ShapeFunctions class.

One of the structures related to all Element classes is called DOF. It describes a single degree of freedom. The degree of freedom is associated with the nodal points of the mesh according to the ProblemDiscretization subclasses. So, in the linear elastostatic problem each nodal point has two associated DOF instances and in the Poisson problem just one DOF instance. The DOF instances are kept in a list inside a container class called DOFS.

Although the VEM matrices are computed algebraically, the FEM matrices in general require numerical integration both inside the element (area integration) and on the edges that lie on the natural boundary (line integration). Thus, we have implemented two classes, AreaIntegrator and LineIntegrator, which contain methods that integrate a given function inside the element and on its boundary. There are several classes related to the numerical integration. IntegrableFunction is a template interface that has a method called apply that must be implemented. This method receives a sample point and must be implemented so that it returns the evaluation of a function at the sample point. We include three concrete IntegrableFunction implementations, one for the body force, another one for the stiffness matrix and the last one for the boundary vector.

4.2 Problem conditions

Fig. 3 presents the classes for the problem conditions used in the linear elastostatic and Poisson problems. The problem conditions are kept in a structure called Conditions that contains the physical properties of the material (Material class), the boundary conditions and the body force. BodyForce is a class that contains two functions pointers that represent the body force in each of the two axes of the Cartesian coordinate system. These two functions must be instantiated by the user to include a body force in the problem. By default, Conditions creates an instance of the None class, which is a subclass of BodyForce that represents the nonexistence of body forces. Material is an abstract class that keeps the elastic constants associated with the material properties (Young’s modulus and Poisson’s ratio) and has an abstract function that computes the material matrix; Material has two subclasses, MaterialPlaneStress and MaterialPlaneStrain, which return the material matrix for the plane stress and plane strain states, respectively.

Figure 3: UML diagram for the Veamy library. Problem conditions

To model the boundary conditions, we have created a number of classes: Constraint is an abstract class that represents a single constraint — a constraint can be an essential (Dirichlet) boundary condition or a natural (Neumann) boundary condition. PointConstraint and SegmentConstraint are concrete classes implementing Constraint and representing a constraint at a point and on a segment of the domain, respectively. Constraints is the class that manages all the constraints in the system, and the relationship between them and the degrees of freedom; EssentialConstraints and NaturalConstraints inherit from Constraints. Finally, ConstraintsContainers is a utility class that contains EssentialConstraints and NaturalConstraints instances. Constraint keeps a list of domain segments subjected to a given condition, the value of this condition, and a certain direction (vertical, horizontal or both). The interface called ConstraintValue is the method to control the way the user inputs the constraints: to add any constraint, the user must choose between a constant value (Constant class) and a function (Function class), or implement a new class inheriting from ConstraintValue.

4.3 Norms of the error

As shown in Fig. 4, Veamy provides functionalities for computing the relative -norm and -seminorm of the error through the classes L2NormCalculator and H1NormCalculator, respectively, which inherit from the abstract class NormCalculator. Each NormCalculator instance has two instances of what we call the NormIntegrator classes: VeamyIntegrator and FeamyIntegrator. These are in charge of integrating the norms integrals in the VEM and FEM approaches, respectively. In these NormIntegrator classes, the integrands of the norms integrals are represented by the Computable class. Depending on the integrand, we define various Computable subclasses: DisplacementComputable, DisplacementDifferenceComputable, H1Computable and its subclasses, StrainDifferenceComputable, StrainStressDifferenceComputable, StrainComputable and StrainStressComputable. Finally, DisplacementCalculator and StrainCalculator (and their subclasses) permit to obtain the numerical displacement and the numerical strain, respectively; and StrainValue and StressValue classes represent the exact value of the strains and stresses at the quadrature points, respectively.

Figure 4: UML diagram for the Veamy library. Computation of the -norm and -seminorm of the error

4.4 Computation of nodal displacements

Each simulation is represented by a single Calculator2D instance, which is in charge of conducting the simulation through its simulate method until the displacement solution is obtained. The procedure is similar to a finite element simulation. The implementation of the simulate method is summarized in Algorithm 2.

Input: Mesh
Initialization of the global stiffness matrix and the global force vector for each element in the mesh do
      Compute the element stiffness matrix Compute the element force vector Assemble the element stiffness matrix and the element force vector into global ones
end for
Apply natural boundary conditions to the global force vector Impose the essential boundary conditions into the global matrix system Solve the resulting global matrix system of linear equations Output: Column vector containing the nodal displacements solution
Algorithm 2 Implementation of the simulate method in the Calculator2D class

The resulting matrix system of linear equations is solved using appropriate solvers available in the Eigen library eigenweb for linear algebra.

5 Polygonal mesh generator

In this section, we provide some guidelines for the usage of our polygonal mesh generator Delynoi delynoiweb .

5.1 Domain definition

The domain is defined by creating its boundary from a counterclockwise list of points. Some examples of domains created in Delynoi are shown in Fig. 5. We include the possibility of adding internal or intersecting holes to the domain as additional objects that are independent of the domain boundary. Some examples of domains created in Delynoi with one and several intersecting holes are shown in Fig. 6.

Figure 5: Domain examples. (a) Square domain, (b) rhomboid domain, (c) quarter circle domain, (d) unicorn-shaped domain

Listing 1 shows the code to generate a square domain and a quarter circle domain. More domain definitions are given in Section 6 as part of Veamy’s sample usage problems.

1std::vector<Point> square_points = {Point(0,0), Point(10,0), Point(10,10), Point(0,10)};
2Region square(square_points);
3std::vector<Point> qc_points = {Point(0,0), Point(10,0), Point(10,10)};
4std::vector<Point> quarter = delynoi_utilities::generateArcPoints(Point(10,0), 10, 90.0, 180.0);
5qc_points.insert(quarter_circle_points.end(), quarter.begin(), quarter.end());
6Region quarter_circle(qc_points);
Listing 1: Definition of square and quarter circle domains
Figure 6: Examples of domains with holes. (a) Square with an inner hole, (b) square with four intersecting holes, (c) unicorn-shaped domain with an inner hole, (d) unicorn-shaped domain with an intersecting hole

To add a circular hole to the center of the square domain already defined, first the required hole is created and then added to the domain as shown in Listing 2.

1Hole circular = CircularHole(Point(5,5), 2);
Listing 2: Adding a circular hole to the center of the square domain

5.2 Mesh generation rules

We include a number of different rules for the generation of the seeds points for the Voronoi diagram. These rules are constant, random_double, ConstantAlternating and sine. The constant

method generates uniformly distributed seeds points; the

random_double method generates random seeds points; the ConstantAlternating method generates seeds points by displacing alternating the points along one Cartesian axis. Fig. 7 presents some examples of meshes generated on a square domain using different rules. We show how to generate constant (uniform) and random points for a given domain in Listing 3.

Figure 7: Polygonal mesh generation on a square domain using different rules. (a) constant, (b) random_double, (c) ConstantAlternating, (d) sine
1dom1.generateSeedPoints(PointGenerator(functions::constant(), functions::constant()), nX, nY);
2dom2.generateSeedPoints(PointGenerator(functions::random_double(0,maxX), functions::random_double(0,maxY)), nX, nY);