A higher-order conformal decomposition finite element method for plane B-rep geometries

Computers and Structures 214 (2019) 15–27

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A higher-order conformal decomposition finite element method for plane B-rep geometries J.W. Stanford, T.P. Fries Institute of Structural Analysis, Graz University of Technology, Lessingstr. 25/II, 8010 Graz, Austria

a r t i c l e

i n f o

Article history: Received 10 August 2018 Accepted 13 December 2018 Available online 29 December 2018 Keywords: Higher-order finite elements Level-set method Fictitious domain method Interface capturing NURBS B-rep geometry

a b s t r a c t A higher-order accurate conformal decomposition finite element method (CDFEM) is proposed, using an explicit boundary description of the domain. A level-set concept is introduced which converts the explicit B-rep geometry based on NURBSs to an implicit representation. This concept also applies to higher-order fictitious domain methods (FDMs) which are often based on implicit descriptions. In the proposed CDFEM, an arbitrary, non-conforming background mesh is introduced and manipulated to capture non-smooth features of the boundary such as corners and to ensure shape-regular elements. Elements that are cut by the boundary are decomposed into conforming sub-elements. The final mesh composed by higher-order Lagrange elements enables higher-order accurate approximations in the context of the FEM. Special emphasis is on the fact that NURBSs feature reduced continuity between knot spans which must be properly considered in the mesh generation, not only at corners. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction In recent years there have been efforts to tighten the link between design (i.e., CAD-software) and analysis (i.e., FEMsoftware). Despite high efforts from industry and the research community, many open topics remain. Many difficulties are related to the different requirements of the geometry description, which makes the direct use of the geometry description of a domain for an analysis, as provided by a typical CAD-tool, rather non-trivial, a fact already recognized by the CAGD-community [1,2]. In general, to use CAD-geometry within numerical simulations, there are three alternatives: The first one would be to directly use the geometry—or the involved NURBS basis functions—for the approximation as in Isogeometric Analysis [3]. This approach has the advantage that results from the numerical simulation can be easily transferred back into the CAD-model (e.g., iteratively as in shape-optimization) and no creation of a separate mesh is needed. However, such an approach strongly relies on a geometry parametrization which is directly suitable for analysis. Challenging are tensor-product patches with collapsed or trimmed edges [4]; or the mere fact that only a parametrization of the geometry’s boundary is available, while a typical FEM-simulation requires a discretization of the whole domain of interest.

E-mail address: [email protected] (J.W. Stanford) URL: http://www.ifb.tugraz.at (J.W. Stanford) https://doi.org/10.1016/j.compstruc.2018.12.006 0045-7949/Ó 2018 Elsevier Ltd. All rights reserved.

The second option is to avoid the reliance on the provided geometry parametrization by generating a body-fitted mesh, i.e., a mesh, whose edges are accurately aligned to match the domain of interest. Apart from accurately and robustly capturing the boundary, additional challenges arise for creating meshes suitable for higher-order methods as discussed in [5]. In general, to construct a conforming mesh, two major options exist [6–8]: One option is to use the parametrization of the boundary to create a triangulation of it and then use either advancing front techniques [9,10] or a Delaunay triangulation to discretize the bulk domain [11,12]. These methods generally produce unstructured meshes. To produce structured meshes, methods like Marching-Cubes [13] or space-trees [14–16] can be used. These first fill the bulk with regular elements and then intersect those elements with the domain boundary. Methods, combining those two approaches, as shown in [17] also exist. For the creation of higher-order elements, a common approach is to first create a linear mesh and then transfer the elements therein to higher-orders. This can be applied to both of solid geometry and B-rep geometry. Special attention needs to be paid to avoid elements with invalid mappings (i.e., where the minimal Jacobian is negative), as discussed in [18–20]. Apart from using the parametrization provided by the CAD geometry, one can also employ boundary value problems to produce higher-order elements, as shown for instance in [21,22] with the Winslow equation, or with linear elasticity in [23]. Recently, optimization-based methods to ensure valid elements were presented in [24,25]. An interesting approach, where the exact bound-

16

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

ary description is used to create conforming elements [26,27] also needs to be mentioned. A third option would be to skip the laborious task of mesh creation at all and use basis functions from a background mesh for approximation. This background mesh does not usually conform to the boundary, rather it only embeds the given geometry. In these methods—collectively known as embedded domain, fictitious domain or immersed boundary methods—the challenge of capturing the boundary is shifted from the meshing-step to the integration of the weak form and the application of boundary conditions. An example is the Finite Cell Method (FCM)—a representative of embedded domain methods commonly applied for solid mechanics—as originally proposed in [28] and recently reviewed in [29]. For the application of boundary conditions, the FCM uses a triangulation of the surface instead of a NURBS-based discretization [30]. For integration purposes, the domain is partitioned using a space-tree subdivision. Recently, however, improved integration schemes were introduced [31,32]. Immersed domain methods may also consider NURBS-based boundary descriptions exactly during integration, as for example shown in [33,34] using special NURBS enhanced elements [26]. An integration scheme based on standard Lagrangian polynomials applicable to fictitious domain methods featuring optimal higher-order accuracy was also presented in [35,36]. A hybrid of classical meshing and embedded domain methods is the conformal decomposition FEM, introduced in a low-order context in [37] and extended to higher orders in [38,39]. The basic idea is to employ an arbitrary, non-conforming background mesh at the start and decompose elements cut by the boundary into conforming sub-elements. Additional measures have to be employed to ensure the shape regularity of the elements and avoid illconditioning. The important consequence is that meshes of any order suitable for a classical p-FEM analysis may be generated fully automatic. The CDFEM has so far only been considered for implicit geometry descriptions, e.g., in [37–39]. We conclude that, on the one hand, the classical FEM relies on suitable meshes based on explicit boundary descriptions such as given, e.g, by NURBSs. On the other hand, one may typically associate higher-order accurate FDMs and the CDFEM with implicit boundary representations based on level-sets. Therefore, it can be highly useful to convert an explicit description to the implicit one. The main contribution of this work is to propose a level-set concept for B-rep geometries based on distance calculations. For each knot span, different level-set functions are introduced: One for the actual distance to the bounded knot span, one to the extended knot span, and one for each end point of the individual knot span. Based on this, a global level-set function may easily be generated. This procedure may even close undesired gaps in the explicit boundary representation often present in the context of trimming, see [40] and the extensive overview in [4]. Based on the generated level-set data, the automatic, higherorder meshing algorithm for the CDFEM falls into the following steps: A background mesh is employed and manipulated such that (i) element nodes are present at the knots where the continuity of the boundary is reduced and (ii) the resulting elements are shaperegular as already shown in [41]. Then, the zero-level-sets are automatically detected and meshed followed by a decomposition of the cut elements following [38,39]. Herein, the higher-order CDFEM is used for the first time for explicit boundary descriptions. However, we emphasize that many parts of the proposed meshing procedure—including the conversion from NURBS geometries to level-sets—are also relevant within the context of fictitious domain methods: for providing integration sub-cells in the cut background elements, for visualization, and possibly for the imposition of boundary conditions as shown in [42].

