An open source hp -adaptive discontinuous Galerkin finite element solver for linear elasticity Engineering Software

Open source codes are a key ingredient to greater research integrity and accountability in computational science and engineering. However, many of these codes have not been developed with modification of the base code as their primary consideration. Existing codes may provide an environment for researchers to quickly test out their ideas under different physical conditions in a high level way but they are not always ideal for those interested in the development of numerical methods. The majority of existing open source discontinuous Galerkin finite element codes are written in C ++ and there is a significant learning curve for junior researchers to adopt, understand and modify the underlying code/routines. This paper presents an open source hp -adaptive discontin- uous Galerkin finite element code written in MATLAB that has been explicitly designed to make it easy for users, especially MSc/PhD-level researchers, to understand the method and implement new ideas within the core code. Although the code is focused on solving problems in linear elasticity, it is straightforward to modify it to solve other physical equations.


Introduction
Open source codes are essential for the health of research in computational science and engineering as they provide a route-in for researchers to start working with numerical methods without having to implement hundreds (often tens of thousands) of lines of computer code. Open source code also promotes collaboration between researchers and transparent inspection of numerical results -resulting in greater research integrity and accountability [1]. There are already available many open source codes to solve specific engineering problems like for example fractures in materials [2,3] or simulating structures using plates and shells [4,5].
The majority of numerical methods for solving partial differential equations (PDEs) can be divided into two categories; those based on meshes partitioning the physical domain and those not using meshes. Meshless methods [6] were invented to mitigate some of the difficulties associated with resolving the geometry of physical domains with meshes. In common with methods based on meshes like finite element methods (FEMs), meshless methods use shape functions to approximate the solutions but not as interpolants as in FEMs, rather than as approximants. Partition of unity methods is another class of methods not based on standard meshes but on a coverage of the physical domain with overlapping patches [7,8]. The category of methods based on meshes is dominated by FEMs and their variations. One of the most prominent classes of FEMs is multiscale methods [9,10] in which the approximated solution is decomposed into coarse and fine-scale solutions. This is achieved by rewriting the PDE model into a pair of coupled macroscopic and microscopic models. Other classes of methods were invented to overcome some of the aspects of classical FEMs. Isogeometric analysis breaks free from the limitations of polygonal and polyhedral shaped elements and polynomials of integer order shape functions [11]. Extended and enriched FEMs allow adding to the approximation spaces functions incorporating specific knowledge of the solution [12,13]. Cut FEMs eliminate the difficulties to fit complex or implicit surfaces with standard meshes [14,15]. A lot of effort has been made in studying ways to control the error in FEMs [16][17][18]. A very popular approach is adaptivity as a way to improve the accuracy of computed solutions in many fields adapting meshes without leaving the FEM framework [19][20][21].
There are of course methods not easily classifiable in these two categories like numerical methods based on Taylor-series [22].
However, the solution of partial differential equations in engineering analysis is dominated by the continuous Galerkin (CG) variant of FEM. There are a huge number of open source codes based on conventional CG finite elements (for example, see [23][24][25][26][27][28], amongst many others). An alternative to the conventional CG formulation is the family of discontinuous Galerkin (DG) finite element methods and there are many reasons that researchers may wish to adopt a DG modelling framework, such as: (i) their ability to easily allow for adaptive mesh refinement, hanging nodes and different polynomial orders in adjacent elements due to the relaxation of the compatibility between elements; (ii) that, unlike CG elements, DG elements can handle excessive aspect ratios without degradation of their predictive capabilities; (iii) that, unlike CG methods, degrees of freedom are not shared between elements, making easier to split a mesh on several computing nodes; (iv) that, unlike CG methods, boundary conditions are imposed weakly, making easier to incorporate complicated boundary conditions into models; (v) more freedom in choosing shape functions since no continuity across elements is needed.
There are several variants of DG finite element methods, categorised by the inter-element communication technique used to couple adjacent elements. In this paper we adopt a symmetric interior penalty DG (SIPG) approach. Despite their advantages, there are far fewer DG finite element open source codes (compared to CG-based methods) and the majority of these codes have not been developed with understanding and modification of their core code as their primary focus. To the best of the authors' knowledge, the following open-source DG codes are freely available (grouped by underpinning implementation language): C++ based: hpGEM [29], DUNE (Python binding) [30], DoGPack (Python/MATLAB post-processing possible) [31], NEKTAR++, FreeFEM [32], MFEM [33] and FEniCS (Python front-end available) [34]; Python based: PyFR [35] and Nutils [36]. It is clear that the majority of these codes are based on C++ for computational efficiency. However, several researchers are interested in the development of methods rather than the analysis of large-scale problems. For these researchers, being able to easily access and understand the underlying source code is essential and these C++ implementations may be off-putting or a complete barrier to using these open source codes to develop their ideas. Python, being open-source to its very core, is a credible alternative and several of the C++ codes have Python front ends and/or processing environments, as identified above. However, in the authors' opinion, within mathematical and technical computing, MATLAB is the dominant tool taught to engineers and scientists whilst at university and MATLAB has been specifically designed with these users in mind [39]. The documentation, help and debugging environment offered by MATLAB is unparalleled and this allows researchers (especially junior scientists/engineers) to quickly develop their scientific ideas whilst minimising the syntax barrier.
In this paper we present a MATLAB-based open-source symmetric interior penalty DG finite element code that has been designed to make the method accessible to early career researchers, such as final year undergraduates, MSc and PhD students as well as post-doctoral researchers moving into this area of research. This point is emphasised by the fact that the first author of this manuscript was a final year MEng student during the development of this paper and associated code. Another point of novelty of this paper is the inclusion of error estimate driven hp-adaptivity, where the code automatically refines the background mesh in terms of element size (h) or polynomial order (p) of specific elements based on the estimate error in each element. This provides a powerful stress analysis tool which is able to achieve extreme accuracy (errors less than 10 − 10 ) with modest computational effort and on resources that are commonly available to early career researchers (standard laptop/desktop PCs). The physical problem focus of the code is linear elasticity, however once the structure of the code is understood, it is straightforward to modify the code to analyse other physical problems, include alternative element formulations, test out different adaptivity routines, etc.
Following this introduction, the layout of the rest of the paper is as follows: Section 2 outlines the considered model problem and key details of the symmetric interior penalty DG finite element method (including the error estimate and adaptivity algorithm), Section 3 explains how to install and test the method, Section 4 details the structure of the program and demonstration problems are show in Section 5. Finally, brief conclusions are included in Section 6. All of the code and associated documentation is available from https://github.com/Robert-Bird/Disco ntinuous-Galerkin-MATLAB.

