A General Procedure to Formulate 3D Elements for Finite Element Applications

: This paper presents a general procedure to formulate and implement 3D elements of arbitrary order in meshes with multiple element types. This procedure includes obtaining shape functions and integration quadrature and establishing an approach for checking the generated element’s compatibility with adjacent elements’ surfaces. This procedure was implemented in Matlab, using its symbolic and graphics toolbox, and complied as a GUI interface named ShapeGen3D to provide ﬁnite element users with a tool to tailor elements according to their analysis needs. ShapeGen3D also outputs ﬁles with the element formulation needed to enable users to implement the generated elements in other programming languages or through user elements in commercial ﬁnite element software. Currently, ﬁnite element (FE) users are limited to employing element formulation available in the literature, commercial software, or existing element libraries. Thus, the developed procedure implemented in ShapeGen3D offers FEM users the possibility to employ elements beyond those readily available. The procedure was tested by generating the formulation for a brick element, a brick transition element, and higher-order hexahedron and tetrahedron elements that can be used in a spectral ﬁnite element analysis. The formulation obtained for the 20-node element was in perfect agreement with the formulation available in the literature. In addition, the results showed that the interpolation condition was met for all the generated elements, which provides conﬁdence in the implementation of the process. Researchers and educators can use this procedure to efﬁciently develop and illustrate three-dimensional elements