The remainder of this paper is outlined as follows: Section 2 defines explicit boundary representations based on NURBSs. The conversion to level-set functions for the implicit description is given in Section 3. In Section 4, the conformal decomposition, tailored to NURBS geometries, is outlined. Numerical examples are shown in Section 5 and show the success of the method. Finally, Section 6 concludes this work. 2. NURBS-based domain description 2.1. Non-Uniform Rational B-Splines (NURBS) A univariate NURBS, as defined in [43] is a piecewise smooth rational polynomial of degree p and is defined as

Pn i¼0 N i;p ðuÞ wi P i CðuÞ ¼ P ; n i¼0 N i;p ðuÞ wi

u0 < u < um

ð2:1Þ

with weights wi associated to control points P i and B-Spline basis functions N i;p ðuÞ. Also part of the definition are knots ui , which are contained in a non-decreasing (i.e., ui 6 uiþ1 ) sequence of U ¼ fu0 ; . . . ; um g, denoted as the knot vector. Unless stated otherwise, herein open knot vectors are assumed; that is the entries u0 and um are repeated ðp þ 1Þ times. A multiplicity k of an entry in the knot vector reduces the curve’s continuity to C pk at those points. For knots that only appear once, this implies a continuity of C p1 at the knots. The regions between knots, denoted as knot spans, are smooth. For brevity we refer to a knot either as an entry ui in the knot vector or to its respective image Cðui Þ, depending on the context. Equivalently, depending on the context, we denote a knot span either as the closed interval ½uj ; ujþ1  or its image Cj ðtÞ ¼ CðuÞ with u 2 ½uj ; ujþ1 . Furthermore, to allow for concise formulation, the restrictions P 0 – P 1 and P n – P n1 are assumed. For simplicity, but without loss of generality, it is assumed that each knot span is being re-parameterized as a rational Beziér curve. This can be achieved by Beziér-splitting as described in [43], which results in new control points P j;i and weights wj;i . A rational Beziér curve is then defined as

Pn i¼0 Bi;p ðtÞ wj;i P j;i Cj ðtÞ ¼ P ; n i¼0 Bi;p ðtÞ wj;i

t 2 ½0; 1

ð2:2Þ

with the mapping

tðuÞ ¼

u  uj ujþ1  uj

uj – ujþ1

ð2:3Þ

and the Bernstein basis functions

Bi;p ðtÞ ¼

p!  ti  ð1  tÞpi : i !  ðp  iÞ!

ð2:4Þ

NURBS (and therefore Beziér curves) posses several properties, of which we make use of the following:  Strong convex hull property: The knot span Cj ðtÞ lies in the convex hull of the control points P jp . . . P j .  Ends are interpolatory: Cðu0 Þ ¼ P 0 and Cðum Þ ¼ P n .  End tangents are parallel to the first/last segment of the control polygon: C 0 ðu0 Þ ¼ k ðP 1  P 0 Þ and C 0 ðum Þ ¼ k ðP n  P n1 Þ, where k is an arbitrary scalar multiplier. 2.2. Domain boundary description Generally a boundary-representation geometry, is a geometry description where the solid domain is implied through given boundary information. In this work, the boundary information in

17

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

^ i can be determined based on a scalar proThe sign of /i and / duct with the curve’s tangent vector:



sðx; tH Þ ¼ hx  Ci ðtH Þ;

0

1

1

0



C0i ðtH Þi;

where h; i denotes the scalar product and Ci ðtH Þ is the closest point to x and C0i ðtH Þ is the tangent vector at that point. A visual comparison between / and /^ can be found in Fig. 2a and b. i

i

In addition, we define two auxiliary level-set functions with

qi;1 ðxÞ ¼ hx  P i;0 ; P i;1  P i;0 i Fig. 1. Exemplified boundary description with @ X ¼ fS1 ¼ fC 1 ¼ fC1 ; C2 ; C3 g; C 2 ¼ fC4 gg; S2 ¼ fC 3 ¼ fC5 ggg.

two dimensions is defined as followed: A boundary is a collection of a single or multiple sets of curves @ X ¼ fSi g. Each set of curves describes a closed loop and consists of one or more NURBS curves S ¼ fC j g, see Fig. 1 for an example. As noted previously, each curve can further be decomposed into knot spans. That is C j ¼ fCk g. The curves C i are to be defined in a way so that the left hand side of the curve’s tangent lies inside of the domain. Further it is noted that the hierarchy @ X ¼ fSi ¼ fC j ¼ fCk ggg can also be flattened to @ X ¼ fCk g. Additionally, it is assumed the boundary description fulfills the basic mathematical needs to describe solid geometries, as discussed in [44,45], which include: (i) non-self intersecting boundaries and (ii) no ‘‘dangling” edges.

3. Level-set-based domain description The task is to automatically generate meshes for the domain explicitly defined by NURBSs as described above. Therefore, it is useful to first convert the explicit boundary representation into an implicit description based on level-sets. 3.1. Implicit domain definition An implicit definition of the domain via means of the level-set function, a concept taken from [46,47], is given by

X ¼ fx j /ðxÞ > 0g;

ð3:1Þ

where the level-set function /ðxÞ is negative outside of X, positive inside X and zero on the boundary @ X. One such function to satisfy these requirements is the signed distance function

/ðxÞ ¼  min kx  xH k;

ð3:2Þ

8xH 2 @ X

where k  k denotes the Euclidean distance and the sign depends on whether x is in- or outside of X.

For each knot span Ci , a level-set function /i will be defined, similarly to Eq. (3.2) simply by replacing xH with Ci ðtÞ:

t 2 R:

ð3:3Þ

It needs to be noted, even though the knot span Ci ðtÞ is a bounded curve as per definition in Eq. (2.2), an evaluation outside is technically possible. In Eq. (3.3) we refer to this unbounded interpretation (t 2 R). In cases, where the signed distance to the bounded curve (t 2 ½0; 1) is sought, we refer to the level-set function as

/^i ðxÞ ¼  min kx  Ci ðtÞk

t 2 ½0; 1:

and

qi;2 ðxÞ ¼ hx  P i;n ; P i;n1  P i;n i;

ð3:6Þ

