Automatic 3D modeling by combining SBFEM and transfinite element shape functions

The scaled boundary finite element method (SBFEM) has recently been employed as an efficient means to model three-dimensional structures, in particular when the geometry is provided as a voxel-based image. To this end, an octree decomposition of the computational domain is deployed and each cubic cell is treated as an SBFEM subdomain. The surfaces of each subdomain are discretized in the finite element sense. We improve on this idea by combining the semi-analytical concept of the SBFEM with certain transition elements on the subdomains' surfaces. Thus, we avoid the triangulation of surfaces employed in previous works and consequently reduce the number of surface elements and degrees of freedom. In addition, these discretizations allow coupling elements of arbitrary order such that local p-refinement can be achieved straightforwardly.


Introduction
The scaled boundary finite element method (SBFEM) is a semi-analytical technique -loosely based on finite elements -that involves a discretization of boundaries of computational (sub-)domains.Roughly speaking, this method aims at transforming a partial differential equation (PDE) in two or three spatial coordinates into a set of ordinary differential equations (ODE) in one coordinate by discretizing all but this one coordinate.In order to apply this idea effectively, a particular coordinate system is usually chosen in which one coordinate ξ points from the origin 1 to the boundary while the remaining one or two coordinate(s) (η, ζ) describe a parametrization of the boundary.
The SBFEM was originally developed to model large and unbounded domains in the context of soil-structure interaction and was inspired by concepts such as similarity [65], cloning [66], as well as the thin layer method [31,33].A detailed description of the underlying formulation can be found in the early papers [57,67,58], as well as the recent textbook [56].Later, it has been noticed that the SBFEM can be employed as a means to construct arbitrary star-convex2 elements [47,45,8,46].Such polygonal elements enable the application of rather flexible meshing procedures, compared to conventional finite elements which are usually restricted to triangular/quadrilateral and tetrahedral/hexahedral shapes.A particular variant of such a meshing paradigm consists in the use of domain decompositions of the quad-/octree type.The most common of this class of decompositions consists of square-/cube-shaped cells (see Fig. 1a for a minimal example).While, from a geometrical viewpoint, these meshes only require quadrilateral/hexahedral subdomains, the SBFEM concept al-lows each side to be divided into an arbitrary number of surface elements on the boundary.Hence, coupling subdomains of different size is straightforward.This idea has been used successfully for the efficient meshing of complex geometries, where flexible local refinement is desirable [43,44,28,25].Such a meshing paradigm is particularly useful for image-based analyses, i.e., in applications where the geometry (the distribution of material parameters, etc.) is provided in a pixel graphics format.A number of approaches have been developed to mesh images, most of which can be divided into two categories (see [68] for a comprehensive literature review): The first type includes meshing based on boundary detection for each region in the image, where a region is assumed to be homogeneous.This can be straightforwardly achieved using techniques such as a marching cubes algorithm [41].Subsequently, these boundaries are utilized to mesh the region inside using conventional meshing algorithms, such as the advancing front [40] or Delaunay tessellation [9].The second type involves techniques that mesh the image directly.These techniques are generally quicker due to the straightforward nature of the meshing process.One very basic example is the pixel-based approach [34], where each pixel is modeled as a quadrilateral/hexahedral finite element.While such a meshing approach is absolutely trivial and fast, it leads to an unnecessarily large number of degrees of freedom.In order to improve the efficiency of the meshing, other direct meshing techniques can be deployed such as the aforementioned quad-/octree meshing structure [68].
Quadtree meshing is performed by recursively dividing an image matrix into four equal-sized cells at a time.A criterion of homogeneity is established based on the difference in the maximum and minimum color intensity within a cell.If this difference exceeds a userdefined threshold, the cell is divided.The maximum and minimum cell size allowed are also specified by the user.As a result, the meshing scheme adaptively refines the regions around different interfaces.Simultaneously, the scheme keeps relatively larger cells within each region where the material can be assumed to be homogeneous.However, since cells of different sizes exist in quadtree meshes, compatibility issues are encountered when using conventional finite elements.This issue was circumvented by utilizing the scaled boundary finite element method (SBFEM) for quadtree meshes, as this method allows the use of arbitrary polyhedral subdomains.A detailed explanation of image-based analysis using the SBFEM can be found in [28,53].
Notwithstanding the success of this method in twodimensional image-based analyses, the situation in threedimensional modeling is somewhat different.The quad-tree-based domain decomposition can be extended to 3D without much further effort [53].In this case, we divide a cubic cell into 8 smaller ones of equal size and refer to this technique as octree decomposition (Fig. 1b).From the viewpoint of the SBFEM, each cubic cell constitutes one subdomain, which implies that we have to discretize the surface of each cube in a finite element sense.Fig. 1b reveals that the surfaces of 3D octree decompositions basically consist of 2D quadtree meshes and hence exhibit the same compatibility issues discussed above.That is, due to the two-dimensional finite element grid that has to be used for the discretization of the boundary, hanging nodes are generated.In this example, the mesh on the front surface of the structure in Fig. 1b is identical to the two-dimensional case in Fig. 1a.In finite element applications, the hanging node problem can be circumvented by either using some kind of constraints (constraint equations, Lagrange multipliers, etc.) or triangular elements.
In previous works, this issue has been tackled by the second option, i.e., triangulating the surfaces of individual cubes to ensure compatibility with adjacent subdomains, see Fig. 1c.This approach has been applied successfully to solve problems related to acoustics [38] and fracture [52].However, this concept may be considered somewhat unsatisfactory as the triangular elements require additional degrees of freedom and are intrinsically less accurate compared to quadrilateral elements.Furthermore, the previous approach does not allow combining elements of different order within the same model, which is one of the major advantages exploited in 2D SBFEM models [28,26].
To improve on the previous formulation, we present in the current paper an approach to discretize the surfaces of octree-based SBFEM models without introducing additional nodes.This is achieved by employing a certain type of transition element that allows us to formulate shape functions on quadrilateral domains with an arbitrary subdivision of each edge.The application of such elements was inspired by the papers on the socalled 'pNh elements' [21,15,63], however, similar formulations exist and date back to the 1970s [18,19,16,3,6].In a recent publication, we presented a detailed review and generalization of these classes of transition elements [13].There, we refer to these general transition elements as xNy-elements to indicate that it is in principle possible to couple arbitrary element families.In the current paper we shall only provide a very brief summary of the formulation (Sect.2.3) and refer the reader to the pertinent literature.The purpose of this work is mainly to demonstrate that these transition elements can in fact be used to solve the hanging node problem in three-dimensional SBFEM models.The motivation for this approach is twofold: -We wish to avoid any triangulation of the surfaces in three-dimensional (octree-based) models.This goal can easily be achieved using said transition elements, since they allow splitting each edge into an arbitrary number of sections.Consequently, we avoid introducing additional elements that are not required by the topological decomposition.Especially when us-ing higher order interpolation, the reduction in degrees of freedom (DOFs) can be very significant.-We wish to allow different element orders within the same model, which is particularly useful in the case of inhomogeneous materials and wave propagation problems.Up to now, conventional Lagrange elements have been used on the surfaces of SBFEM subdomains, which required the element order to match between adjacent subdomains.Using the proposed approach generally allows connecting elements of arbitrary orders.