The method
In this section, we detail the numerical analysis aspects of the implemented method. The model problem and numerical methods implemented in the code are analysed in [40] where in-depth analysis can be found.

Model problem
The problem implemented in the code is linear elasticity with a variety of boundary conditions: Let Ω be a bounded polygonal domain in where n = (n x , n y ) is the unit vector perpendicular to the boundary of Ω and pointing out and n ‖ is the tangential unit vectors to the boundary, t(u) is the traction component of the stress, i.e. The set Γ D may be empty. In this case, problem (1) may not have a unique solution if rigid motions are in the kernel [41]. Therefore, the natural choice of space to enforce a unique solution for problem (1) is where R is the space containing all rigid motions. In case the prescribed boundary conditions define a problem with rigid motions in the kernel, the well-definiteness of the problem can be enforced using the average boundary condition to remove rigid motions [42]. Such a case happens for example when only Neumann boundary conditions are used, see problem 3 in Section 3.2. The code automatically applies the average boundary condition if needed.

Symmetric interior penalty discontinuous Galerkin method
In this section we introduce our DG method used in the code to solve problem (1). We assume that the computational domain Ω can be partitioned into a shape regular mesh T of triangular and affine elements and we denote with K a generic element of T . We allow for irregular meshes with a maximum of one hanging node per edge. Due to our assumption that the initial mesh is shape regular, allowing at most one hanging node per edge implies that the diameters of the elements in all refined meshes are of bounded variation for any pair of neighbouring elements. We denote E (T ) and E int (T )⊂E (T ) the set of all edges of the mesh T and the subset of all interior edges respectively and by E BC (T )⊂E (T ) the subset of all boundary edges. The set E BC (T ) is partitioned in the three subsets E D (T ), E N (T ) and E T (T ) that are the sets containing the edges forming the three portions of the boundary Γ D , Γ N and Γ T . We define h K and h E to be the diameter of the element K and the length of the edge E respectively.
In the implementation, we allow to vary the polynomial degrees across elements under the constraints that the variation is at most equal to one for any pair of neighbouring elements. This is enforced during mesh adaptation, too. For each element K of the mesh T we associate a polynomial degree p K ≥ 1 and we introduce the degree vector p = { p K : K ∈ T } and we define p min as the minimum of p K on the mesh T . For any E ∈ E (T ), we introduce the edge polynomial degree p E by where + and -denote the two elements sharing the face E. Hence, for a given partition T of Ω and a degree vector p on T , we define the hp-version DG finite element space by where P pK (K) is the space of polynomials of degree at most p K and R is the set of rigid motions. We also use + and -to denote all quantities related to the elements. Given an edge E ∈ E int (T ) shared by two elements K + and K − , we define the jump [[⋅]] operator and the average {⋅} operator on vectors and tensors as: where n + K = (n + x , n + y ) the outward unit normal on the boundary ∂K + of an element K. Such definitions are modified on edges along the boundary of the domain, i.e. E ∈ E BC (T )⊂E (T ), along the boundary we set Thus, the DG approximation problem (1) reads as follows: Find u h ∈ V p (T ) such that where the bilinear form where γ is the penalty constant. As proved in [43], the penalty term adjust for any variation in h or p in the mesh. This means that the penalty constant γ is mesh independent and the only condition that γ must satisfy is to be big enough as stated in [43]. In the authors' experience, γ = 10 is big enough for most problems. For more complicated PDEs, the definition of the penalty term can be modified in view of [44] to take into account varying coefficients in the second order term of the PDE, see for example [45][46][47][48][49][50]. The natural norm for problem (5) is the DG norm: where ‖ ⋅ ‖ 0,K and ‖ ⋅ ‖ 0,E are respectively the L 2 -norm on an element K and on an edge E.