where P i;0 ; P i;n denote the first or last control point of Ci . These auxiliary level-sets give the distance to the tangent plane spanned by the first and last segment of the control polygon, as illustrated in Fig. 2c and d. 3.2.1. Numerical computation of the level-set function To compute the level-set functions in Eqs. (3.3) and (3.4), one can make use of the fact that if xH ¼ Ci ðtH Þ is the closest point on the curve to x, then the scalar product of the difference vector x  xH with the curve’s derivative at that point vanishes:

f i ðx; tÞ :¼ hx  Ci ðtÞ; C0i ðtÞi ¼ 0:

ð3:7Þ

Note, that this is a necessary, but not a sufficient condition. Solving this nonlinear equation for all t is done by a standard Newton–Raphson scheme. With ft g being the set of all t that fulfill the above equation, it is easy to compute the level-set function

/i ðxÞ ¼ sðx; t Þ  kx  Ci ðt Þk; where

t ¼ argmin kx  Ci ðtÞk t

t 2 ft g:

^ i as well. To reliably Note, that the same process also works for / find all roots of Eq. (3.7) with the Newton–Raphson method, good start values and corresponding ranges need to be determined a priori. Literature suggests a number of heuristic methods for this task. While [48–50] use properties involving the control polygon to limit the range of the start value, [51] provides a criterion for the uniqueness of a root within a knot span. Algorithms presented in the literature try to avoid computational cost by reducing the number of start values to a single one, however they all depend on x. So, hardly information can be re-used when computing /i ðxÞ for a large number of points. Therefore an alternative algorithm has been proposed by the authors in [52]. It pre-evaluates a knot span at select points and re-uses this information to compute a piecewise linear approximation of f i ðx; tÞ f i ðx; tÞ to obtain start values. An alternative to iterative methods would be to bring Ci ðtÞ into an implicit form Ci ðtÞ ! ci ðxÞ so that t is being eliminated, where ci ðxÞ is an algebraic function, and thus can be used directly as a level-set function. For more details on the subject of implicitation, the interested reader is referred to [53–55]. It should be noted, that cðxÞ does not in general feature the signed-distance property which complicates its use in the meshing algorithm discussed below. h

3.2. Signed distance to knot spans

/i ðxÞ ¼  min kx  Ci ðtÞk

ð3:5Þ

ð3:4Þ

3.3. Signed distance to domains composed of multiple NURBS curves In the previous subsection the creation of one level-set function /i corresponding to one knot span Ci was discussed. For a complete

18

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

Fig. 2. Different level-sets associated with the same knot span (red). The solid lines (gray) represent iso-values of values P 0, whereas the dashed ones represent values < 0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

domain description via level-sets, these functions may be combined into a single one with

/ðxÞ ¼ min fj/^i jg  sign ðminfj/^i jgÞ:

ð3:8Þ

In the rare case when multiple /i are of the same magnitude but different sign the computation of the sign needs special attention. This typically can occur near corners with an enclosing angle a < p=2. In such situations, a comparison of tangential vectors of each curve or curvatures will yield the desired sign as described in [52]. As an alternative, ray-casting techniques may also lead to the desired sign. Note that the resulting level-set function is not yet suitable for decomposing the background mesh, since its smoothness depends on the smoothness of the boundary itself, which is reduced at the knots as mentioned above.

4. Conformal decomposition of background mesh The NURBS geometry is embedded into a, typically rectangular, background mesh, so that X XBG . The task is to automatically generate a conforming mesh composed of higher-order elements. This is achieved by (i) first moving nodes of the background mesh onto the knots of the boundary representation and (ii) decompose cut elements into conforming sub-elements. The first step is of utmost importance for an accurate approximation due to continuity considerations, as shown later in Section 5.1. After the manipulation of the background mesh, the decomposition of a cut element may be realized with respect to only one of the level-set functions /i , individual for each knot span. The resulting elements are then combined into a mesh if they are inside the domain, which is being tested using the sign of /ðxÞ. Fig. 3 shows a schematic overview.

4.1. Background mesh Starting point is a background mesh with (bi-) linear elements. The size of XBG is determined based on the convex hull property of NURBS. As first step, the mesh is being locally refined such that all background elements do not contain more than one knot each. This allows the meshing of geometries with complex boundaries with a reasonable amount of elements. It is also useful and easily implemented to further refine background elements containing a knot (e.g. at a re-entrant corner) or based on the curvature of the knot spans to improve the geometry representation. 4.1.1. Consideration of knots A straightforward method to ensure sufficiently smooth boundaries is to move select corner nodes of the background mesh onto knots. To reduce the impact on the mesh quality induced by this operation, we found it useful to also move nodes within a given radius R in the same direction. The node movement vector for a given knot can be given as d ¼ d0  ð1  r=RÞ where d0 denotes the movement vector that moves a selected node on the given knot, and r the distance of a node to the knot. In case a node is affected by multiple movements, their average is taken. An alternative approach to move nodes of the background mesh without affecting element regularity, is to solve a pseudo elasticity problem where the necessary nodal movements are applied as Dirichlet boundary conditions. There has been extensive research e.g. in the field of fluid–structure interaction on how to best preserve the element regularity during this operation [56–60]. In this publication, an element is regular if its mapping from the reference element yields a non-negative Jacobian for every point in the reference domain. This condition is tested by evaluating the Jacobian on a fine sub-grid. Alternatively, one could also transfer the isoparametric mapping to a Bernstein basis. Then, by the

Fig. 3. Overview of the conformal decomposition process: (a) B-rep description, (b) level-set / to describe the domain, (c) smooth level-sets /i used for decomposition (gray: zero-level-set of /i , black: zero-level-set of /), (d) manipulated background mesh, (e) decomposed, conforming mesh.

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

non-negativeness of the Bernstein basis functions, it is sufficient to check the associated coefficients for positiveness, as shown for example in [19]. 4.1.2. Avoiding ill-shaped elements after decomposition Since edges of the background mesh may initially be arbitrarily close to the domain boundary, it has proven useful to move nodes that are close to the boundary away from it using a similar scaling function as shown above, see also [41]. The displacement vector for a node xi is



di ¼

dmax 

 dmax /ðxi Þ r/ðxi Þ ; dmin kr/ðxi Þk kr/ðxi Þk

ð4:1Þ