Theory
The overarching goal of this contribution is to combine the concepts of SBFEM, octree-based mesh generation, and transition elements.This combination results in an efficient and robust numerical tool which can be applied to a wide range of problems in engineering and physics.To this end, we provide a brief descriptions of the fundamental principles underlying the SBFEM in Sect.2.1, while octree-meshes are discussed in Sect.2.2.
Considering the formulation of transition elements, we sketch the basic ideas in Sect.2.3 and refer to the literature concerning transfinite elements for an in-depth derivation.

SBFEM
The formulation of the Scaled Boundary Finite Element Method is presented for the case of linear elasticity, governed by the following partial differential equation for the displacement field u = u(x, y): Here, ρ denotes the mass density and D is the elasticity matrix, which, for isotropic material with shear modulus G and Poisson's ratio ν, reads L is the differential operator which may be rewritten for later use as with b1 =   1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 By adjusting the differential operator, as well as the material parameters, other common linear partial differential equations are obtained and can be solved by the SBFEM in a similar fashion [2,38].
The SBFEM constitutes a semi-analytical approach, in the sense that the governing partial differential equation (in three coordinates) is discretized in two directions while remaining analytical in the third coordinate.In order to employ this idea effectively, a particular coordinate transformation is usually applied such that the two discretized directions (denoted as η, ζ) define a parametrization of a domain's boundary, while the analytical coordinate ξ describes the direction from the origin to the boundary. 3A few key steps of the procedure are summarized in this section for easier reference.The complete formulation including a detailed derivation of the basic coordinate transformation and the semi-discretization can be found in [57].For dynamic problems, we apply the solution procedure proposed in [55], which is based on a continued fraction expansion of the dynamic stiffness matrix.A concise summary of the procedures applied to dynamic problems is presented in [26].
In general, the coordinate transformation employed in the SBFEM is written as where x, ŷ, ẑ define a three-dimensional Cartesian coordinate system and x, y, z are the values of these Cartesian coordinates on the boundary of the domain, parametrized by the local coordinates η, ζ. 4 Without loss of generality, the origin of the coordinate system is chosen to be inside the domain.Furthermore, the domain is assumed to be star-convex.The boundary is usually interpolated using a set of two-dimensional shape functions N(η, ζ) (such as Lagrange shape functions or splines); thus, the coordinates and their derivatives can be written as and analogously for ŷ, ẑ.The vectors x n , y n , z n denote the nodal values of the coordinates (or, more generally, the coefficients of shape functions if the chosen basis is not node-based).However, the reader may note that, for the current contribution, an interpolation of the boundary is not required as all surfaces of the octree-meshes addressed here are squares, for which the coordinate transformation could easily be written in closed form.Nevertheless, we stick to this more general formulation for later reference.With these transformations, the Jacobian determinant is obtained as where J is the determinant of the matrix The differential operator L transformed into the local coordinates reads The displacement field is interpolated on the boundary as where we assume for simplicity that the same shape functions are employed as for the interpolation of the geometry (isoparametric elements).Consequently, the approximated strain-displacement relationship reads with Applying the method of weighted residuals to obtain a weak form of the governing equations in the coordinates η, ζ, leads to the ordinary matrix differential equation with the coefficient matrices where d denotes the number of dimensions (2 or 3).Many details on the properties and the available solution procedures of similar differential equations are discussed in [32].For the static case (ω = 0), the differential equation is usually re-written as with and Here, q n denotes the vector of nodal forces and I is the identity matrix.A set of eigenfunctions of the semidiscrete differential operator in Eq. ( 16) can be obtained by performing an eigenvalue decomposition of the matrix Z. Alternatively -to improve accuracy in cases where (near)-parallel eigenvectors exist -a blockdiagonal Schur factorization can be employed instead, which leads to a decomposition of the form Details are discussed in [54].In the above decomposition, Ψ 11 , Ψ 21 denote the modes that lead to finite displacements at the origin, and S 11 is the corresponding Schur block.The solution for the case of a bounded domain is then obtained as On the boundary (ξ = 1), we obtain which we use to compute a stiffness matrix such that In the dynamic case (ω > 0), the differential equation in displacements ( 14) is converted into one for the dynamic stiffness S(ω) on the boundary: Assuming an approximation of the form a mass matrix is obtained from the solution of the Lyapunov equation: To enhance the accuracy for higher frequencies, highorder stiffness and mass matrices can be computed of the form 0 , S 0 , ..., S where the matrices 1 are derived from a continued fraction expansion of the dynamic stiffness matrix and the matrices X (i) are introduced for preconditioning.The details of this approach can be found in [55,1].

