1 Introduction
When the Galerkin weak formulation of a boundaryvalue 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 nonpolynomials or highorder 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 TalischiPaulino: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 BeiraoDaVeigaBrezziCangianiManziniMariniRusso: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 finitedifference 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 GainTalischiPaulino:2014 ; ManziniRussoSukumar:2014 ; BeiraodaVeigaLovadinaMora: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 highorder 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, objectoriented 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 builtin 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. GainTalischiPaulino: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 wellestablished free and open source objectoriented 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 50line 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 builtin 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 MATLABbased polygonal mesh generator.

It is designed using the objectoriented 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 twodimensional linear elastostatic problem and the twodimensional Poisson problem.
We are also aware of the MATLAB Reservoir Simulation Toolbox lie:MRS:2016 , which provides a module for first and secondorder virtual element methods for Poissontype 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 objectoriented C++ implementation of the VEM.
The remainder of this paper is structured as follows. The model problem for twodimensional linear elastostatics is presented in Section 2. Section 3 summarizes the theoretical framework of the VEM for the twodimensional linear elastostatic problem. Also in this section, the VEM element stiffness matrix for the twodimensional Poisson problem is given. The objectoriented implementation of Veamy is described and explained in Section 4. In Section 5, some guidelines for the usage of Veamy’s builtin 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 onedimensional 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:
(1) 
where
is the Cauchy stress tensor and
is the gradient operator. The gradient of the displacement field can be decomposed into its symmetric () and skewsymmetric (
) parts, as follows:(2) 
where
(3) 
is the strain tensor, and
(4) 
is the skewsymmetric gradient tensor that represents rotations. The Cauchy stress tensor is related to the strain tensor by
(5) 
where is a fourthorder 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:
(6) 
which leads to the standard form of presenting the weak formulation: find such that
(7a)  
(7b)  
(7c) 
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 closedform 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 meshdependent. 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 twodimensional 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.
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
(8) 
where is defined through the mean value of a function over the element nodes given by
(9) 
where is the number of nodes of coordinates that define the polygonal element^{1}^{1}1Eq. (9) in fact defines any ‘barred’ term that appears in this paper.; is a secondorder tensor and thus can be uniquely expressed as the sum of a symmetric and a skewsymmetric tensor. Let the symmetric and skewsymmetric tensors be denoted by and , respectively. The spaces of rigid body modes and constant strain states over are defined, respectively, as follows:
(10) 
(11) 
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:
(12) 
for extracting the rigid body modes,
(13) 
for extracting the constant strain states, and
(14) 
for extracting the linear polynomial part. And since , the projection operators satisfy the relation
(15) 
We know by definition that the space includes linear displacements. This means that . Thus, any can be decomposed into three terms, as follows:
(16a)  
(16b) 
that is, into their rigid body modes, their constant strain states and their additional nonpolynomial or highorder functions, respectively.
The explicit forms of the projection operators that are defined through (12)(14) are given in Ref. GainTalischiPaulino:2014 and are summarized as follows: let the cellaverage of the strain tensor be defined as
(17) 
where the divergence theorem has been used to transform the volume integral into a surface integral. Similarly, the cellaverage of the skewsymmetric gradient tensor is defined as
(18) 
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
(19) 
where and are the rotation and translation modes of , respectively. And the projection of onto the space of constant strain states is given by
(20) 
Hence, by (15) the projection of onto the space of linear displacements is written as
(21) 
The projection operator satisfies some important energyorthogonality conditions that are invoked when constructing the VEM bilinear form. The proofs can be found in Ref. GainTalischiPaulino:2014 . The energyorthogonality conditions are given next.
The projection satisfies:
(22a)  
(22b) 
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:
(23) 
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 energyorthogonality condition (22b) have been used.
The first term on the righthand 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 nonpolynomial or highorder 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 twodimensional Cartesian space and the skewsymmetry of ^{3}^{3}3Note that and .. The projection (19) can be written as follows:
(24) 
where the basis for the space of rigid body modes is:
(25) 
Similarly, on considering the symmetry of , the projection (20) can be written as
(26) 
where the basis for the space of constant strain states is:
(27) 
On each polygonal element of edges with nodal coordinates denoted by , the trial and test displacements are locally approximated as
(28) 
where and are assumed to be the canonical basis functions having the Kronecker delta property (i.e., Lagrangetype 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:
(29) 
and the components of the basis for the space of constant strain states:
(30) 
The discrete version of the projection to extract the rigid body modes is obtained by substituting (28) and (29) into (24), which yields
(31) 
where
(32) 
(33) 
and
(34) 
with
(35) 
and
(36) 
Similarly, substituting (28) and (30) into (26) leads to the following discrete version of the projection to extract the constant strain states:
(38) 
where
(39) 
with
(40) 
and
(41) 
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:
(42) 
3.5 VEM element stiffness matrix
The decomposition given in (3.3) is used to construct the approximate meshdependent 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
(43) 
where its righthand 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 BeiraoDaVeigaBrezziCangianiManziniMariniRusso:2013 ; BeiraodaVeigaBrezziMarini:2013 :
For all and for all in

Consistency: and
(44) 
Stability: two constants and , independent of and of , such that
(45)
The discrete version of the VEM element bilinear form (43) is constructed as follows. Substitute (3.4) into the first term of the righthand 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 righthand side of (43). This yields
(46) 
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
(47) 
for plane strain condition, and
(48) 
for plane stress condition, where is the Young’s modulus and is the Poisson’s ratio.
The first term on the righthand 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 righthand 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 BeiraoDaVeigaBrezziCangianiManziniMariniRusso:2013 ; BeiraodaVeigaBrezziMarini:2013 ; GainTalischiPaulino:2014 ). Herein, we adopt given by GainTalischiPaulino:2014
(49) 
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:
(50) 
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 Lagrangetype 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
(51) 
and can be computed exactly using a trapezoidal rule, which gives
(52) 
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 BeiraoDaVeigaBrezziCangianiManziniMariniRusso:2013 ; BeiraodaVeigaLovadinaMora:2015
(53) 
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 cellaverage . Thus, the body force part of the discrete VEM element linear form can be computed as follows BeiraoDaVeigaBrezziCangianiManziniMariniRusso:2013 ; BeiraodaVeigaBrezziMarini:2013 ; artioli:AO2DVEM:2017 :
(54) 
where
(55) 
Hence, the VEM element body force vector is given by
(56) 
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 twonode onedimensional element, the VEM element traction force vector can be computed on an element edge lying on the natural (Neumann) boundary, as follows:
(57) 
where
(58) 
and .
3.7 norm and seminorm of the error
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 twodimensional 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
(61) 
(62) 
(63) 
(64) 
On using the preceding matrices, the projection matrix is and the final expression for the VEM element stiffness matrix is written as
(65) 
where is the identity () matrix and has been used in the stability stiffness as this represents a suitable choice for in the Poisson problem BeiraoDaVeigaBrezziCangianiManziniMariniRusso:2013 .
4 Objectoriented 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 objectoriented 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 twodimensional linear elastostatic problem using threenode 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 twodimensional 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 (http://www.netlib.org/numeralgo/) as the na51 package. In addition, a website (http://camlab.cl/software/veamy/) is available, where the software is maintained. After downloading and uncompressing the software, the main directory “Veamy2.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.
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 threenode triangle is the lowestorder 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.
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 threenode 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 threenode 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.
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.
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.
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.
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.
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.
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.