Error estimator
This section defines the error estimator implemented in the code. In [40] the reliability and efficiency proofs of the error estimator are presented.
The error estimator is defined as where the three terms under the sum are defined as where f h and g N,h are the L 2 projections of f and g N onto the finite element space and where g D,h and g T,h are piecewise polynomial approximated of traces of functions in H 1 (Ω) as in [51].

Adaptivity
The code supports arbitrary high order triangular elements as defined in Section 2.2.3 in [52]. The hp-adaptive strategy used here was originally proposed in [53] and consequently used in [54]. The elements are chosen for either h or p refinement using Algorithm 1. The marking strategy in Algorithm 1 uses two threshold values δ 1 and δ 2 , with δ 2 ≥ δ 1 , to determine what elements to refine in h and what elements in p. All the elements satisfying η 2 K , are marked for h − refinement and all elements satisfying The remaining elements are not marked for refinement at all. To only apply h-adaptivity set δ 2 = δ 1 ∕ = 0. Similarly, to only apply p-adaptivity, set δ 2 = 1 and δ 1 ∕ = 0. To uniformly refine in h set δ 2 = δ 1 = 0 and to uniformly refine in p set δ 2 = 1 and δ 1 = 0.

Installation
The only MATLAB add on that SIPG Linear Elastic uses is the Symbolic Math Toolbox™, which can be downloaded from mathworks.com. To begin the installation, clone or download the repository from GitHub. Discontinuous-Galerkin-MATLAB is a MATLAB directory environment that contains main functions, placed in the top directory, that can be run directly, the path_add.m script to add folders to the MATLAB search path, and folders containing routines that perform the analysis. The package has the functionality to support mesh generation with the MESH2D, a MATLAB based Delauney mesh generator for two-dimensional geometries [55,56] which can be downloaded from the MathWorks® website. If the user wishes, this package can be integrated into the program by copying the MESH2D file into the SIPG Linear Elastic directory. Care should be taken to ensure that the folder containing the package is added to the MATLAB search path, so it is advised that the user edits line 28 of path_add.m as appropriate. The three problems described in Section 5 can be run without MESH2D. The code analytical_problem.m needs MESH2D. Table 1 describes the subfolders in SIPG_linear_elastic_v1.0.
The reference manual can be generated automatically using the M2HTML toolbox [57]: m2html('mfiles', 'SIPG_dir', 'htmldir','doc_dir', 'recursive','on'); where SIPG_dir is the relative path to the folder Discontinuous-Galerkin-MATLAB and doc_dir is the relative path to the folder where the user wants the reference manual to be saved. All functions in the Discontinuous-Galerkin-MATLAB code are integrated with the MATLAB help system. The installation consists in running the script path_add.m to add the necessary folders to the MATLAB search path.