Octree mesh
While the concept of the SBFEM can generally be applied to any star-convex subdomain, the use of the quadtree/octree decomposition has become particularly popular.As already discussed in the introduction, octreemeshes can be obtained very easily by starting with a coarse discretization consisting of few cubes and recursively dividing cubes into eight smaller cubes according to some criteria.These criteria may be based on whether the cube is intersected by a boundary or whether the variation of material parameters within the cube exceeds a predefined threshold.For further details, we refer again to [53].Also, the reader may note that such a meshing procedure will not lead to smooth boundaries of complex geometries due to the 'stair-case'-like approximation of interfaces.These geometry errors are acceptable in many applicationsparticularly in image-based analyses where the available initial geometry description is a pixel-based image.In other cases, it is desirable to improve the description of the boundaries by applying smoothing algorithms, cutting the cubes that are intersected by the boundary [39], or combining the octree decomposition with fictitious domain concepts [25].Such approaches are the subject of other publications and will not be discussed here in detail.In the current contribution, we restrict ourselves to pure octree-meshes since this is the part where our approach differs from previous works.
The proposed application of transition elements changes the way the meshes on the surfaces of the SBFEM subdomains are handled.For now, let us assume that the octree decomposition is 'balanced', which means that adjacent subdomains differ by at most one refinement level.Consequently, each edge of each cube is either split into two segments or not.Taking into account the symmetries under rotation, there exist only six distinct meshing patterns on the cube surfaces.These patterns are depicted in Fig. 2 and differ in the number of sides that are split in order to ensure conformal coupling to the adjacent subdomain.Note that the first and last pattern in each row all consist of standard quadrilateral elements; thus, in the proposed approach, there are only five different types of surface elements to be considered. 5As will be shown in the following sec-tions, we can construct shape functions with nodes at the desired positions.For comparison, Fig. 2 also displays the subdivision of the patterns using conventional triangular and quadrilateral elements as has been done in previous works.
The difference in both approaches becomes prominent when applying higher-order interpolations.As an example, Fig. 3 shows the same patterns with additional nodes using a polynomial order of 3. Note that the transition elements do not require internal nodes to achieve complete polynomials up to an order of 3.This is in contrast to the Lagrange elements employed previously.In that sense, the shape functions are similar to that of serendipity finite elements.Furthermore, it should also be noted that, when using transition elements, there is no need to restrict the decomposition to balanced meshes.In fact, subdomains of arbitrary sizes can be coupled straightforwardly.This would be rather tedious when using conventional elements, since adequate triangulations would have to be defined for each of the (arbitrarily large number of) different patterns.