which will ensure all nodes closer than dmax have a minimum distance dmin to the domain boundary. In order to avoid undesired element mappings, di is applied in increments. In the examples shown in Section 5 the parameters were chosen as dmin ¼ 0:02 h and dmax ¼ 2 h. 4.2. Conformal decomposition 4.2.1. General procedure Conformal decomposition is the process where background elements intersecting the boundary are decomposed in a way that the resulting sub-elements are conforming to the boundary. The background-element-wise decomposition for a single, smooth level-set function follows closely the outline in [35,38,36] and is sketched as follows: 1. Evaluate the level-set function at all element nodes. 2. Determine if the background element is cut by comparing signs of the level-set values within the element. 3. Identify the topological cut situation as depicted in Fig. 4. 4. Depending on the cut situation, find the roots of the level-set function along the element edges. These intersection points span the interface element (the blue line in Fig. 5), which is a 1d-element that captures the boundary. 5. Determine the final position of the inner nodes of the interface element, by projecting them from the initially straight interface element onto the zero-level-set. 6. Define a mapping based on transfinite interpolation and the interface element for the sub-elements into the cut background element. For an overview of this procedure, see Fig. 5. In case no valid cut situation is found during step (3), the background element will be adaptively refined as needed.

19

4.2.2. Ensuring valid topological cut situations It is clear that, since the domain boundary may cut background elements in an arbitrary fashion, there may be an infinite number of topological cut situations (i.e., configurations in which @ X cuts the boundary of a background element). However, assuming sufficiently small elements with respect to the curvature of the boundary, the number of cut situations can be limited to the six shown in Fig. 4. These include all situations where @ X cuts each elementedge maximum once, the element boundary is crossed exactly twice and where the element boundary is crossed at nonadjacent corners. Other cut situations are denoted invalid and are resolved by adaptively refining the mesh, until only topological cut situations are encountered. As seen from Fig. 6, this does not introduce hanging nodes and therefore poses no difficulties in the following mesh generation. The described procedure uses standard Lagrangian elements with equidistant nodal spacing. In a related work [61] the authors have shown that Chen-Babuška elements can also be employed. For the sake of completeness it needs to be noted, that since element corner nodes are being moved onto knots, the definition of valid cut cases from [36,35,38] needs to be extended to allow intersection at non-adjacent element corners. Cut cases used in this work are listed in Fig. 4. 4.2.3. Consideration of knots With the background mesh prepared in a way that all knots coincide with element corners, it may seem sufficient to use the global, generally non-smooth level-set from Eq. (3.8) for the composition procedure. However for an accurate approximation, the level-set function—rather than only its zero-isoline—must be sufficiently smooth in the interior of an element. Therefore the smooth level-set functions /i ðxÞ of the knot spans must be used for this task instead. For obtaining regular elements it is desirable to only decompose elements which are actually cut by the knot span. Not every element cut by the unbounded /i ðxÞ ¼ 0 is also cut by the bounded Ci ðtÞ, as shown in Fig. 7. Therefore, to determine whether a given element is intersected by Ci ðtÞ, an intersection test is carried out by testing if

n o minf/i g  maxf/i g < 0 for i j qi;j > 0 ;

ð4:2Þ

where the index j denotes the corresponding end of the knot span. The ðÞi are numerical values obtained by evaluating the three levelsets /i ; qi;1 ; qi;2 on a sufficiently dense sub-grid within the elements. The limitation given for i essentially excludes values that would indicate a ‘‘false intersection” and prevents the decomposition of background elements. For a finite element analysis, the resulting set of conforming (sub)-elements needs to be transformed to a mesh with a concept of degrees of freedom and a connectivity matrix implying the (unique) element nodes. 5. Numerical results 5.1. Influence of boundary continuity To demonstrate the importance of moving nodes the background mesh onto the knots, or rather the consequence of not doing so, a series of numerical studies are shown. Based on the well-known problem of the plate with hole under uniaxial traction, a series of domains were derived where the arc-segment (which is usually C 1 between its endpoints) is replaced by splines which fea-

Fig. 4. Topological cut cases used in this work.

ture different continuity (C 0 ; C 1 ; C 2 ; C 1 ) on the boundary at a selected point as depicted in Fig. 8. A numerical study was then

20

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

Fig. 5. Overview of the general conformal decomposition process.

Fig. 6. Background mesh width edges cut twice and three times. (a) Initial background mesh, (b) Refined background mesh to resolve invalid cut situations, (c) decomposed higher-order mesh. Note that due to high-curvature, the mesh was refined further during decomposition.

Fig. 7. Illustration for the partial decomposition.

Fig. 8. Used B-rep geometries with C 0 ; C 1 ; C 2 ; C 1 -continuous point in the arc-like segment. The reduced continuities were achieved by successive knot insertion followed by a slight movement of a control point nearby.

conducted to determine the convergence behavior for geometry and boundary value problem approximations on those four different domains without considering this knot in the meshing procedure. That is, no node of the background mesh is moved to this

knot and influence on the accuracy is investigated. For comparison, a series of studies was conducted where the knot was properly considered. In all eight studies, the other five corners were properly considered. For the boundary value problem, an elastic shell was

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

chosen with Dirichlet boundary conditions prescribed on the whole boundary. In this and the subsequent studies, errors were measured in the following norms:

R   1 dx  Aref  eA ¼ X Aref eC ¼ k/ðxÞkL2 ðCÞ   R 1=2 e : r dx  U ref  X eU ¼ : U ref

21

ð5:1Þ ð5:2Þ ð5:3Þ

Fig. 9. Convergence results for the study of continuity and moving nodes onto knots. Each column in the grid represents the same C k variation, whereas each pair of rows represent the same measure of error ei . The dashed lines indicate the slope of p þ 1 to the corresponding color, with p being the polynomial degree of the basis functions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

22

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

That is, eA measures the error in the area, eC in the boundary representation, and eU in the energy norm, where Aref and U ref were determined trough an overkill solution using an extremely fine mesh for each of the four geometries. It can be shown, as for instance in [62], that all three measures should converge with Oðh

pþ1

Þ. The results of the study are summarized in Fig. 9. From 0

the convergence plots it is seen that for the C continuous case, 2

the rate of convergence is bounded to a slope of h for all three

error measures in the log–log plots shown, which is only optimal for p ¼ 1. A similar behavior can be observed for the C 1 and C 2 continuous situations: For the C 1 case only the series with p  2 converge optimal. All series with higher polynomial degree do not display an increased rate of convergence. Similar observations can be made in the C 2 case, where no series converges faster than the p ¼ 3 series. In contrast the C 1 case and the options where the kink has been properly considered do not show this behavior of

Fig. 10. The spanner geometry.

Fig. 11. Results of the convergence study for the spanner.

Fig. 12. Geometries used for the plate studies.

23

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