Introduction
The finite element method (FEM) is a popular numerical technique for solving partial differential equations (PDE) in engineering and mathematics.Traditional areas in which this method has been applied include structural analysis, heat transfer, fluid flow, mass transport, electromagnetic potential, and augmented reality [1][2][3][4][5].Under this method, the domain is subdivided into smaller parts, called elements, composed of multiple nodes.The number of nodes within an element determines the order of the interpolation function that represents the approximate solution of the PDE.The interpolation function is defined in terms of nodal basis functions, also known as shape functions.The error associated with the interpolation is minimized using the calculus of variation, resulting in a system of algebraic equations.The variational formulation also requires the integration of finite element variables that are a function of shape functions and other properties, such as shape function derivatives, over the element's domain.These integrals are approximated using a quadrature rule, which requires selecting appropriate integration points and corresponding weights to obtain an exact integration.Finally, each element's equations are combined into a larger system of equations that models the entire problem.
The shape functions and integration quadrature are essential parts of the FEM formulation of an element.Generally, the nodal shape functions are defined in a reference element in local coordinates.If the same shape functions used to interpolate displacements in a reference element are used to map local and global coordinates, they are termed isoparametric elements.The nodal functions and integration quadrature are available for several three-dimensional elements in the literature and are incorporated into commercial FEM software.These elements are generally lower order and limited to four specific elements: hexahedrons from the Lagrangian and Serendipity families, tetrahedrons, prisms, and pyramids [6].For instance, ABAQUS [7] only offers elements of quadratic order to model solids, including a hexahedron from the Serendipity family, also known as the 20-node brick element, tetrahedrons, and prism elements.
Researchers generally use elements only up to the second order of the geometries mentioned above when modeling three-dimensional solids.However, there is a need to employ elements of higher order in many applications.For instance, transition elements are extensively used for mesh refinement [8], where relying only on quadratic elements brings up hanging node problems.Unconventional elements, such as the xNy-element concept [9], are developed to tackle this issue.Transition elements are also employed in contact problems.Buczkowski [10] developed 22 and 21-node elements, and Smith et al. [11] developed 14-node hexahedral isoparametric elements to analyze contact problems by modifying the 8 and 20-node hexahedral elements.In addition to transition elements, there is also a need for higher-order three-dimensional elements as they offer higher accuracy in calculating Lagrangian solid dynamics.For instance, the spectral finite element method [12] uses the interpolation function of high-order Lagrange polynomials to capture high-frequency wave propagation that benefits fields such as structural health monitoring [13][14][15][16] and seismology [17].
Formulating shape functions for elements with arbitrary numbers of nodes in 3D requires defining the nodal distribution, interpolation functions, and integration quadrature.The method for determining the shape functions [3] requires incorporating a polynomial function as the interpolation function in each nodal point.These calculations can be carried out by hand for one-dimensional and lower-order two-dimensional elements.For three-dimensional cases, the procedure is extremely laborious.Another aspect in generating the element formulation is verifying that the shape functions fulfill the following conditions [3]: (a) Interpolation: take a unit value at node i and zero at all other nodes; (b) Local support: vanish over any element boundary that does not include node i; and (c) Interelement compatibility: No discontinuities at the element boundary, i.e., for solid mechanics applications, the elements should deform along the common surface without openings and overlaps.If two elements share a surface, the difference in the value of the shape functions at common nodal points between the two elements should be zero.If the element is a hexahedron or tetrahedron, the need to define the nodal distribution can be avoided by following specific procedures available in the literature.For the hexahedron element of the Lagrangian family, the procedure is explained in [3], whereas in the case of the Serendipity family, Arnold and Arnold [18] provided the theory to obtain shape functions for any given order.In the case of higher-order tetrahedrons, Silvester [19] proposed a simple and elegant procedure based on four-digit node numbering to obtain nodal coordinates that generate a uniform nodal distribution, which is observed to cause Runge oscillations [20].Some alternative nodal distributions to mitigate this effect are Fekete points, Gauss-Legendre points, and Chebyshev points.
There are several methods to obtain the integration quadrature, such as the Gaussian quadrature rule [21] and Legendre-Gauss-Lobatto nodes [22] for hexahedrons.In the case of higher-order tetrahedrons, the standard collapsed coordinate system that maps the quadrature to the tetrahedron from a reference hexahedron [23] is a popular choice.Another approach is based on the integration scheme used for the finite cell method, where a high-order tetrahedron is subdivided into multiple sub-tetrahedrons [20], and numerical integration is carried out on each of them using Gaussian quadrature.
Several existing libraries provide the nodal shapes for higher-order elements that have been obtained a priori.For example, ELEMENTS [24] is a library written in c++ that focuses on supporting novel element spaces, such as Serendipity element space and any order spectral elements.Nekter++ [25] is a FEM solver that supports variable polynomial order in space and heterogeneous polynomial order within two and three-dimensional elements.Scilab [26] is another platform that offers spectral elements in 2D.Some other libraries, such as deal.ii[27], Dune [28], and MFEM [29] are also available.These codes do not generate shape functions nor determine the appropriate integration quadrature for a custom element and do not provide a visual representation.
No methodology exists in the literature for formulating elements of arbitrary order in an automated way.Moreover, if a new element is generated, checks must be performed to guarantee that the interpolation and location support conditions are satisfied.Such a task is difficult to perform for a 3D element as they contain more than three boundaries or surfaces where the local support condition needs to be satisfied for each shape function.In addition, an approach must be established to verify the interelement compatibility of elements in a multielement mesh.Hence, a general procedure is developed in this paper to perform these tasks programmatically.This procedure is capable of generating shape functions and integration quadrature for: (1) arbitrary or custom elements with a given nodal distribution and interpolation function; (2) hexahedrons of the Serendipity and Lagrange families; and (3) tetrahedrons of any given element order.The integration quadrature includes integration points and weights for: (1) volume; and (2) surface integration.The normal vector and parametric coordinates for each element surface are also obtained to enable users to implement the traction vector in FEM applications.The interpolation and local support conditions are checked qualitatively through computer graphics, while interelement compatibility is verified by establishing expressions that ensure nodal shape agreement from adjacent elements at common nodes.
An attractive feature of the proposed procedure is that it can be automated using symbolic and graphical tools, readily available in programming platforms such as Matlab, Python, and Mathematica.In this paper, Matlab was used to implement the procedure, but other users can program the procedure on the platform of their choice.In addition, a Graphical User Interfaces (GUI), named ShapeGen3D (ShapeGen 3D installation file link: https://github.com/sjdyke-reth-institute/ShapeGen3D/blob/main/ShapeGen3D.mlappinstall, accessed on 19 September 2023), was developed to expedite the use of the procedure.GUIs have facilitated the implementation of novel techniques developed by FE code developers, such as in [30,31].Hence, the development of ShapeGen3D can serve as an interactive approach for developing 3D elements and providing finite element users the flexibility to adopt elements of arbitrary order in their simulations.
The remainder of this paper is organized as follows.Section 2 presents the generalized procedure for generating element properties and explains the implementation of the GUI interface, ShapeGen3D.Section 3 obtains the formulation for several 3D elements through the use of ShapeGen3D.The numerical examples are relevant to applications that can be potentially impacted by this work, including the generation of tailored elements designed to reduce the computational cost associated with the analysis (Section 3.3), transition elements that connect fine to coarse mesh regions, such as in contact problems (Section 3.4), and high-order spectral elements to simulate wave propagation behavior (Section 3.5).Section 4 discusses the potential impact of this development in the computational mechanics field.Finally, Section 5 concludes with some final remarks.

Materials and Methods
This section illustrates the general procedure to formulate and implement an element within a multielement FE mesh through computer implementation.The formulation of the element includes three stages, which are (1) the generation of shape functions, (2) integration quadrature, and (3) verification by checking interpolation and local support conditions.If stages 1 and 2 are completed, the element can be incorporated into a multielement mesh if the interelement compatibility is satisfied between the adjacent elements of different types.Hence, stage 3 checks the elements' compatibility.The procedure was outlined in such a way as to make computer implementation possible.This procedure was implemented using a Matlab app with a GUI interface.
Figure 1 breaks down the steps to formulate an element (stages 1 to 3).At the input step, a distinction is made between hexahedron and tetrahedron elements and arbitrary elements.Arbitrary elements, such as transition elements, do not have any specified nodal coordinates and shape functions.Hence, an approach was developed that utilizes (1) nodal coordinates and (2) the monomial basis of the interpolation function to calculate the shape functions, as described in Section 2.1.The procedure for determining the shape functions of hexahedrons and tetrahedrons of elements followed existing literature (see Sections 2.3 and 2.4).Next, in stage 2, the volume and surface integration quadrature are obtained.The integration quadrature in volumetric space includes only integration points and weights, whereas the surface integration quadrature includes integration points, weights, and parametric coordinates for each element's surface.As arbitrary elements do not have any reference quadrature, a method similar to the finite cell method [32,33] was implemented, as described in Section 2.2.Following this approach, the integration quadrature for any element type can be determined for both volume and surface integrations in an automated Next, in stage 2, the volume and surface integration quadrature are obtained.The integration quadrature in volumetric space includes only integration points and weights, whereas the surface integration quadrature includes integration points, weights, and para-metric coordinates for each element's surface.As arbitrary elements do not have any reference quadrature, a method similar to the finite cell method [32,33] was implemented, as described in Section 2.2.Following this approach, the integration quadrature for any element type can be determined for both volume and surface integrations in an automated way, allowing computer implementation.Moreover, the formulation of standard quadratures for hexahedrons and tetrahedrons was implemented as outlined in Sections 2.3 and 2.4, respectively.Stage 3 is a verification step to check the interpolation and local support conditions.This verification was conducted through a graphical representation of the shape function in Matlab (see Section 2.6).The locations where the shape functions yield a 0 value are presented as a void (blank space).
For the final implementation stage, if an element is to be incorporated into a multielement mesh, the interelement compatibility condition between two adjacent elements must be checked.In this paper, a method was developed that can be used to check this condition for each of the integration points of the adjacent surfaces numerically.Such a task is highly labor-intensive to perform by hand for 3D elements, but it can be automized, as presented in Section 2.5. Figure 2 shows the overview of this stage, which requires the input of the surface pairs that will form the interelement boundary.be checked.In this paper, a method was developed that can be used to check this condition for each of the integration points of the adjacent surfaces numerically.Such a task is highly labor-intensive to perform by hand for 3D elements, but it can be automized, as presented in Section 2.5. Figure 2 shows the overview of this stage, which requires the input of the surface pairs that will form the interelement boundary.The layout of this methodology was essential for developing GUI software that automizes the generation of arbitrary 3D element formulation.Finally, Section 2.7 shows the implementation of the GUI interface, which facilitates the execution of the procedure outlined in Figures 1 and 2.

Computation of Shape Functions for a Given Nodal Distribution and Interpolation Function
The generation of shape functions of the development stage is discussed here.Assume that an element defined by a node-set   = [      ] consists of  nodes and an interpolation function   (, , ).If   (, , ) is a linear combination of arbitrary constants   and a function with a monomial basis,   (, , ), then the interpolation function can be written as This function,   , can also be expressed in terms of the shape function terms,   , as where   is a variable associated with node .Plugging Equation (1) into Equation ( 2) The shape functions can be obtained for a given set of known function   by manipulating Equation (3) to solve for parameters   .In matrix form, this operation becomes x The layout of this methodology was essential for developing GUI software that automizes the generation of arbitrary 3D element formulation.Finally, Section 2.7 shows the implementation of the GUI interface, which facilitates the execution of the procedure outlined in Figures 1 and 2.

Computation of Shape Functions for a Given Nodal Distribution and Interpolation Function
The generation of shape functions of the development stage is discussed here.Assume that an element defined by a node-set P i = x i y i z i consists of n nodes and an interpolation function f n (x, y, z).If f n (x, y, z) is a linear combination of arbitrary constants c i and a function with a monomial basis, F i (x, y, z), then the interpolation function can be written as This function, f n , can also be expressed in terms of the shape function terms, N i , as where u i is a variable associated with node i. Plugging Equation (1) into Equation ( 2) The shape functions can be obtained for a given set of known function F i by manipulating Equation (3) to solve for parameters c i .In matrix form, this operation becomes x Plugging c i in Equation ( 3) and isolating u i , yields the shape functions [4].If the interpolation function f n is unable to form a non-invertible matrix F with node set P, then the node-set and interpolation function are not valid.If invertible, shape functions can be formulated that automatically satisfy the interpolation and partition of unity condition.This method cannot be used to generate shape functions for rational polynomials, such as in the case of Wachspress basis functions [34] and Radial basis functions [35].This approach for determining shape functions was implemented using the symbolic toolbox in Matlab.The steps implemented in the Matlab routines are summarized in Figure 3. Plugging   in Equation ( 3) and isolating   , yields the shape functions [4].If the interpolation function   is unable to form a non-invertible matrix  with node set , then the node-set and interpolation function are not valid.If invertible, shape functions can be formulated that automatically satisfy the interpolation and partition of unity condition.This method cannot be used to generate shape functions for rational polynomials, such as in the case of Wachspress basis functions [34] and Radial basis functions [35].This approach for determining shape functions was implemented using the symbolic toolbox in Matlab.The steps implemented in the Matlab routines are summarized in Figure 3.

Determination of Integration Quadrature for an Arbitrary Element
The integration quadrature for an arbitrary element is determined following an approach similar to the finite cell method [32,33], which involves dividing the element into a finite number of tetrahedron cells and mapping the integration points from a reference tetrahedron to each cell.First, the Delaunay triangulation [36] is carried out for the given set of node , which subdivides the arbitrary element into multiple tetrahedrons in the  space.Each of the tetrahedrons is termed as   , where the subscript  indicates the element number with nodal coordinates (  ,   ,   ), and  can take values from 1 to 4. Next, a reference element, consisting of a standard orthogonal tetrahedron in  space with vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1) is formulated.An integration quadrature, i.e., integration points and weights, up to order 4 is used from [37] for such refer-