Transition elements based on the xNy element concept
As already discussed in the previous sections, transition elements are being used in the context of finite element models to realize local mesh refinement.In the framework of the SBFEM, transition elements can be used to discretize each face of a hexahedral SBFE which is connected to a number of smaller SBFEs having the same or similar shape functions.This versatility would, in principle, allow us to connect different types of SBFEs based on various sets of shape functions, comprehensively discussed in Ref. [27].In this paper, our focus is, however, on high order shape functions based on Lagrange polynomials that are defined on a non-equidistant grid of points [30,48].Typically, a Gauss-Lobatto-Legendre (GLL) nodal distribution is chosen to generate the high order Lagrange shape functions.This particular set of nodes includes the element (interval) boundaries at ±1, while the interior nodes are the roots of the Lobatto polynomials of order p − 1 if shape functions of order p are constructed [10].

Projection: Fundamental idea
The idea of transition elements, as implemented in this paper, can be traced back to the works of Gordon and co-workers [16,17,19,20].In their articles, the transition element was referred to as transfinite element due to the special type of mapping that was applied.In the  remainder of the article, we will stick to the term transition element.It is worth mentioning that the same approach to deriving shape functions for transition elements can also be used to achieve an accurate approximation of the geometry in high order FEMs.In the p-FEM, this kind of geometry description is commonly referred to as blending function method [59,35,5] and is widely used to incorporate the exact geometry of the structure stemming, for example, from computer-aideddesign (CAD) software.More recent applications of the transition element technology can be found in Refs.[61,62,49,50,51] In the current contribution, the main goal is to couple SBFEs of different sizes, where each boundary of each element is discretized using two-dimensional quadrilateral elements.Therefore, the task is to develop a transition element that features piece-wise polynomial shape functions such that it can couple to multiple smaller elements.At this point, only the basic idea is sketched and we refer to the pertinent literature [16,17,19,20,61,62,49,50,51] and the references cited therein.In order to achieve conforming elements, we need to interpolate a function Ξ(ξ, η) over the reference domain To this end, we project the arbitrary bivariate function Ξ(ξ, η) onto a different (carefully selected) space of bivariate functions [19].Such a projection will be denoted as P[Ξ(ξ, η)] and is composed of two projectors P ξ [Ξ(ξ, η)], P η [Ξ(ξ, η)] that interpolate Ξ(ξ, η) along the local directions ξ and η.Additionally, a mixed projector P ξη [Ξ(ξ, η)] is introduced to ensure that redundant terms cancel out.Hence, the projection is written as where the mixed projection operator is defined as Note that the individual projection operators are both linear and idempotent where f (ξ, η) and g(ξ, η) are (continuous) bivariate functions and the subscript s ∈ {ξ, η}.