limited convergence rates. This study clearly shows the importance of considering all points with reduced continuity between the knot spans (not only the corners) in the proper meshing. 5.2. Linear elastic spanner To test the applicability of the proposed method for engineering examples, the spanner shown in Fig. 10a was studied as a plane stress problem in linear elasticity. For the definition of the boundary value problem in terms of partial differential equations, the reader is referred to Section A.1 or to standard literature such as [63,64]. The geometry was modeled using three NURBS curves of degree p ¼ 2, with the geometry described in Table 1 in the appendix. The material parameters are as follows: Shell thickness of t ¼ 5 mm, Young’s modulus E ¼ 2:1 105 N=mm2 , Poisson ratio m ¼ 0:3, equally distributed line load with a resultant of F ¼ 500 N. The exact area was computed as Aref ¼ 20008:0118044056510 mm2 and the total stored elastic energy was determined as U ref ¼ 16:51041118  1 108 Nmm. Boundary conditions were as depicted in Fig. 10b and the resulting deformations and von Mieses stresses in Fig. 10c. The results in Fig. 11 show the convergence of the measures eA and eU as defined in Eqs. (5.1) and (5.3) as well as the estimated condition number j of the stiffness matrix. The condition number 2

grows mostly with Oðh Þ as expected [62], while the measured errors converge optimal for orders up to p ¼ 4. The deviation from this rate in the range of larger element sizes indicates (a) a higher number of and (b) less shape-optimal elements due to adaptive refinements. See Fig. 10b for such a mesh. Nevertheless, the results show that this adaption-driven behavior vanishes for smaller h and optimal rates of convergence can be achieved until the measured error seems to be dominated by the effects of machine precision. 5.3. Plates The following study demonstrates the applicability of our method to Kirchhoff–Love plates and Reissner–Mindlin plates. These are more advanced test cases because the condition number

4

of Kirchhoff–Love plates is known to scale with Oðh Þ [62] as the governing equations involve derivatives up to 4th order, leading to higher requirements on the mesh quality. For the Kirchhoff–Love plate, a mixed formulation as in [65] is used, where the primary unknowns are the deflection w, the bending and twisting moments mx ; my ; mxy , thus avoiding the requirement of C 1 -continuous shape functions. As for the Reissner–Mindlin theory, the implementation is based on [66], yielding the deflection w and the rotations of the cross sections /x and /y as primary unknowns. For both formulations, the governing equations are found in the appendix in Section A.2. Two different geometries depicted in Figs. 12a and b are studied, each with both plate theories, resulting in a total number of four studies. The NURBS-curves for both geometries are documented in Tables 2 and 3 in the appendix respectively. Chosen material properties are: E ¼ 3:5 107 kN=m2 ; m ¼ 0:3; , ¼ 5=6 and the plate thickness is h ¼ 0:3 m for the Kirchhoff–Love formulation and h ¼ 0:5 m for the Reissner–Mindlin variant. The surface traction f z and Dirichlet boundary conditions were determined based on manufactured solutions given below. To facilitate the determination of approximation errors, manufactured solutions were applied to all four resulting test cases. For the Kirchhoff–Love plate, it is sufficient to define the deflection as in:

wðx; yÞ ¼ ½sinð0:2 xÞ þ 2 cosð0:2 x yÞ þ 1  105

ð5:4Þ

The other quantities required can be derived using the governing equations. Likewise, for the Reissner–Mindlin plate it is sufficient to define the cross-sectional rotations as in:

  /x ðx; yÞ ¼ lxx sin 2 lxx p  104   /y ðx; yÞ ¼ cos 3 lyy p  104

lx ¼ 5 ly ¼ 3

ð5:5Þ

Note, that the functions defined in Eqs. (5.4) and (5.5) are not compatible with each other. For both plate theories, the boundary conditions are applied appropriately. For geometries A and B, the coarsest mesh and the displacement fields of the two solutions

Fig. 13. Coarsest mesh and displacement fields for both manufactured solutions. As the used implementation does not support mixed meshes with simplex and quadrilateral elements, the meshes were converted to pure simplex meshes after decomposition.

24

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

Fig. 14. Displacement errors

ew;L2 and estimated condition number j of stiffness matrix for plate study.

are depicted in Fig. 13a and b respectively. The quantity measured in the conducted convergence study is the displacement error in the relative L2 -norm:

ew;L2 ¼

wðxÞ  wh ðxÞ L kwðxÞkL2 ðXÞ

2 ðXÞ

:

ð5:6Þ

From the results in Fig. 14 it can be seen that for the Kirchhoff– Love formulation, higher-order convergence rates are achieved. It can be argued, that the initially (for large h) optimal rates are deteriorated due to the high condition number, which for fourth-order 4

BVPs is known to scale with Oðh Þ [62]. In case of the Reissner– Mindlin plate, optimal rates of convergence are achieved for both 2

proposed method does not require a structured background mesh, can also be applied to fictitious domain methods and is able to achieve optimal higher-order rates of convergence.

Appendix A. Strong and weak formulations of the boundary value problems under consideration A.1. Linear elasticity The equilibrium of forces for the case of linear elasticity is given as

r  r ¼ f ;

ðA:1Þ

geometries and the condition number is increasing with Oðh Þ. The smooth increase of condition numbers in all four test cases, indicates that the proposed method is able to produce meshes suitable for the analysis of plates.

with vector of body forces f and stress tensor r, which can be related to the strain tensor e trough a linear-elastic material law

6. Conclusions

r ¼ k trðeÞ I þ 2le;

In this work, a framework to compute level-set data from B-rep geometries parameterized with NURBS curves was presented and applied to a higher-order conformal decomposition finite element method. In the first step the NURBS-based description is converted into a level-set representation consisting of multiple smooth levelsets using the signed distance function. Each smooth level-set corresponds to a smooth segment of the boundary (i.e., a knot span or Beziér segment). The geometry is embedded into a higher-order background mesh. With the help of these level-sets, during the conformal decomposition, the cut background elements are replaced by conforming sub-elements which are accurately aligned to the domain boundary. In numerical studies it has been shown that optimal rates of convergence can be achieved even for non-trivial geometries. The

ðA:2Þ

where k; l are the Lamé material parameters, which for a planestress state (as used in the examples) can be computed from the Young’s modulus E and Poisson ratio m:

l¼

E ; 2 ð1 þ mÞ

k¼

Em : ð1 þ mÞ ð1  2mÞ

ðA:3Þ

The linear strain tensor can be computed from the displacements u as

e ¼ 1=2½ru þ ðruÞ| :

ðA:4Þ

Following the standard procedure of the FEM, this is converted to a weak from and, upon introducing test and trial functions, the discrete weak form following the Galerkin method is obtained as: Find uh 2 H1 such that

25

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

R X

rðuh Þ : eðwh ÞdX ¼

R

h

Xw

 f dX þ

^ u ¼u r  n ¼ ^t

R CN

wh tdC

on CD ; on CN ;

h

ux ¼ w;x þ cxz

8wh 2 H10 ; ðA:5Þ

uy ¼ w;y þ cyz :

ðA:10Þ