Determination of Integration Quadrature for an Arbitrary Element
The integration quadrature for an arbitrary element is determined following an approach similar to the finite cell method [32,33], which involves dividing the element into a finite number of tetrahedron cells and mapping the integration points from a reference tetrahedron to each cell.First, the Delaunay triangulation [36] is carried out for the given set of node P, which subdivides the arbitrary element into multiple tetrahedrons in the xyz space.Each of the tetrahedrons is termed as S e , where the subscript e indicates the element number with nodal coordinates (x n , y n , z n ), and n can take values from 1 to 4. Next, a reference element, consisting of a standard orthogonal tetrahedron in ξηΓ space with vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1) is formulated.An integration quadrature, i.e., integration points and weights, up to order 4 is used from [37] for such reference tetrahedron element.The shape functions of the reference 4-noded tetrahedron element are Finally, the integration quadrature of the global cells or tetrahedrons is calculated.Each integration point coordinate (ξ i , η i , Γ i ) from the reference tetrahedron in ξηΓ space is mapped to each S e in the xyz space following Equation (7) [3] The weight of the integration points for the global element w x i , where superscript x indicates xyz space, is mapped by multiplying the weights of integration points of the reference element w ξ i with the determinant of Jacobian, as Following the isoparametric mapping, the Jacobian matrix can be determined as For linear tetrahedrons, i.e., the reference element, The integration quadrature on each surface is also determined to allow the computation of equivalent nodal flux.First, the surfaces of an element must be identified.This step can be accomplished by determining the convex hull [36], which is the smallest convex set that bounds all the nodes of an element.Hence, the convex hull of the Delaunay triangulation is determined to obtain the integration points and weights at the surfaces.Each surface s of the set is a triangle with three vertices (a n , b n , c n ), where n is from 1 to 3. The integration point coordinates for a surface s are determined as, The normal vector for the surface, pointing outside the element interior, can be determined by performing the following cross-product operation as, where v indicates the second norm of vector v.The parametric coordinates of a surface s, s n x , s n y can be determined as The area of the surface is the weight of the integration, which can be determined as