Projection: Definition of operators
The methodology sketched in the previous subsection is deployed to construct conformal finite element spaces ensuring that elements of different sizes can be coupled by introducing an approach to derive piece-wise polynomial shape functions.To achieve this goal, a simple definition of the projection operators suffices.As already introduced in the early works of Gordon and Hall [17,19], a linear Lagrange interpolation polynomial is appropriate to define the projectors as where ξ i and η i denote the nodes that are used to define the Lagrange polynomials.The functions ψ i (s) are identical to the linear Lagrange polynomials N 1 i (s) and in this context they are referred to as linear blending functions The definition of the mixed projection operator immediately follows from the successive application of Eqs. ( 33) and ( 34) to the function Ξ(ξ, η i ) Looking at Eq. ( 37), it is immediately clear that the mixed projector interpolates the function only at discrete nodal locations, whereas the individual projectors take the element edges into account.This procedure ensures that elements are compatible at the edges.
Based on the presented concept, we are now able to derive the shape functions for arbitrary transition elements.To this end, the function Ξ(ξ, η) is chosen such that it represents the piece-wise polynomial function that enables coupling elements of different sizes.Recall that the goal is to couple three-dimensional SBFEsubdomains of different sizes.In this case, it suffices to ensure the edges feature conformal shape functions as long as serendipity elements of degree p ≤ 3 are used.These elements do not include interior/bubble shape functions that are necessary to the polynomial completeness [11].If higher order elements or tensor product formulations [12] are deployed, we also have to account for the interior/bubble shape functions to ensure a conformal coupling of three-dimensional SBFEsubdomains.Note that the shape functions connected to the interior of the surface elements do not need to be included in the projection process, i.e., the bubble functions are identical to those of standard finite elements.
One possibility is to only add those hierarchic bubble modes that are known from the p-version of FEM [59] to complete the polynomial.

Numerical examples
In this section, we present the results of several numerical experiments that we performed in order to validate the proposed approach.We begin by analyzing simple geometries -namely a rectangular cuboid and a cube -and conduct standard patch tests, as well as a modal analysis.We then proceed to demonstrate the applicability of the approach to more complex geometries.
In all examples, we assume linear elastic behavior of the material and small deformations.We denote by Ω the computational domain and Γ u , Γ q are the parts of the boundary where Dirichlet and Neumann boundary conditions are applied, respectively.Hence, the general problem statement may be written as where u Γ and σ Γ denote the boundary conditions and L is the differential operator as given in Eq. (3).