Testing
After installation, the user should call the function tests.m in the command window. The function tests.m is based on the unit test framework used in Matlab. This is used to check that the program is working as expected on the user's computer by performing a set of four simulations, listed in Table 2, and thereby testing critical routines. By comparing the values from this simulation with a set of expected values stored in .txt files, the function checks that all routines are performing nominally.
Apart from the global stiffness matrix calculation DG_algorithm.m and the volumetric integral of the body force force_integration_vol.m, other routines common to all test problems are the error estimation algorithms error_calc_ele.m and the calculation of the L 2 and DG norms. Depending on the boundary conditions defined, different subroutines within the error estimation and the calculation of the DG norm will be executed and tested by tests.m.
The primary role of test problem 1 is to verify that the mesh adaptivity algorithms are functioning correctly. To achieve this, the program solves the problem with the displacement solution shown in Table 2 with five different mesh adaptivity methods: uniform h, uniform p, adaptive h, adaptive p and adaptive hp. Test problems 2, 3, and 4 establish that the force surface integral functions associated with each type of boundary condition are performing as expected. When executing a problem with only Neumann boundary conditions, the problem is indeterminate and rigid body motion is possible. To prevent rigid body motion, such that all displacements and rotations are zero, for instance, test problem 3, it is necessary to generate Lagrangian force components to cause the system to act in a deterministic manner. This is performed in the UV_ave_bc.m and rot_ave_BC.m, which are additional algorithms whose performance is checked by test problem 3.

Mesh data structures and flags
We now consider the format used in the program. The initial mesh of all simulations must be a conforming mesh, i.e. no hanging node must be present. The coordinates of the nodes are stored in the variable coord and element topology information is stored in the variable etpl.mat, discussed below.
The elements used are triangular with an anticlockwise local node numbering convention, and the global node index begins at one. Global nodal coordinates are stored in the variable coord, an [x,y] list with dimensions nodes × 2, where nodes is the number of nodes in the mesh. Nodes are shared between elements if they are not hanging nodes. However, in DG methods, degrees of freedom are not shared between elements. Element topology information is stored in the MATLAB structure array etpl, which initially contains the fields mat, poly and tree; an additional field, ed_recalc, is calculated at the start of each adaptivity loop and is not required as an input. A description of each field, their input format and array dimensions is shown in Table 3.
In the mat field, the element faces are defined by stepping Algorithm 1. hp-refinement strategy: anticlockwise around the local nodes, indexed by their global node number (row position in coord). Elements in the mesh are indexed according to their row position in etpl.mat, and are assigned a polynomial order in etpl.poly. As for the first generation of elements every element is active, etpl.tree should be initialised by setting all rows as [1,0,0]; storing information about the history of elements using this variable facilitates adaptation of the code to include mesh derefinement, which may be of interest to advanced users. The final data structure describing the mesh is the element face topology matrix etpl_face, which has the format [element 1,element 2,local face 1, local face 2,normal x,normal y,face type]. In the code it is enforced that element 1 is an element for which the face is a complete face. In case of hanging nodes, element 2 with an hanging node in the middle of the face. This variable stores information about face connectivity for all faces in the mesh, illustrated in Fig. 1.
The local face value ranges between one and three, indexed in an anticlockwise direction around each element; normal x, normal y describe the outward facing normal vector orientation; face type is a flag used to indicate whether the face is internal, and, if it is external, the boundary condition to be imposed. Table 4, below, outlines the flagging convention used in the face type column of etpl_face. If a face is external then the value in the relevant row for element 2 in etpl_face is set to zero, as is the value of local face 2.
The mesh data structure format that SIPG Linear Elastic uses is compatible with the mesh generation algorithm MESH2D. Consequently, if additional information is desired by the user about the data structures discussed in this section, this can be found in the MESH2D documentation.