Formulation of Nodal Coordinates and Integration Quadrature for Hexahedron Elements
A standard orthogonal hexahedron in xyz space is considered to formulate hexahedron elements with eight vertices as ( . If c n is an arbitrary constant, the interpolation function for Lagrange family elements can be written as, where m x , m y , and m z are the polynomial order along x, y, and z axis.For the Serendipity family, the interpolation function for such element can be obtained as, Many studies and codes [37,38] follow a uniform nodal distribution that exhibits Runge oscillations.To avoid this numerical artifact, this development follows Legendre-Gauss-Lobatto nodes [22] where the node coordinates for each of the three axes, x P, y P, z P, is the solution of Equations ( 17)-( 19) between +1 and −1, respectively.
where, L m x indicates the Legendre polynomial of the degree m x .Each of these vectors of Equations ( 17)-( 19) comprise a set of (m x + 1), m y + 1 , and (m z + 1), respectively.To obtain the nodal coordinates in three dimensions, a mesh grid of size m y + 1 × (m x + 1) × (m z + 1) is generated using the obtained one-dimensional grids.
For the Lagrange family, all the mesh grid points were considered, whereas only the points on the element edges were considered for the Serendipity family elements.
For integration quadrature, Legendre-Gauss-Lobatto quadrature and Gauss-Legendre quadrature were computed for a given element order along the x, y, and z-axis.Similar to node generation, one-dimensional reference grids were generated at first.For the Legendre-Gauss-Lobatto quadrature [22], the weights x w, for each component of a one-dimensional reference grid in the x axis can be determined as For the Gauss-Legendre quadrature [21], the coordinate of integration points in one dimension, Ix P can be determined as the solutions of Equation ( 21) between +1 and −1 as The corresponding weight components of x w can be determined as Considering that y w and z w are reference grids for the weights along y and z axes, the weights w for the integration points in 3D can be determined using the tensor product [39] as presented in Equation ( 23).The components of the product w are then assigned to the corresponding points.
where " * " indicates elementwise multiplication.For the surfaces, the interpolation order for each surface is determined before computing the integration quadrature.If a surface s contains node i, the number of differentiations along each axis of the shape function that produces a constant can be set as the interpolation order.This operation can be performed for all the nodes contained by s, where the maximum order is chosen for the surface.For s m x , this condition can be defined using a constant C as presented in Equation (24), which is carried out using the symbolic toolbox of Matlab.
The coordinates and weights of the integration points and parametric coordinates for each of the six surfaces for hexahedron elements, x = 1, x = −1, y = 1, y = −1, z = 1, and z = −1 are defined according to Table 1.
Table 1.Integration points weights and parametric coordinates of the element's surfaces.

Formulation of Nodal Coordinates and Integration Quadrature for Tetrahedron Elements
This section presents the methodology implemented in the general procedure to formulate tetrahedron elements.The Lobatto triangle grid for orthogonal polynomials proposed by Luo et al. [40] is implemented to generate element properties for the standard orthogonal tetrahedron in xyz space, which consists of four vertices as (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1).This nodal distribution does not exhibit Runge oscillations.According to this method, if an m th degree interpolating polynomial is used, the interpolation function can be expressed as, c n x i y j z k where, (i where c n is an arbitrary constant.A one-dimensional reference grid comprising a set of m − 1 degrees is generated as, where x P is the solution of Equation ( 26) between the range +1 and −1 as, The L m is the Legendre polynomial, which can be written as, Then, using the grid, the i node coordinate, x i y i z i can be determined as,
For elements with order ≤ 4, the integration quadrature is used from [37], whereas the integration quadrature for elements with order > 4 is determined using the finite cell method.According to this method, the element is subdivided into multiple sub-tetrahedrons, and the integration points and weights for each of these sub-tetrahedrons are then determined based on the Gaussian quadrature.The decomposition of the tetrahedron into 8 and 64 sub-tetrahedrons is shown in Figure 4.The method adopted for the arbitrary element is followed to determine the integration quadrature at the surface with the parametric coordinates presented in Table 2.

Interelement Compatibility Check between Two Elements for a Chosen Surface Pair
The nodal compatibility condition requires that the nodal values of a variable evaluated at the common nodes of adjacent elements must be equal.This condition will be satisfied if all the nodes of the interelement surface have the same coordinates and satisfy the interpolation condition.The interelement compatibility condition will be satisfied if the value of a variable is equal for all the points on the interelement surface.If the interelement condition is satisfied, the nodal compatibility condition is satisfied by default.This paper develops a new approach to check the interelement compatibility condition between two elements for a given interelement surface using the symbolic toolbox of Matlab.The approach consists of forming a two-element mesh and ensuring the nodal shapes match at the integration points of the interelement surface.
Assume two elements, e 1 and e 2 , are formulated in the coordinates x = x y z , and that the compatibility condition needs to be checked between surface p for element 1 and surface q for element 2. To form an interelement surface Γ by p and q, element 2 must go through coordinate transformation (x → X ) such that the normal surface vectors cancel each other, and the surface centroids have the same coordinates as presented in the equations below A coordinate translation followed by a rotation of the element nodal coordinates is performed to satisfy Equations ( 34) and (35).Figure 5 presents the elements in the local coordinate x, and Figure 6 shows element 2 in the X coordinates after the transformation is conducted.For this purpose, the parametric coordinates of the surfaces of both elements are utilized.A coordinate system, 1 e, with − p n z as its third component, indicates that this normal vector points inward to the element as elements for a given interelement surface using the symbolic toolbox of Matlab.The approach consists of forming a two-element mesh and ensuring the nodal shapes match at the integration points of the interelement surface.Assume two elements,  1 and  2 , are formulated in the coordinates  = (    ), and that the compatibility condition needs to be checked between surface  for element 1 and surface  for element 2. To form an interelement surface  by  and , element 2 must go through coordinate transformation ( → ) such that the normal surface vectors cancel each other, and the surface centroids have the same coordinates as presented in the equations below A coordinate translation followed by a rotation of the element nodal coordinates is performed to satisfy Equations ( 34) and (35).Figure 5 presents the elements in the local coordinate , and Figure 6 shows element 2 in the  coordinates after the transformation is conducted.For this purpose, the parametric coordinates of the surfaces of both elements are utilized.A coordinate system,  1 , with −    as its third component, indicates that this normal vector points inward to the element as For element 2, a coordinate system can be constructed as For element 2, a coordinate system can be constructed as 2 e = q n x q n y q n z (37) According to Equation ( 34), element 2 needs to be rotated such that the q n z aligns with − p n z .The resulting transformation of the x coordinates is Equation (38) facilitates the formation of an interelement boundary between p and q.The rotation matrix, R θ , about the q n z axis, and the coordinate translation is performed to satisfy Equation (35) as where 2 X Rq is the centroid of the surface q after coordinate rotation and R θ Next, the shape functions are expressed in terms of the global coordinate X following a similar approach as The variable substitution from x to X is performed using the Symbolic toolbox in Matlab.For an element to undergo such coordinate transformation, x can be mapped (x→X) using isoparametric formulation as, After the coordinate transformation, the compatibility condition can be checked following the procedure presented below.A variable u at x of element 1 can be determined as where 1 N j is shape function j of element 1 and x j / ∈ p indicates that node j does not fall in the surface p.Similarly, a variable u at x of an element 2 can be determined as The interelement compatibility condition holds if the variable 1 u(x) = 1 u(X), on the interelement surface Γ.Following Equations ( 42)-(44) can be written as If the local support condition is satisfied, the first and third components of Equation (45) will be zero as Using Equations ( 46) and (47) and performing the assembly by setting 2 u m = 1 u k = Γ u k , the compatibility conditions hold if For every x ∈ Γ, and Γ u k Equation (46) will be satisfied if Equation ( 49) can be integrated numerically for each integration point as Similarly, Equations ( 46) and ( 47) can be checked numerically as Hence, the compatibility condition is satisfied if Equations ( 50)-( 52) are satisfied for a given element surface pair between two elements.
Next, the shape functions are expressed in terms of the global coordinate  following a similar approach as The variable substitution from  to  is performed using the Symbolic toolbox in

Generation of Graphical Representation of Element Properties
Element properties, including (a) shape functions, N i (b) differentiation of N i along the x, y, and z-axis, and (c) integration points and weights, are plotted using the Matlab graphics toolbox.These plots facilitate the understanding of the element interpolation function and help to check the interpolation and local support conditions.
The Matlab graphics toolbox is utilized to create a 3D graphical representation of elements' properties, such as shape function value.For the hexahedron element, the element is subdivided into a 3D grid to perform plotting.For tetrahedrons and arbitrary elements, each surface is plotted.The fidelity of the plot can be increased by discretizing the reference element further by tetrahedron partitioning, and subsequently, each point's coordinates can be mapped to S e following Equation (7). Figure 7 presents an example to illustrate the process with a pyramid element.In step 1, the element is subdivided into two tetrahedrons, S 1 and S 2 .In steps 2 and 3, a reference tetrahedron is generated and subdivided into multiple grids to increase the graphics fidelity, and finally, in step 4, each grid point is mapped to the global tetrahedron.
ment is subdivided into a 3D grid to perform plotting.For tetrahedrons and arbitrary el-ements, each surface is plotted.The fidelity of the plot can be increased by discretizing the reference element further by tetrahedron partitioning, and subsequently, each point's coordinates can be mapped to   following Equation (7). Figure 7 presents an example to illustrate the process with a pyramid element.In step 1, the element is subdivided into two tetrahedrons,  1 and  2 .In steps 2 and 3, a reference tetrahedron is generated and subdivided into multiple grids to increase the graphics fidelity, and finally, in step 4, each grid point is mapped to the global tetrahedron.

Implementation of GUI Interface
The procedure is compiled as a GUI software: ShapeGen3D, and organized into four segments, as presented in Figure 8, to enable users to navigate, formulate, and verify the accuracy of the formulation.Segment 1 is the input for the software that lets users choose the element family of interest and the desired interpolation order.Segment 2 includes the execution command to run the software and obtain the element formulation.Segment 3 allows users to store the element properties.Segment 4 represents the results and post-

Implementation of GUI Interface
The procedure is compiled as a GUI software: ShapeGen3D, and organized into four segments, as presented in Figure 8, to enable users to navigate, formulate, and verify the accuracy of the formulation.Segment 1 is the input for the software that lets users choose the element family of interest and the desired interpolation order.Segment 2 includes the execution command to run the software and obtain the element formulation.Segment 3 allows users to store the element properties.Segment 4 represents the results and postprocessing blocks that will enable users to observe and manipulate plots that can be used to verify the accuracy of the generated formulation visually.Detailed descriptions for each of these segments are provided below.
A short description and a screenshot of the sub-options in Segment 1 are presented in Table 3 and Figure 9, respectively.The element order and integration quadrature type must be defined for developing hexahedron and tetrahedron elements.On the other hand, in the case of arbitrary or custom elements, the nodal distribution must be inputted by users, and the software will determine the integration quadrature.
Table 3. Segment 1 sub-options and description.

Seg 1: Input Options Sub-Options Description
Hexahedron or Brick element Polynomial degree along x m x of Equation ( 17) Polynomial degree along y m y of Equation ( 18) Polynomial degree along z m z of Equation (19) Lagrangian Family Checking this option will formulate a Lagrangian element.Otherwise, a Serendipity element will be formulated

Seg 1: Input Options Sub-Options Description
Tetrahedron element Element order m in Equation ( 27)

Number of Sub-Tetrahedral
The number of sub-tetrahedrons used for integration quadrature following the finite cell method.It is active when an element order greater than 4 is used.
Integration order for reference tetrahedron The integration order for each sub-tetrahedron.It is active when an element order greater than 4 is used.

Arbitrary or Custom element Input node coordinates
A table where the user can input the nodal distribution of the desired element and cartesian coordinates (x, y, z) for each node.

Input interpolation function
A table where the user can input F i (x, y, z) from Equation ( 3), where i = [1 d] and d is the number of nodes.

Integration quadrature
If the software detects an arbitrary element from the isoparametric hexahedron family, it will determine the selected integration quadrature from one of the selected options: (1) Gauss-Legendre and (2) Gauss-Lobatto.

Integration order for reference tetrahedron
If the element is not from an isoparametric hexahedron family, the integration quadrature is determined following the procedure presented in Section 2.2 with the integration order selected in this option processing blocks that will enable users to observe and manipulate plots that can be used to verify the accuracy of the generated formulation visually.Detailed descriptions for each of these segments are provided below.
A short description and a screenshot of the sub-options in Segment 1 are presented in Table 3 and Figure 9, respectively.The element order and integration quadrature type must be defined for developing hexahedron and tetrahedron elements.On the other hand, in the case of arbitrary or custom elements, the nodal distribution must be inputted by users, and the software will determine the integration quadrature.

Integration order for reference tetrahedron
If the element is not from an isoparametric hexahedron family, the integration quadrature is determined following the procedure presented in Section 2.2 with the integration order selected in this option Check interelement compatibility Used to illustrate two consecutively generated elements Segment 2 represents the "Execution" block.It includes the "Generate shape function" button, which executes the code with the input given in Segment 1 and generates shape functions, integration quadrature, and graphical representation of the element properties.The code can automatically detect isoparametric hexahedrons and tetrahedron Segment 2 represents the "Execution" block.It includes the "Generate shape function" button, which executes the code with the input given in Segment 1 and generates shape functions, integration quadrature, and graphical representation of the element properties.The code can automatically detect isoparametric hexahedrons and tetrahedron elements.There are several checkboxes, such as the "Do not guess" command.If this box is checked, the software avoids detecting the element type and assumes it is an arbitrary element.This segment also includes a display window that shows messages and the progress of the execution of the code.If the "Check interelement compatibility" box is selected in Segment 1, this option will be activated in Segment 2. The software will compute the nodes with the same coordinates between two consecutive elements and subtract the corresponding shape functions.If the difference is 0 at an interelement boundary, that boundary satisfies the compatibility condition between two elements, i.e., there are no discontinuities at the boundary.
Segment 3 is a block that outputs results in the Matlab workspace to allow users to implement the elements in other programming languages or through user elements in commercial finite element software, such as ABAQUS and ANSYS.Users can also output shape functions and derivatives of the shape functions evaluated at the integration points by ticking the "Shape matrix" and "Derivatives of Shape matrix" options.A description of the workspace files is presented in Table 4.These files are readily available for FEM calculations.a, b, d are node numbers and integers.n j (n x j , n y j , n z j ) is the normal vector pointing outward of the domain on the surface containing node j.

Integration_Points_Coordinates_Weights
Stores integration points and weights inside the element g j (g x j , g y j , g z j ) where g x j : x coordinate of integration point j.w j : weight of integration point j.

Integration_Points_Coordinates_Weights_Vectors on_Surface
Stores integration points and weights on surfaces and surface normal vector pointing outward of the domain for all surfaces g j (g x j , g y j , g z j ) where g x j : x coordinate of integration point j.w j : weight of integration point.j. n j (n x j , n y j , n z j ) is the normal unit vector pointing outward of the domain on the surface at integration point j.

Parametric_Coordinates
Stores parametric coordinates for each surface of the element x j and n y j are the parametric x and y coordinates, respectively.S represents the surface number.

Shape_Matrix
Stores the value of each shape function evaluated at each integration point N i | g j : shape function i evaluated at integration point j.
dx_Shape_Matrix Stores the value of the first derivative of each shape function evaluated at each integration point : shape function derivative i evaluated at integration point j.
Finally, segment 4 represents a results and post-processing block that shows the results obtained by the software and allows users to manipulate graphics.There are many components, such as the "Shape function, N" spinner block, which is a control block where the user can select to plot the number of shape function properties, such as their value and derivatives.This segment also allows users to plot other element properties, such as node number, integration points, and weights, as shown in the bottom left part of the software.It also includes a display window that shows the shape function property in equation format and other user-controlled options to manipulate the graphics options, such as the view angle, axis, image resolution (Peels), and transparency (Fill).Increasing the "Peels" will increase the overall image resolution.If "Fill" is chosen as 0, the software will not show any 0 value and will represent it as a blank space.

Numerical Examples
This section shows the properties of three different 3D elements that were obtained following the procedure developed in this paper.The results were generated through the GUI software: ShapeGen3D.First, the element properties, i.e., nodal coordinates, shape functions, and integration quadrature, of the popular 20-node brick element were obtained and compared against the literature to verify the accuracy of the formulation.Second, a 21-node hexahedron transition element similar to a 20-node quadratic hexahedron element but with one node on the xy face at z = 1 was formulated.The interpolation, local support condition, and compatibility with the 20-node brick element were checked.Third, a custom element was formulated to replicate the Germain-Lagrange plate [41] bending, showing how the methodology implemented in the software can help solve challenging problems efficiently.The fourth example shows the development of a 43-node transition element from a fourth to second-order Lagrangian element.Finally, two high-order spectral elements were formulated to provide a complete overview.The steps followed in ShapeGen3D to generate the elements' shape functions and other finite element properties for each of the examples are presented in Table A9.

20-Node Brick Element
A 20-node brick element follows serendipity element formulation with m x = 2, m y = 2, and m z = 2.The software determines the nodal coordinates, shape functions (Stage 1), and integration quadrature (Stage 2).The obtained results were saved as the files explained in Table 4 and are shown in the Appendix A, see Tables A1-A5, which includes: (1) nodal coordinates; (2) corresponding shape functions; (3) integration points and weights; (4) integration points on the surface and weights; and (5) nodes on the surface and surface normal vector pointing outward of the domain for all surfaces, respectively.The nodal coordinates, shape functions, and integration quadrature match the formulation found in the literature [37], thus verifying the accuracy of the software results.Figure 10 represents the input to the ShapeGen3D and the obtained results with nodes as black dots and integration points for volume integration as red dots.

21-Node Brick Element
Next, the 21-node brick element was formulated to illustrate the development of a custom element.For this example, the node coordinates and the interpolation function generated for the 20-node element were imported to "Arbitrary or Custom element" by pressing the "Import previous element" button.Then, one node was added manually with coordinate (0,0,−1), and several trials for adding a polynomial term associated with coefficient c 21 were made until an invertible F matrix (Equation ( 5)) was obtained by the software.The additional polynomial basis that was added to the interpolation function was x 2 y 2 (z + 1), which yields 2x 2 y 2 at the z = 1 surface that contains the 21st mid-surface node and vanishes at the z = −1 surface without a mid-surface node.The interpolation function is presented in Equation (53) as, The input to ShapeGen3D is presented in Figure 11, which shows the input coordinate of the 21st node and the polynomial term with coefficient c 21 .The obtained nodal coordinates, corresponding shape functions, and nodes on the surface and surface normal vector pointing outward of the domain for all surfaces are presented in Tables A6-A8, respectively.Other element variables, such as integration points and weights, match those of the 20-node brick element shown in Table A8.explained in Table 4 and are shown in the Appendix, see Tables A1-A5, which includes: (1) nodal coordinates; (2) corresponding shape functions; (3) integration points and weights; (4) integration points on the surface and weights; and (5) nodes on the surface and surface normal vector pointing outward of the domain for all surfaces, respectively.The nodal coordinates, shape functions, and integration quadrature match the formulation found in the literature [37], thus verifying the accuracy of the software results.Figure 10 represents the input to the ShapeGen3D and the obtained results with nodes as black dots and integration points for volume integration as red dots.

21-Node Brick Element
Next, the 21-node brick element was formulated to illustrate the development of a custom element.For this example, the node coordinates and the interpolation function generated for the 20-node element were imported to "Arbitrary or Custom element" by pressing the "Import previous element" button.Then, one node was added manually with coordinate (0,0,−1), and several trials for adding a polynomial term associated with coefficient  21 were made until an invertible  matrix (Equation ( 5)) was obtained by the software.The additional polynomial basis that was added to the interpolation function was  2  2 ( + 1), which yields 2 2  2 at the  = 1 surface that contains the 21st midsurface node and vanishes at the  = −1 surface without a mid-surface node.The interpolation function is presented in Equation ( 53 The input to ShapeGen3D is presented in Figure 11, which shows the input coordinate of the 21st node and the polynomial term with coefficient  21 .The obtained nodal coordinates, corresponding shape functions, and nodes on the surface and surface normal For stage 3, each shape function can be analyzed to check the accuracy of the calculations performed by the software.The element must satisfy two conditions: (1) the interpolation condition; and (2) local support.As stated in Section 2, the interpolation condition was satisfied when an invertible matrix F of Equation ( 5) was obtained.The graphical illustrations of each shape function can be used to check if the local support condition is satisfied.For example, Figure 12 plots shape functions at nodes 5 and 21, N 5 and N 21 , where the shape functions are shown to vanish at the face that does not contain the nodes, and thus, satisfying the local support condition.lations performed by the software.The element must satisfy two conditions: (1) the interpolation condition; and (2) local support.As stated in Section 2, the interpolation condition was satisfied when an invertible matrix F of Equation ( 5) was obtained.The graphical illustrations of each shape function can be used to check if the local support condition is satisfied.For example, Figure 12 plots shape functions at nodes 5 and 21,  5 and  21 , where the shape functions are shown to vanish at the face that does not contain the nodes, and thus, satisfying the local support condition.Next, on stage 4, the compatibility with the 20-node brick element was checked.For this purpose, the "Check compatibility" button was selected.Figure 13 shows the results of the intercompatibility check of the first surface pair between the two elements, which was performed according to the procedure described in Section 2.5.The green lamp next to the compatibility text in Figure 10 indicates that the interelement compatibility condition was satisfied.A red lamp would indicate that the intercompatibility condition was not met.This condition can also be verified qualitatively using the plots of shape functions corresponding to the connecting nodes, which show the same shape function profile at the element boundary.By checking all of the surface pairs, it is concluded that any surface can be the interelement surface between the 20 and 21-node elements, except the  face at  = 1.Next, on stage 4, the compatibility with the 20-node brick element was checked.For this purpose, the "Check compatibility" button was selected.Figure 13 shows the results of the intercompatibility check of the first surface pair between the two elements, which was performed according to the procedure described in Section 2.5.The green lamp next to the compatibility text in Figure 10 indicates that the interelement compatibility condition was satisfied.A red lamp would indicate that the intercompatibility condition was not met.This condition can also be verified qualitatively using the plots of shape functions corresponding to the connecting nodes, which show the same shape function profile at the element boundary.By checking all of the surface pairs, it is concluded that any surface can be the interelement surface between the 20 and 21-node elements, except the xy face at z = 1.
Next, on stage 4, the compatibility with the 20-node brick element was checked.For this purpose, the "Check compatibility" button was selected.Figure 13 shows the results of the intercompatibility check of the first surface pair between the two elements, which was performed according to the procedure described in Section 2.5.The green lamp next to the compatibility text in Figure 10 indicates that the interelement compatibility condition was satisfied.A red lamp would indicate that the intercompatibility condition was not met.This condition can also be verified qualitatively using the plots of shape functions corresponding to the connecting nodes, which show the same shape function profile at the element boundary.By checking all of the surface pairs, it is concluded that any surface can be the interelement surface between the 20 and 21-node elements, except the  face at  = 1.

43-Node Hexahedron Element
A 43-node hexahedron element was formulated following a similar procedure as the 21-node, such that the interpolation function provided a solution to the plate bending equation.According to the Kirchhoff-Love plate theory [41] for isotropic plates, the displacement along the z-axis, , follows the differential equation as, where  is the constant distributed load, and  is the flexure rigidity.If the time derivative is canceled, it forms the Germain-Lagrange plate equation.A polynomial solution of this differential equation can be written as, A 3D element was formulated keeping the lowest order as quadratic along the z-axis such that the interpolation function contained the polynomial basis presented in Equation

43-Node Hexahedron Element
A 43-node hexahedron element was formulated following a similar procedure as the 21-node, such that the interpolation function provided a solution to the plate bending equation.According to the Kirchhoff-Love plate theory [41] for isotropic plates, the displacement along the z-axis, w, follows the differential equation as, 2ρh D where q is the constant distributed load, and D is the flexure rigidity.If the time derivative is canceled, it forms the Germain-Lagrange plate equation.A polynomial solution of this differential equation can be written as, A 3D element was formulated keeping the lowest order as quadratic along the zaxis such that the interpolation function contained the polynomial basis presented in Equation (55).First, a 36-node hexahedron element of the Serendipity family was formulated with the options m x = 4, m y = 4 and m z = 2.This interpolation function of this Serendipity element does not contain the x 2 y 2 term, hence, it cannot replicate the solution of Equation (54) exactly.Thus, the element was imported as an arbitrary element, and seven additional nodes with coordinates (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), and (0, 0, 0) were added.Finally, seven polynomial basis as x 2 y 2 z 2 , x 2 y 2 z, x 2 yz 2 , xy 2 z 2 , x 2 y 2 , z 2 y 2 and x 2 z 2 were added to formulate the 43-node element that has the polynomial interpolation function as presented in Equation (56).The screenshot of this element formulation is presented in Figure 14.(55).First, a 36-node hexahedron element of the Serendipity family was formulated with the options   = 4,   = 4 and   = 2.This interpolation function of this Serendipity element does not contain the  2  2 term, hence, it cannot replicate the solution of Equation (54) exactly.Thus, the element was imported as an arbitrary element, and seven additional nodes with coordinates (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), and (0, 0, 0) were added.Finally, seven polynomial basis as  2  2  2 ,  2  2  ,  2  2 ,  2  2 ,  2  2 ,  2  2 and  2  2 were added to formulate the 43-node element that has the polynomial interpolation function as presented in Equation (56).The screenshot of this element formulation is presented in Figure 14.

43-Node Transition Element from Fourth to Second-Order Lagrangian Element
The fourth-order Lagrangian element is a 125-node element (  =   =   = 4) with five nodes along each axis.Such high-order elements can replicate wave propagation and are widely used in the structural health monitoring field [13][14][15][16].The second-order Lagrangian element is a 27-node element (  =   =   = 2) with three nodes along each axis.In this example, a transition element was developed to bridge these two elements.The transition element is considered to be of second order along the -axis.
The  = 1 surface of the transition element should match the  = ∓1 surface of the 125-node element, and  = −1 surface should match the 27-node element's surface.Hence, there will be 25 nodes on the  = 1 surface and nine nodes on the  = −1 surface.As the order of the transition element is two, instead of a 125-node element, a 75node element with   =   = 4 and   = 2 was considered as the reference

43-Node Transition Element from Fourth to Second-Order Lagrangian Element
The fourth-order Lagrangian element is a 125-node element m x = m y = m z = 4 with five nodes along each axis.Such high-order elements can replicate wave propagation and are widely used in the structural health monitoring field [13][14][15][16].The second-order Lagrangian element is a 27-node element ( m x = m y = m z = 2 with three nodes along each axis.In this example, a transition element was developed to bridge these two elements.The transition element is considered to be of second order along the z-axis.
The z = 1 surface of the transition element should match the z = ∓1 surface of the 125-node element, and z = −1 surface should match the 27-node element's surface.Hence, there will be 25 nodes on the z = 1 surface and nine nodes on the z = −1 surface.As the order of the transition element is two, instead of a 125-node element, a 75-node element with m x = m y = 4 and m z = 2 was considered as the reference element.Upon formulation of the 75-node element, the nodes at the z = −1 surface are removed to replicate the surface of a 27-node element.
The minimum order of the element is 2. Hence, all the monomial basis functions of order 2 are kept.As the element needs to transition from order 4 to 2 along the z-axis, the monomial basis that has order m z = 1 was kept and manipulated as (z + 1), such that those terms vanish at the z = −1 surface, as presented in Table 5.The remaining 43 terms form the interpolation function of the 43-node transition element.
Factors with m z = 1 x 3 z x 3 (z + 1) x 4 z x 4 (z + 1) x 3 yz x 3 y(z + 1) x 4 yz x 4 y(z + 1) x 2 y 3 z x 2 y 3 (z + 1) x 4 y 2 z x 4 y 2 (z + 1) x 2 y 4 z x 2 y 4 (z + 1) x 4 y 3 z x 4 y 3 (z + 1) x 3 y 4 z x 3 y 4 (z + 1) x 3 y 3 z x 3 y 3 (z + 1) x 4 y 4 z x 4 y 4 (z + 1) The generated element satisfies the local support conditions presented in Figure 15 for two shape functions, nodes 6 and 35.The inter-element compatibility conditions were also checked, and compatibility was observed between the z = 1 surface of the transition element and the fourth-order element and between the z = −1 surface of the transition and the second-order element.Figure 16 presents a graphical overview of the compatibility check.The red hollow diamond represents the nodes of the second-order (Figure 16a) and fourth-order (Figure 16b) elements, whereas the black dots represent the transition element.It is evident that the node coordinates match, and the shape functions agree at the common nodes.

Spectral Elements
Finally, the shape functions were generated for a hexahedron element of the Lagrangian family with 108 nodes to provide a complete overview of the methodology.In ShapeGen3D, m x = 5, m y = 5, and m z = 2 were provided as input parameters.The shape function for node 70 of the element is presented in Figure 17a.The red spot of the plot indicates the location of this node.It can be observed that the shape function satisfies the interpolation and local support conditions as it vanishes to 0 over surfaces that do not contain the node.The same behavior is observed for other nodes.This spectral element was used for simulating guided waves by Soman et al. [15].also checked, and compatibility was observed between the  = 1 surface of the transition element and the fourth-order element and between the  = −1 surface of the transition and the second-order element.Figure 16 presents a graphical overview of the compatibility check.The red hollow diamond represents the nodes of the second-order (Figure 16a) and fourth-order (Figure 16b) elements, whereas the black dots represent the transition element.It is evident that the node coordinates match, and the shape functions agree at the common nodes.Factors with   = 1  3   3 ( + 1)  3   3 ( + 1)  4   4 ( + 1)  3   3 ( + 1)  3   3 ( + 1)  4   4 ( + 1)  4   4 ( + 1)

Spectral elements
Finally, the shape functions were generated for a hexahedron element of the Lagrangian family with 108 nodes to provide a complete overview of the methodology.In Shape-Gen3D,   = 5,   = 5, and   = 2 were provided as input parameters.The shape function for node 70 of the element is presented in Figure 17a.The red spot of the plot indicates the location of this node.It can be observed that the shape function satisfies the interpolation and local support conditions as it vanishes to 0 over surfaces that do not contain the node.The same behavior is observed for other nodes.This spectral element was used for simulating guided waves by Soman et al. [15].
Additionally, a 5th-order tetrahedron element (56 nodes) with input  = 5 (Equa- Additionally, a 5th-order tetrahedron element (56 nodes) with input m = 5 (Equation ( 13)) was formulated.The shape function for the 50th node of the element is presented in Figure 17b.All shape functions for this element were also observed to satisfy the interpolation and local support conditions.Thus, analysts using the spectral finite element method can use ShapeGen3D to determine the formulation for arbitrary noded elements.

Discussion
This development reduces implementation effort in formulating and verifying arbitrary 3D elements for FEM applications.It provides FEM users the possibility to employ elements beyond those available in commercial FEM software, FEM libraries, and openaccess codes.The proposed procedure facilitates the formulation of higher-order elements, which can be shown to solve challenging problems with fewer degrees of freedom in some applications.The software was developed as part of a NASA-sponsored space technology research institute, Resilient Extra-Terrestrial Habitat Institute (RETHi), which aims to propel space exploration forward by developing new knowledge, technologies, and techniques.ShapeGen3D was used to develop the 21-node brick elements presented in this paper to capture meteoroid impact loading on a space habitat model.The model was embedded within a Modular coupled virtual testbed (MCVT) [42][43] that integrates several sub-systems of a space habitat to evaluate their interdependence and the effectiveness of safety controls when a habitat is exposed to external disturbances.ShapeGen3D allowed developers to easily adjust the elements of the structural model until the functional requirements were met.A 43-node brick element was also developed with Shape-Gen3D to mimic wave propagation due to meteoroid impact and assist with generating a structural health monitoring algorithm for detecting damage.
Upon publication, this development is expected to benefit the computational mechanics community beyond the RETHi group.The higher-order and transition element formulation that the software can generate will allow researchers to employ elements that can increase the accuracy and reduce the computational time when dealing with simulations that involve contact, mesh transition regions, and Lagrangian solid dynamics.Furthermore, ShapeGen3D can become an educational tool to illustrate three-dimensional FEM concepts.

Conclusions
This work presents the development of a generalized procedure for formulating nodal coordinates, shape functions, and integration quadrature of higher-order hexahedrons, tetrahedrons, and arbitrary or custom elements in 3D.This procedure also derived expressions for checking the interelement compatibility of adjacent elements, which enables incorporating newly formulated elements within multielement meshes.The procedure was incorporated into GUI software, ShapeGen3D, to facilitate its implementation.

Discussion
This development reduces implementation effort in formulating and verifying arbitrary 3D elements for FEM applications.It provides FEM users the possibility to employ elements beyond those available in commercial FEM software, FEM libraries, and openaccess codes.The proposed procedure facilitates the formulation of higher-order elements, which can be shown to solve challenging problems with fewer degrees of freedom in some applications.The software was developed as part of a NASA-sponsored space technology research institute, Resilient Extra-Terrestrial Habitat Institute (RETHi), which aims to propel space exploration forward by developing new knowledge, technologies, and techniques.ShapeGen3D was used to develop the 21-node brick elements presented in this paper to capture meteoroid impact loading on a space habitat model.The model was embedded within a Modular coupled virtual testbed (MCVT) [42,43] that integrates several sub-systems of a space habitat to evaluate their interdependence and the effectiveness of safety controls when a habitat is exposed to external disturbances.ShapeGen3D allowed developers to easily adjust the elements of the structural model until the functional requirements were met.A 43-node brick element was also developed with ShapeGen3D to mimic wave propagation due to meteoroid impact and assist with generating a structural health monitoring algorithm for detecting damage.
Upon publication, this development is expected to benefit the computational mechanics community beyond the RETHi group.The higher-order and transition element formulation that the software can generate will allow researchers to employ elements that can increase the accuracy and reduce the computational time when dealing with simulations that involve contact, mesh transition regions, and Lagrangian solid dynamics.Furthermore, ShapeGen3D can become an educational tool to illustrate three-dimensional FEM concepts.

Conclusions
This work presents the development of a generalized procedure for formulating nodal coordinates, shape functions, and integration quadrature of higher-order hexahedrons, tetrahedrons, and arbitrary or custom elements in 3D.This procedure also derived expressions for checking the interelement compatibility of adjacent elements, which enables incorporating newly formulated elements within multielement meshes.The procedure was incorporated into GUI software, ShapeGen3D, to facilitate its implementation.In this paper, several elements were formulated to verify the accuracy of the implementation.ShapeGen3D allows users to verify that the nodal shape conditions are met and checks the element's compatibility with other elements through a graphical interface.This advanced feature was illustrated through the 21-hexahedron element, where compatibility with a 20-node brick element was checked.Additionally, 43-node elements were formulated to capture the behavior of a Germain-Lagrange plate and establish a transition between Lagrangian elements of different order.Finally, a fifth-order tetrahedron element was generated to simulate guided waves.Thus, this work can have potential impacts on the accuracy and efficiency of FEM runs, as it can be used to generate elements that minimize the number of elements and degrees of freedoms used in traditional analyses, develop transition elements that connect fine to coarse mesh regions, and generate high-order spectral elements that accurately capture waves propagating through solid bodies.The numerical results illustrate that the generalized procedure captured in ShapeGen3D can significantly benefit research and development in FEM and numerical methods.

Figure 3 .
Figure 3. Algorithm for shape function determination.

Figure 3 .
Figure 3. Algorithm for shape function determination.

Figure 5 .
Figure 5. Elements 1 (a) and 2 (b) in x coordinate.Bold black rectangle indicates surface p and q.

Figure 6 . 2 is
Figure 6.Element 2 after transformation in  coordinate.According to Equation (34), element 2 needs to be rotated such that the    aligns with −    .The resulting transformation of the  coordinates is  =  1   2  (38) Equation (38) facilitates the formation of an interelement boundary between  and  .The rotation matrix,   , about the    axis, and the coordinate translation is performed to satisfy Equation (35) as  = (  1     2

Figure 6 .
Figure 6.Element 2 after transformation in X coordinate.

Figure 9 .
Figure 9. Overview of the options in Segments 1 and 2: (a) Hexahedron or Brick elements, (b) Tetrahedron element, (c) Arbitrary or Custom elements, and (d) Check interelement compatibility.

Figure 9 .
Figure 9. Overview of the options in Segments 1 and 2: (a) Hexahedron or Brick elements, (b) Tetrahedron element, (c) Arbitrary or Custom elements, and (d) Check interelement compatibility.

Figure 12 .
Figure 12.Nodal shape function for (a) Node 5 and (b) Node 21 evaluated in the space of the 21-node hexahedron element.

Figure 13 .
Figure 13.Shape function comparison between the 20 and 21-node brick elements.

Figure 13 .
Figure 13.Shape function comparison between the 20 and 21-node brick elements.

Figure 17 .
Figure 17.Evaluation of shape functions over the volume of the element for (a) 70 of a higherorder hexahedron and (b) 50 of a higher-order tetrahedron element.

Figure 17 .
Figure 17.Evaluation of shape functions over the volume of the element for (a) N70 of a higher-order hexahedron and (b) N50 of a higher-order tetrahedron element.

Table 2 .
Parametric coordinates of element surface.

Table 2 .
Parametric coordinates of element surface.

Table 4 .
Output files to Matlab workspace.(n x i , n y i , n z i ) where n x i : x coordinate of node i .

Table 5 .
Factors of interpolation function of 75 and 43 node elements.

Table 5 .
Factors of interpolation function of 75 and 43 node elements.