Static analyses and patch tests
We begin by conducting linear and higher order patch tests to check the validity of the proposed approach and evaluate the rate of convergence under optimal conditions.These tests are performed for a cuboid with a width and height of a = b = 2 and a length of L = 4, see Fig. 4. Hence, the computational domain is given as In all tests, Dirichlet boundary conditions are applied to the two faces indicated in Fig. 4, while all other surfaces are traction-free.The material is isotropic and homogeneous with the following properties: Young's modulus: Despite the fact that the material is homogeneous, we choose to divide the domain into two regions, where different element sizes are employed.Thus, transition elements are incorporated at the interface between the two regions (see Fig. 5 for examples).The element order is varied between p = 1 and p = 3.We perform hrefinement by dividing each subdomain into eight in each refinement step.The element size h is defined as the edge length of the largest element in the mesh.To assess the accuracy of the computed results, we evaluate the L 2 norm of the relative error in displacements with respect to an analytical solution.We study the following three cases: Uniaxial tension We assume a state of uniaxial tension, such that the exact solution of this problem is and consequently, the following Dirichlet boundary conditions are applied: Hence, the numerical errors should be negligible, as long as the elements on the surface of each subdomain are capable of representing a linear variation of the displacement field exactly.Figure 5a shows the computed errors when performing h-refinement with varying element order.The results show that the proposed approach passes the linear patch test with error levels in the order of 10 −13 .
Beam bending As a 2 nd order patch test, we analyze the cuboid (with the same properties as before) in a state of pure bending (with a unit radius of curvature).The analytical solution of the displacement field is given in [60] (Chapter 10) as Hence, for the case of ν = 0, our Dirichlet boundary conditions read Results are presented in Fig. 5c.Since the exact solution is a quadratic function of x, y, z, the errors are negligible when using an element order of p ≥ 2. On the other hand, if an element order of p = 1 is chosen, the results are not exact.In this case, the numerical convergence rate -computed based on the last two points of the graphs -is obtained with a value of 2.3.
Cantilever beam Finally, we perform a 3 rd order patch test by applying boundary conditions that result in a cubic variation of the displacement field.The analytical reference solution of the displacement field ) is based on a cantilever beam which is weakly fixed on one end and subject to a shear force F applied at the other end (at z = 4 in our example).The analytical solution can be found in [4]:  According to this analytical solution, Dirichlet boundary conditions are applied at the two ends of the cantilever as Again, all other surfaces are traction-free.As expected, the error is negligible when using an element order of p = 3, see Fig. 5e.For linear and quadratic elements, we obtain a numerical convergence rate of 1.7 and 3.3, respectively.
Comparison against triangulation In addition to the above validation against analytical solutions, we compare the proposed approach to a previous implementation, in which surfaces are subdivided into triangular and rectangular elements rather than using transition elements (see the explanations in Sect.2.2).Such a discretization is depicted in Fig. 6a and can be compared directly with the corresponding mesh in Fig. 5f.For the sake of brevity, we present here only the results for the  cantilever beam.A comparison of the error in displacements using both approaches is presented in Fig. 6b.
The results obtained based on the proposed discretization (denoted as 'xNy' in the figure) are of course identical to the ones in Fig. 5e and are included here to facilitate the comparison.It can be seen that the application of triangulated surfaces ('tri') 6 results in very similar error levels for the same element size.Only when using an element order of p = 3 we observe a significant increase of the error levels -particularly for large element sizes -when using triangular elements.Plotting the same errors against the number of degrees of freedom (Fig. 6), the triangulated mesh obviously requires more DOFs to achieve the same error level as the transition elements.This is due to the fact that the quadtree patterns lead to a larger number of interior nodes when applying the triangulation (cf.Figs. 2 and 3).While in the case of linear elements (in the current example) only a few additional nodes are required at the centroids of triangulated surfaces, for higher-order elements the difference in the number of DOFs can be quite significant.
For instance, when comparing the finest discretization with p = 3, the numbers of DOFs differ by a factor of approximately 2.6 (292,539 vs. 110,715 DOFs).

Modal analysis -cube
We continue by performing a modal analysis in order to validate the proposed approach for dynamic problems.The computational domain is a cube of width 8 The material parameters are chosen as Young's modulus: E = 1 Poisson's ratio: ν = 0 Mass density: ρ = 1 All surfaces are traction-free.The mesh is defined by (8/h) 3 cubic subdomains, where h denotes the subdomains' side length.To ensure that transition elements are present in the mesh, an SBFE-subdomain at one of the corners of the cube is divided into eight, see Fig. 7 for examples of the mesh.A numerical reference solution has been computed by utilizing the conventional spectral element method on a uniform grid of 8   plotted in Fig. 7 with respect to the element size h for different element orders p.When using p = 1, 2, 3, the computed numerical convergence rate is obtained as 1.8, 3.8, and 4.5, respectively.Fig. 8: Modal analysis of a cube -relative error of the eigenfrequencies.