Creation of data structures and flags
The data structures and their respective fields required to solve the SIPG problem are etpl, coord and etpl_face. They are either hardcoded as. txt files, as in the case of problems 1-to-4 in Table 2 or generated with MESH2D. There are two stages in creating the data structures if MESH2D is used, this is highlighted with the analytical_problem.m example. First, the mesh input data for MESH2D is defined, in this case in analy-tical_problem_generator.m. Second the input data is used by MESH2D to create an element topology which is then restructured into coord, etpl and etpl_face; this happens with the routine seed_mesh_square.m.
The input data required for MESH2D is, node, edge and BC, and in this Once the input data has been defined by the user, seed_mesh_square.m is called. The first operation of which is to call the routine refine2, of MESH2D, which takes node and edge as in inputs and outputs a mesh in the form of an element topology and coordinate matrix, respectively stored in etpl.mat and coord. Additionally, the list conn is also outputted of size [number of external edges × 2], each row is 1 edge and contains 2 global node numbers. It is then extended by the finding_the_boundary.m routine to have a third column that contains the face flags for each edge. finding_the_boundary.m operates on the basis that the node numbers for each external edge in BC also exist in conn. Hence there exists a set of edges in conn that link between the nodes defining a boundary edge in BC, these edges thus all have the same face flags defined in BC. Now that the element topology of the initial conforming mesh has been created, the fields poly and tree of etpl can be defined as in the previous section. etpl.tree is initialised by setting all rows to [1,0,0] and etpl.poly is initialised by setting the rows to [row number,initial polynomial order]. The next stage to seed_mesh_square.m is creating etpl_face with two steps, both of which are called by etpl_face_square_func.m. First, create etpl_face for all internal faces and second, for all external faces. Both are dependent on etpl.mat and coord created by MESH2D, whilst the latter also depends on conn.
The creation of the internal faces is performed by etpl_face_all.m and is described with Algorithm 2. It loops over all the elements in the mesh nels and their local face numbers, 1-to-3, on lines 2 and 3. If the element el and its local face number f have not been stored in etpl_face, controlled by the if statement on line 4, the face counter face_count is increased by 1 and etpl_face_index is updated to indicate the element, el and face f have been seen by the algorithm. Then, el, f and the outward normal n to the face are stored in etpl_face on line 8. Next, a search of etpl.mat is performed to find an neighbour element el_n, this is defined as an element which shares two nodes with el. If el_n exists, then it and the local face number that contains the shared nodes are stored in etpl_face, additionally etpl_face_index is updated to indicate the element, el_n and face f_n have been seen.
Algorithm 2 will find all faces in the domain and store them in etpl_face, however the last column of the etpl_face will not be defined for the external faces. This occurs in etpl_face_ext_find_square.m, described by Algorithm 3, and uses conn to assign the face flag to the external faces of etpl_face.
The last data structure to be created is ed which stores for each element its degrees of freedom and steers the local element stiffness matrices into the global matrix. It is created during the accumulation of  It is important to highlight three aspects of the DOF storage in ed. First, elements which are not active will have only zeros in their row, as controlled by the if statement. Second, the DOF only increases with active elements, this ensures that no rows and columns with zero values will appear in the global stiffness matrix when the local element matrices are steered into the global matrix. Last, ed is reset to zero after each adaptive step to account for different elements that have become active or changed polynomial order. This therefore means it is possible   for the DOF for an element to change during an adaptivity step.