Hence, a relationship between shear forces and shear deformation can be given as

"

#

"

w;x  ux

#

^ are prescribed displacements on the Dirichlet boundary CD where u and ^t are prescribed tractions on the Neumann boundary CN with the restrictions C ¼ CN [ CD and CN \ CD ¼ £. For further details, the reader is referred to standard literature such as [67].

with

A.2. Plate formulations

cS ¼

The governing equations for the plate models considered in Section 5.3 are given as follows: The equilibrium of (vertical) forces and moments can be stated as

where , is a shear correction factor (set to 5=6). The strong form can derived directly by substituting (A.10) and (A.11) into (A.6), which results in:

qx;x þ qy;y ¼ f z ; mxx;x þ mxy;y ¼ qx ;

ðA:6Þ

mxy;x þ myy;y ¼ qy ; where qx ; qy are the shear forces and mxx ; myy the bending moments in the respective directions. The torque moment is designated as mxy and f z is the vertical surface traction. Partial derivatives are written as ðÞ;k ¼ @ðÞ=@k. The rotation of cross-sections—lines originally perpendicular to the middle-surface—are denoted as ux ; uy and can be related to the moments as

2

2

3

3

ux;x þ muy;y 6 7 6 7 m ¼ 4 myy 5 ¼ cB 4 mux;x þ uy;y 5; 1m mxy ðux;y þ uy;x Þ 2 mxx

ðA:7Þ

3

Eh 12 ð1  m2 Þ

A.2.1. Kirchhoff-Plate plate formulation For the Kirchhoff–Love plate, lines orthogonal to the middle surface remain straight and orthogonal. Hence, the rotation of these lines can be directly related to the change in plate deflection w

uy ¼ w;y :

ðA:8Þ

The final system of equations is gained by subsequently plugging Eq. (A.8) into Eqs. (A.7) and (A.6). After re-arranging this can be written as

mxx;xx þ 2 mxy;xy þ myy;yy ¼ f z ; 1 w;xx þ c ð1 m2 Þ ðmxx  m myy Þ ¼ 0; B

w;yy þ c

B

1 ð1m2 Þ

¼ cS

w;y  uy

;

ðA:11Þ

Eh, ; 2 ð1 þ mÞ

cS ðw;xx  ux;x Þ þ cS ðw;yy  uy;xy Þ ¼ f z ; cB ðux;xx þ muy;xy Þ  cB cB

1m 2

1m 2

ðux;yy þ uy;xy Þ  cS ðw;x  ux Þ ¼ 0;

ðux;xy þ uy;xx Þ  cB ðuy;yy þ mux;xy Þ  cS ðw;y  /y Þ ¼ 0; ðA:12Þ

containing w; /x ; /y as unknowns. The weak form is then obtained following the standard FEM procedure based on the Galerkin method. Appendix B. NURBS geometries used in the numerical studies Tables 1–3.

P 1;i

where E is the Young’s modulus, m the Poisson’s ratio and h the plate thickness.

ux ¼ w;x ;

qy

Table 1 Control points, weights and knot vectors for the three quadratic NURBS curves defining the spanner. Coordinates are in mm.

with

cB ¼

q¼

qx

ðmyy  m mxx Þ ¼ 0;

ðA:9Þ

1 w;xy þ cB ð1 mÞ mxy ¼ 0;

containing w and the entries in the symmetric moment tensor m as the four unknowns. The continuous weak form is arrived by standard methods and the discrete weak form by the Galerkin method. It is noted, that this results in a mixed formulation keeping the deflection w and m as primary unknowns. Further details can be found in [65,68]. A.2.2. Reissner–Mindlin plate formulation For the Reissner–Mindlin formulation, shear deformations add up to the rotation lines orthogonal to the mid-surface

xi

yi

55.5 21.0 18.7 27.9 32.9 45.2 68.8 146.7 224.6 255.1 283.9 301.5 316.6 331.8 331.8 331.8 316.7 301.5 284.0 255.3 225.0 144.6 64.3 37.4 19.0 15.0 10.5 0.7 11.6 53.4 71.1

21.1 45.6 30.8 27.4 19.1 1.1 1.1 1.1 1.1 1.1 8.9 15.0 4.3 6.5 25.0 43.6 54.4 65.2 59.1 49.1 49.1 49.1 49.1 49.1 68.8 73.1 76.7 84.5 82.4 75.4 36.9

wi

ui

1.00 0.88 1.00 0.94 1.00 0.87 1.00 1.00 1.00 0.99 1.00 0.89 1.00 0.89 1.00 0.89 1.00 0.89 1.00 0.99 1.00 1.00 1.00 0.92 1.00 1.00 1.00 0.91 1.00 0.88 1.00

0.0 0.0 0.0 77.9 77.9 96.7 96.7 139.8 139.8 295.6 295.6 356.0 356.0 390.3 390.3 424.6 424.6 458.9 458.9 493.2 493.2 553.4 553.4 714.0 714.0 764.9 764.9 776.4 776.4 800.0 800.0 877.9 877.9 877.9 (continued on next page)

26

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27

Table 1 (continued)

Table 2 (continued)

P 2;i

wi

xi

yi

71.1 48.1 25.1 18.0 12.4 9.9 0.0 2.5 9.5 32.5 55.5

36.9 43.0 49.2 51.1 46.4 27.5 0.0 6.9 8.8 14.9 21.1

P 3;i xi

xi 1.00 1.00 1.00 0.89 1.00 0.82 1.00 0.89 1.00 1.00 1.00

0.0 0.0 0.0 47.6 47.6 61.1 61.1 112.2 112.2 125.7 125.7 173.3 173.3 173.3

wi

ui

25.1 1.1 1.1 1.1 25.1 49.1 49.1 49.1 25.1

1.00 0.71 1.00 0.71 1.00 0.71 1.00 0.71 1.00

0.0 0.0 0.0 37.7 37.7 75.4 75.4 113.1 113.1 150.8 150.8 150.8

wi

12.1 9.5 6.4

1.8 0.0 0.0

1.00 0.95 1.00

xi

yi

14.5 11.7 10.4 13.1 14.5

1.8 3.8 1.8 0,07 2 1.8

P 2;i xi 6.4 0.3 0.3 0.3 0.3 4.3 4.3 9.7 13.4 18.9 15.7 18.9 16.1 12.1

wi

ui

1.00 1.00 1.00 1.00 1.00

0.0 0.0 3.3 5.7 9.0 11.3 11.3

wi

ui

P 1;i

wi

ui

1.00 1.00 1.00 1.00 1.00 1.00 1.00

0.0 0.0 0.0 0.0 22.5 22.5 22.5 45.0 45.0 45.0 45.0

wi

ui

1.00 0.96 1.00 0.96 1.00 0.96 1.00

0.0 0.0 0.0 20.0 20.0 39.1 39.1 59.4 59.4 59.4

wi

ui

1.00 0.71 1.00 0.71 1.00 0.71 1.00 0.71 1.00

0.0 0.0 0.0 12.6 12.6 25.1 25.1 37.7 37.7 50.3 50.3 50.3

yi 13.5 10.2 8.6 0.0 8.6 10.2 13.5

6.1 8.4 11.2 24.0 11.2 8.4 6.1

P 2;i xi

yi

6.1 4.2 13.3 16.0 13.5 4.3 6.1

13.5 14.3 9.5 0.1 9.4 14.3 13.5

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

0.0 0.0 6.0 10.0 11.0 17.4 21.4 27.7 33.1 39.6 46.2 51.9 55.9 60.9 65.7 65.7

P 3;i xi 8.0 8.0 0.0 8.0 8.0 8.0 0.0 8.0 8.0

yi 0.0 0.0 4.0 5.0 11.3 11.3 5.0 5.0 10.3 6.5 1.8 0.5 4.6 1.8

0.0 0.0 0.0 6.1 6.1 6.1

Table 3 Three NURBSs for plate geometry B (of respective orders p ¼ 3; 2; 2). Coordinates are in m.

Table 2 Three NURBSs for plate geometry A (of respective orders p ¼ 1; 1; 2). Coordinates are in m. P 1;i

ui

yi

xi

yi

319.8 319.8 295.8 271.8 271.8 271.8 295.8 319.8 319.8

P 3;i

ui

yi 0.0 8.0 8.0 8.0 0.0 8.0 8.0 8.0 0.0

References [1] Shephard MS, Beall MW, O’Bara RM, Webster BE. Toward simulation-based design. Finite Elem Anal Des 2004;40(12):1575–98. [2] Riesenfeld RF, Haimes R, Cohen E. Initiating a CAD renaissance: Multidisciplinary analysis driven design: Framework for a new generation of advanced computational design, engineering and manufacturing environments. Comput Methods Appl Mech Engrg 2015;284:1054–72. [3] Hughes T, Cottrell J, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Engrg 2005;194(39-41):4135–95.

J.W. Stanford, T.P. Fries / Computers and Structures 214 (2019) 15–27 [4] Marussig B, Hughes TJR. A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects. Archive Comp Mech Engrg 2017:1–69. [5] Dey S, Shephard MS, Flaherty JE. Geometry representation issues associated with p-version finite element computations. Comput Methods Appl Mech Engrg 1997;150(1):39–55. [6] Frey P, George P-L. Mesh generation: Application to finite elements. Chichester: John Wiley & Sons; 2008. [7] Lo DS. Finite element mesh generation. Boca Raton, FL: CRC Press; 2014. [8] Thompson NW Joe F, Soni Bharat, editors. Handbook of grid generation. Boca Raton, FL: CRC Press; 1998. [9] Peraire J, Vahdati M, Morgan K, Zienkiewicz O. Adaptive remeshing for compressible flow computations. J Comput Phys 1987;72(2):449–66. [10] Löhner R, Parikh P. Generation of three-dimensional unstructured grids by the advancing-front method. Internat J Numer Methods Fluids 1988;8 (10):1135–49. [11] Frederick CO, Wong YC, Edge FW. Two-dimensional automatic mesh generation for structural analysis. Internat J Numer Methods Engrg 1970;2(1):133–44. [12] Watson DF. Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes. Comput J 1981;24(2):167–72. [13] Lorensen WE, Cline HE. Marching cubes: A high resolution 3d surface construction algorithm. ACM SIGGRAPH Comput Graphics 1987;21(4):163–9. [14] Yerry MA, Shephard MS. Automatic three-dimensional mesh generation by the modified-octree technique. Internat J Numer Methods Engrg 1984;20 (11):1965–90. [15] Baehmann PL, Wittchen SL, Shephard MS, Grice KR, Yerry MA. Robust, geometrically based, automatic two-dimensional mesh generation. Internat J Numer Methods Engrg 1987;24(6):1043–78. [16] Shephard MS, Georges MK. Automatic three-dimensional mesh generation by the finite octree technique. Internat J Numer Methods Engrg 1991;32 (4):709–49. [17] Schroeder WJ, Shephard MS. A combined octree/Delaunay method for fully automatic 3-d mesh generation. Internat J Numer Methods Engrg 1990;29 (1):37–55. [18] Luo X, Shephard MS, Remacle J-F, O’Bara RM, Beall MW, Szabó BA, et al. Pversion mesh generation issues. In: 11th International meshing roundtable. p. 343–54. [19] Luo X-J, Shephard MS, O’Bara RM, Nastasia R, Beall MW. Automatic p-version mesh generation for curved domains. Eng Comput 2004;20(3):273–85. [20] Sherwin SJ, Peiró J. Mesh generation in curvilinear domains using high-order elements. Internat J Numer Methods Engrg 2002;53(1):207–23. [21] Knupp PM. Winslow smoothing on two-dimensional unstructured meshes. Eng Comput 1999;15(3):263–8. [22] Fortunato M, Persson P-O. High-order unstructured curved mesh generation using the Winslow equations. J Comput Phys 2016;307:1–14. [23] Xie ZQ, Sevilla R, Hassan O, Morgan K. The generation of arbitrary order curved meshes for 3d finite element analysis. Comput Mech 2012;51(3):361–74. [24] Toulorge T, Geuzaine C, Remacle J-F, Lambrechts J. Robust untangling of curvilinear meshes. J Comput Phys 2013;254:8–26. [25] Gargallo-Peiró A, Roca X, Sarrate J. A surface mesh smoothing and untangling method independent of the CAD parameterization. Comput Mech 2013;53 (4):587–609. [26] Sevilla R, Fernández-Méndez S, Huerta A. NURBS-enhanced finite element method (NEFEM). Internat Numer Methods Engrg 2008;76(1):56–83. [27] Sevilla R, Rees L, Hassan O. Mesh generation for the 2d NURBS-enhanced finite element method. In: ECCOMAS Congress 2016 - Proceedings of the 7th European congress on computational methods in applied sciences and engineering, vol. 1, 2016, p. 444–53. [28] Parvizian J, Düster A, Rank E. Finite cell method. Comput Mech 2007;41 (1):121–33. [29] Schillinger D, Ruess M. The finite cell method: A review in the context of higher-order structural analysis of CAD and image-based geometric models. Arch Comput Methods Eng 2014;22(3):391–455. [30] Düster A, Parvizian J, Yang Z, Rank E. The finite cell method for threedimensional problems of solid mechanics. Comput Methods Appl Mech Engrg 2008;197(45):3768–82. [31] Kudela L, Zander N, Kollmannsberger S, Rank E. Smart octrees: Accurately integrating discontinuous functions in 3D. Comput Methods Appl Mech Engrg 2016;306:406–26. [32] Stavrev A, Nguyen LH, Shen R, Varduhn V, Behr M, Elgeti S, et al. Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method. Comput Methods Appl Mech Engrg 2016;310:646–73. [33] Marco O, Sevilla R, Zhang Y, Ródenas JJ, Tur M. Exact 3d boundary representation in finite element analysis based on Cartesian grids independent of the geometry. Internat J Numer Methods Engrg 2015;103 (6):445–68. [34] Marco O, Ródenas JJ, Navarro-Jiménez JM, Tur M. Robust h -adaptive meshing strategy considering exact arbitrary CAD geometries in a Cartesian grid framework. Comput Struct 2017;193:87–109.

27

[35] Fries T, Omerovic´ S. Higher-order accurate integration of implicit geometries. Internat J Numer Methods Engrg 2016;106(5):323–71. [36] Fries TP, Omerovic´ S, Schöllhammer D, Steidl J. Higher-order meshing of implicit geometries - part I: Integration and interpolation in cut elements. Comput Methods Appl Mech Engrg 2017;313:759–84. [37] Noble DR, Newren EP, Lechman JB. A conformal decomposition finite element method for modeling stationary fluid interface problems. Internat J Numer Methods Fluids 2010;63(6):725–42. [38] Omerovic´ S, Fries T-P. Conformal higher-order remeshing schemes for implicitly defined interface problems. Internat J Numer Methods Engrg 2016;109(6):763–89. [39] Fries T-P. Higher-order conformal decomposition FEM (CDFEM). Comput Methods Appl Mech Engrg 2018;328:75–98. [40] Sederberg TW, Finnigan GT, Li X, Lin H, Ipson H. Watertight trimmed NURBS. In: ACM SIGGRAPH 2008 Papers, SIGGRAPH ’08. New York, NY, USA: ACM; 2008. p. 79:1–8. [41] Fries T, Schöllhammer D. Higher-order meshing of implicit geometries, Part II: Approximations on manifolds. Comput Methods Appl Mech Engrg 2017;326:270–97. [42] Stanford JW, Fries T-P. A higher-order conformal decomposition FEM for NURBS-based geometries. In: Proceedings of the 6th European conference on computational mechanics (ECCM 6) and 7th European conference on computational fluid dynamics (ECFD 7), Glasgow, UK. [43] Piegl L, Tiller W. The NURBS book. Berlin: Springer; 1997. [44] Sakkalis T, Shen G, Patrikalakis N. Representational validity of boundary representation models. Comput Aided Des 2000;32:719–26. [45] Requicha AA. Representations of rigid solid objects. In: Computer aided design modelling, systems engineering, CAD-systems. Berlin: Springer; 1980. p. 1–78. [46] Osher S, Fedkiw RP. Level set methods: An overview and some recent results. J Comput Phys 2001;169(2):463–502. [47] Sethian JA. A fast marching level set method for monotonically advancing fronts. Proc Nat Acad Sci 1996;93(4):1591–5. [48] Ma YL, Hewitt W. Point inversion and projection for NURBS curve and surface: Control polygon approach. Comput Aided Geom Des 2003;20(2):79–99. [49] Oh Y-T, Kim Y-J, Lee J, Kim M-S, Elber G. Efficient point-projection to freeform curves and surfaces. Comput Aided Geom Des 2012;29(5):242–54. [50] Selimovic I. Improved algorithms for the projection of points on NURBS curves and surfaces. Comput Aided Geom Des 2006;23(5):439–45. [51] Chen X-D, Yong J-H, Wang G, Paul J-C, Xu G. Computing the minimum distance between a point and a NURBS curve. Comput-Aided Des 2008;40(10):1051–4. [52] Steidl JW, Fries T-P. Automatic conformal decomposition of elements cut by NURBS. In: Proceedings of the VII ECCOMAS congress. [53] Hoffmann C. Implicit curves and surfaces in CAGD. IEEE Comput Graph Appl 1993;13(1):79–88. [54] Hoffmann CM. Geometric and solid modelling: An introduction. San Francisco, CA: Morgan Kaufman; 1989. [55] Sederberg T, Anderson D, Goldman R. Implicit representation of parametric curves and surfaces. Comput Vision Graphics Image Process 1984;28(1):72–84. [56] Masud A, Hughes TJR. A space-time galerkin/least-squares finite element formulation of the navier-stokes equations for moving domain problems. Comput Methods Appl Mech Engrg 1997;146(1):91–126. [57] Brackbill JU, Saltzman JS. Adaptive zoning for singular problems in two dimensions. J Comput Phys 1982;46(3):342–68. [58] Stein K, Tezduyar T, Benney R. Mesh moving techniques for fluid-structure interactions with large displacements. J Appl Mech, ASME 2003;70(1):58. [59] Stein K, Tezduyar TE, Benney R. Automatic mesh update with the solidextension mesh moving technique. Comput Methods Appl Mech Engrg 2004;193(21-22):2019–32. [60] Xu Z, Accorsi M. Finite element mesh update methods for fluid–structure interaction simulations. Finite Elem Anal Des 2004;40(9-10):1259–69. [61] Fries T-P. Higher-order accurate integration for cut elements with ChenBabuška nodes. In: SEMA SIMAI Springer Series. Berlin: Springer; 2016. p. 245–69. [62] Brenner S, Scott R. The mathematical theory of finite element methods. In: Texts in applied mathematics. Berlin: Springer; 2008. [63] Hughes TJR. The finite element method: Linear static and dynamic finite element analysis. New York, NY: Dover; 2000. [64] Zienkiewicz OC, Taylor RL, Fox D. The finite element method for solid and structural mechanics. 7th ed. Oxford: Butterworth-Heinemann; 2014. [65] Boffi D, Brezzi F, Fortin M. Mixed finite element methods and applications. Springer Series in Computational Mathematics, vol. 44. Berlin: Springer; 2013. [66] Bathe K-J, Dvorkin EN. A four-node plate bending element based on Mindlin/ Reissner plate theory and a mixed interpolation. Internat J Numer Methods Engrg 1985;21(2):367–83. [67] Zienkiewicz OC, Taylor RL, Zhu JZ. The finite element method: Its basis and fundamentals. Oxford: Butterworth-Heinemann; 2013. [68] Ciarlet PG, Raviart PA. A mixed finite element method for the biharmonic equation. In: Mathematical aspects of finite elements in partial differential equations. Amsterdam: Elsevier; 1974. p. 125–45.