Modal analysis -Crane tower
To demonstrate the applicability of the proposed approach to more complex geometries, the model of a crane tower as depicted in Fig. 9  The total height of the model is 22.5 m.We apply fixed boundary conditions at the bottom of the base, as indicated in Fig. 9.The mesh has been created automatically by means of an octree decomposition of the voxelbased model, leading to 7158 subdomains.As can be inferred from Fig. 9, the decomposition in this case leads to subdomains of four different sizes and therefore, the side length of the largest element is 8 times larger than that of the smallest elements.A rapid and consistent transition from large to small element sizes is easily achieved by the octree meshing algorithm.Figure 10 shows the mode shapes of the first four modes and lists the corresponding eigenfrequencies.We found that to compute these modes, it suffices to employ quadratic elements on the surfaces of each subdomain, which results in a total of 615,960 DOFs.

Structure under self-weight loading
As a final example, we study the behavior of the structure depicted in Fig. 11 under the influence of selfweight loading.The model of a castle is based on a sample provided in the software MagicaVoxel [14] which we slightly modified and placed on a base of a different material.We scaled the model such that the total height of the structure is 33 m.The base is a homogeneous layer with a thickness of 8 m.The material parameters     The acceleration of gravity is assumed as g = −9.81m/s 2 .As a proof of concept, this example represents a situation where it can be useful to apply a different element order to the two distinct materials and thus better capture the larger deformations in the softer material.This feature will be particularly interesting for more complicated problems, such as wave propagation through layered soils, which has been addressed previously in two dimensions [28].In those cases, it can also be useful to adjust the element order not only based on the material properties but also on the size the subdomain.
In the current rather simple example, we found it sufficient to use elements of order 3 in the base material, while the comparably rigid material of the castle is discretized using linear elements only.The resulting mesh is depicted in Fig. 11a.The deformed geometry due to self-weight can be seen in Fig. 11b, where the colors indicate absolute values (magnitude) of displacement.

Conclusion
The combination of the three-dimensional SBFEM with transition elements allows a consistent discretization of an octree decomposition.Each surface of each subdomain is discretized by quadrilateral elements.Hence, we avoid any unnecessary subdivision into smaller surface elements as had been done in previous publications.The numerical results demonstrate that the computational models created using the proposed technique pass the linear as well as higher order patch tests.Compared to the previous meshing paradigm (involving triangulation), there is no loss of accuracy due to the transition elements, while the number of surface elements and degrees of freedom is reduced.We also demonstrated that the proposed approach allows coupling elements of different interpolation orders straightforwardly.Thus, we are now able to implement local p-refinement, which could previously only be exploited in two-dimensional SBFEM models.

Fig. 1 :
Fig. 1: Discretization in the SBFEM: (a) In 2D, a quadtree-mesh is used to discretize the computational domain straightforwardly.(b) In 3D, surfaces need to be discretized, which feature a structure similar to the quadtree meshes in 2D.(c) In previous work, these surfaces have been triangulated using standard finite elements.Using the transition elements discussed in this paper, the surface meshes can be handled directly without further subdivision (see part (b) of the figure).

Fig. 2 :
Fig.2: Mesh patterns on the surfaces of a balanced octree decomposition.Transition elements allow us to use quadrilateral elements in all cases (top), while conventional elements require a subdivision of the surfaces to treat hanging nodes (bottom).

Fig. 3 :
Fig. 3: Mesh patterns for the case of third order interpolation: Transition elements (top) and Lagrange elements (bottom).
Error vs. degrees of freedom.
is analyzed.The model is based on a stereolithography (STL) file obtained from the online repository Thingiverse [29]. 7From the STLfile, we generated a voxel-based model with the help of the online converter Voxelizer [64].We slightly modified the design -mainly to remove unconstrained parts as well as to assign two different colors that correspond to different materials in the context of the image-based analysis.The modifications have been made using the software MagicaVoxel [14].The material parameters of the tower and the base (as indicated by the different colors) are assumed as material 1 material 2 Young's modulus E: 70 GPa 20 GPa Poisson's ratio ν: 0.35 0.3 mass density ρ: 2.7 kg/m 3 2.4 kg/m 3

Fig. 10 :
Fig. 10: Modal analysis of a crane tower -Mode shapes and eigenfrequencies of the first four modes.
(a) Exemplary mesh.(b) Displacement field due to self-weight.

Table 1 :
Modal analysis of a cube -Reference solution for the first ten nonzero eigenfrequencies.

Table 2 :
Eigenfrequencies of the crane tower, computed using varying element order p.