Modification of data structures during adaptivity
During the adaptive step the mesh is adapted through changes in the size and polynomial order of elements in the mesh. When the changes occur, the element topology structure etpl, node coordinate list coord and face connectivity and boundary matrix etpl_face need to be updated. The changes to these datasets occur in mesh_refine_topology_marking.m, the structure of which is provided by Algorithms 5 and 6.
mesh_refine_topology_marking.m is governed by the for loop on line 2 of Algorithm 5, it loops over the possible element ages in the mesh from oldest to youngest; a younger element is the result of more refinements compared to an older element. This ensures that during the refinement process a maximum of one hanging node can exist on a face. The second loop on line 3 loops over the elements that are marked for h-refinement, within this loop new nodes, elements and faces are created. The first stage of the algorithm updates etpl and coord with homogeneous h-refinement on element el_ref. For homogeneous refinement midside nodes are required, these are either created resulting in updates to coord, or retrieved as they may exist already, on line 5 using the coor-d_check_new.m routine 2 . The node numbers for the vertices are then stored in Nn1,Nn2 and Nn3, the midside nodes Nn4,Nn5 and Nn6, which are combined to create four new elements on line 9 where nel is the number of elements in the mesh. The assignment of the node numbers in these way ensures that the local face numbers of the new elements coincide with el_ref. Elements nel+1-to-nel+3 use a combination of vertex nodes and midside nodes, nel+4 is a formed from only midside nodes. Next on line 10 new elements are assigned their parent polynomial order, and the active element tree etpl.tree is updated on line 11. This is followed by an updating the face connectivity matrix for all faces corresponding to el_ref on line 12 with Algorithm 6. The update to etpl_face is split into two steps in Algorithm 6, the first step is defined with the for loop which spans lines 4-15, corresponding to the elements nel+1-to-nel+3 as these elements will have faces which interact with element and face information that already exists. The second step occurs with the for loop on lines 16-18 since the faces of element nel+4 are shared only with the newly created elements and hence new face information has to be created. The first section of the algorithm starts by finding the faces that correspond to el_ref and storing these in etpl_face_current_el. The local faces j=1:3 of el_ref are then search for in etpl_face_current_el, for a j they stored in etpl_face_current_face. If there is only one row in etpl_face_current_face, then the new elements that share face j with el_ref are assigned the face information of etpl_face_-current_face and stored in etpl_face. If there are two rows of etpl_face_-current_face then the face j of el_ref has a hanging node, each new element  2 When using the coord_check_new.m routine it is important to highlight that the node check is performed using element topology rather nodal position on j is then assigned the row of etpl_face_current_face. The row which they are assigned is found by finding which row contains an opposite element that it share two nodes with. Last the face information for nel+4 is found and stored in etpl_face on lines 16-18, this achieved using a similar an algorithm similar to Algorithm 2.
When storing new information in etpl_face it is important to highlight that all the new element numbers go into the first column, and all new local face numbers are stored in the third column. This ensures that the element and face information of the smaller element on the face is always in the same columns.

Description of example codes
In this section, we briefly describe the structure of the code analy-tical_problem.m. The definition of the problem is performed setting the variables in Table 5. Examples of such procedures are: analytical_problem_generator.m, analytical_problem_1_generator.m, analytical_problem_2_generator.m and ana-lytical_problem_3_generator.m.
In the code, two ways to set up the mesh are used. In analytical_problem.m, the routine seed_mesh_square is used which calls MESH2D. In analytical_problem_1_generator.m, analytical_problem_2_generator.m and ana-lytical_problem_3_generator.m, the routines squaremesh, Lmesh and crackmesh are used. Such routines read the mesh from files.

How to set up a custom problem
The easiest way to set up a new problem is to make a copy of ana-lytical_problem.m and modify the code accordingly. The two main steps are: 1. Define the new problem defining all the variable in Table 5. 2. Define the mesh to use either loading it from file or generating using MESH2D, see Section 4.4. The user may also decide to write routines to import meshes saved in existing formats.
Optionally, the postprocessing section of the code may be changed to display the information in the requested format.

Example problems
In this section, we present three example problems that demonstrate the capabilities of the package. These examples serve as a guide to problem set up by showing how the problems outlined in [40] can be analysed. These are a problem with a smooth solution on a unit square domain, a non-smooth problem on an L-shaped domain, and a crack in a plate problem, performed by executing the .m scripts example_problem_1. m, example_problem_2.m and example_problem_3.m respectively. These problems do not use the set up scripts for use with MESH2D discussed above to generate the mesh.

Problem with an analytical solution on unit square domain
The first example, example_problem_1.m, considers a small strain linear elastic problem on (x, y) ∈ Ω = (0, 1) 2 where x and y are in meters, see Fig. 2. The problem acts in plane strain with Young's modulus E Y = 5 2 Pa and a Poisson's ratio ν = 1 4 . As discussed in [40], the right-hand side f of Equation (1) should be selected such that the exact analytical solution is For this problem, the boundary condition is homogeneous Dirichlet (zero displacement) along the whole boundary. Fig. 2 shows the mesh used for this example.
In the file analytical_problem_1_generator.m all parameters to run this example are defined. For example, on line 73, the type of boundary conditions are defined. The value 2 corresponds to homogeneous ]. An overview of the variables used in setting up the script is given in Table 5.
The program determines the analytical body force from the analytical displacement solution, and then imposes the analytical body force across the domain. These tasks are performed on lines 105 through 123 using the MATLAB symbolic toolbox. This is done in three steps: 1. Differentiating the displacement with respect to area to find the strain components across the domain 2. Multiplying the strain vector by the elastic stiffness matrix D e to find the stress across the body 3. Differentiating the stress with respect to area to find an expression for the body force To perform the simulations, the user should enter the command example_problem_1. The program will output a table detailing the progress of the mesh generation algorithm, for each adaptivity loop information about the time taken to calculate important parts of the DG algorithm, and the error estimate in the command window. Additionally, for each simulation the program will produce three figures: the first shows a colour coded plot of the boundary conditions applied to the external faces, the second is a log-log plot of the error estimate and the DG and L 2 norm convergence, and the third is a surface plot of the von Mises stress across the domain.
An important difference between SIPG Linear Elastic and the program used in [40] is that the former initialises the simulation with an unstructured mesh, and the latter initialises with a structured mesh. Despite this, comparing the convergence plot for the square problem produced with SIPG Linear Elastic, Fig. 3, to Fig. 1 of [40] it is clear that the magnitude of the error estimator, η err , the error in the DG norm, ‖ |u -u h ‖ |, are similar for a given number of degrees of freedom, and the convergence patterns also show similarities.

Non-Smooth problem on an L-Shaped domain
Next, we consider example_problem_2.m, a linear elastic problem on (x, y) ∈ Ω = (− 0.5, 0.5) 2 /([0, 0.5] ×[− 0.5, 0]) where x and y are in meters, acting in plane strain with the same material constants as in example 1. The exact solution u is chosen such that the problem is singular at (x, y) = (0, 0), The relevant .m script for this problem is analytical_problem_2_generator.m, and the domain is shown in Fig. 4a. The structure of the code in analytical_problem_2_generator.m is similar to the code in analytical_pro-blem_1_generator.m which is discussed in the previous section. Here, following the same process as in [40], we apply the non-homogenous Dirichlet boundary condition defined by u to all faces and therefore all external faces in etpl_face have the flag -2 associated with them and the symbolic expressions for u,v are set to equal u. We set loop_end=[13 9 3] to perform 13 adaptivity loops for the hp-adaptive simulation, 9 for adaptive-h refinement, and 3 adaptivity loops for the pure-h refinement simulation. The material constants and δ values are the same as the previous example. Again, as this is a problem with a known analytical displacement solution, the analytical body force is calculated using the steps outlined in the previous section and then imposed across the domain.
Calling the poly_plot function found in the plotters folder of SIPG Linear Elastic with the values of etpl and coord at the end of the first Fig. 3. Log-linear convergence of the error estimator, error in the DG and L 2 norms against mesh refinement using different adaptive strategies for the unit square domain problem.
simulation plots a grey-scale colormap of polynomial order distribution across the refined mesh, as shown in Fig. 4b. After the first simulation with hp refinement the mesh has been refined in h extensively around the point x = 0, y = 0. This is because the error estimate for elements near the singularity surpass the threshold set by the δ 2 parameter, and so are refined in h; the polynomial order, controlled using the δ 1 parameter, is highest for elements close to the singularity due to the large errors near this location. However, for the bulk of the mesh (where the solution is smooth) the polynomial order remains low, meaning computational effort is used more efficiently to focus on problem areas. The convergence of the error estimator, η err , the error in the DG norm, ‖ |u -u h ‖ |, and the error in the L 2 norm ‖ u -u h ‖ are shown in Fig. 5, below. Despite beginning with an unstructured mesh, the results of the simulation performed with SIPG Linear Elastic are remarkably similar to those reported in [40], highlighting the ability of the program to analyse non-smooth problems effectively.

Crack in a plate problem
The final example, example_problem_3.m, considers a problem with a stronger singularity, a crack in a plate with the same material properties   To set up this problem we consider the choice of variables in the script nonanalytical_problem_3_generator.m. The initial mesh, shown in Fig. 6a, is unstructured and finer than that used in [40].
Here, we wish to impose the following boundary conditions: nonhomogenous Neumann between nodes 1 and 2, homogenous mixed between nodes 2 and 3, homogenous Dirichlet between nodes 3 and 4, and homogenous Neumann elsewhere. The etpl_face data for this example, stored in crack_EtplFace.txt, therefore provides an insight into how the flagging system in column seven of etpl_face can be used to define a more complex set of boundary conditions. For this example, d_2=[0.3 0.2 0], d_1=[0.07 0.2 0] (NOTE: The first values are different to the value reported in the error estimator paper) and loop_end= [15 15 5].
Unlike the previous examples where the right hand side f could be determined from the analytical displacement solution, for problems without an analytical displacement solution across the domain f must be specified manually. Hence, inputs for the symbolic expressions for the non-homogenous Dirichlet boundary conditions, u and v, nonhomogenous Neumann boundary conditions, traction and body force to be imposed across the domain, fx and fy are required. A final point to note is that a zero input to these variables should be expressed as, for example, traction=[0*X;0*Y]and not[0;0], to conform with the format expected by the symbolic toolbox. Following these guidelines, a zero input is used for u,v,fx and fy, and the traction is specified by setting traction=[0*X; -1].
The convergence of the error estimator for the three adaptive strategies defined by the δ parameters for this problem is shown in Fig. 7. Note that for problems with a non-analytical displacement solution, the DG and L 2 norms will not converge and are therefore not shown in Fig. 7.
The roughly linear convergence of the hp-adaptive strategy in Fig. 7 indicates exponential convergence of the error estimator, whereas the rate of convergence of the error estimator for the h-adaptive and h-uniform strategies decreases as the number of degrees of freedom increases, a result also reported in [40]. This example therefore shows that, by choosing an appropriate adaptive strategy, SIPG Linear Elastic can also be used to analyse non-linear problems with strong singularities.

Conclusions
The paper presented an open source MATLAB code to solve linear elasticity using discontinuous Galerkin finite elements with hp-adaptivity. MATLAB was chosen as the platform for the project because it is one of the best platforms for prototyping numerical methods.
The code has been consciously designed for developers and researchers to be used as a starting point for their projects. The modularity of the code facilitates modifications like implementing a different operator or changing the definition of the error estimator. Customisation can be done without touching the algorithm responsible for the hp-adaptivity.
As shown in literature, hp-adaptivity can improve accuracy for a variety of PDE problems but it is still not commonly used due to the difficulties in implementing it correctly. hp-adaptivity has been included in the package to allow the users to effortlessly adopt it for their projects. As shown in the examples in Section 5, the code is capable to deliver very high accuracy which is better than what many other packages can deliver and much more than what is requested from a prototype. Given this, the present package can be used from prototyping a new numerical method to producing the results to include in publications spanning through the entire research process.