Loibl, M.

Marussig, B.

Teixeira, Guilherme Henrique

Scipedia

The 9th European Congress on Computational Methods in Applied Sciences and EngineeringECCOMAS Congress 20243–7 June 2024, Lisboa, PortugalCOMPARISON OF INTEGRATION METHODS FOR CUTELEMENTSG. H. TEIXEIRA1, M. LOIBL2, B. MARUSSIG11 Institute of Applied MechanicsGraz University of Technology,Technikerstraße 4/II, 8010, Graz, Austriae-mail: teixeira@tugraz.at2 Institute of Mechanics and Statics,University of the Bundeswehr Munich.Werner-Heisenberg-Weg 39, 85577, Neubiberg, Germanye-mail: michael.loibl@unibw.deKey words: Integration, Cut Elements, IsoGeometric AnalysisSummary.Using an interface inserted in a background mesh is an alternative way of constructing acomplex geometrical shape with a relative low meshing efforts. However, this process mayrequire special treatment of elements cut by the interface. Our study focuses on comparingthe integration of cut elements defined by implicit and parametric curves. We investigate theefficiency and robustness of open-source tools such as Algoim [5](a library for quadrature onimplicitly defined geometries) and Ginkgo [2](a library for isogeometric analysis on Booleanoperations with a parametric description) with numerical examples computing the area definedby the interface and benchmarks for 2D elasticity problem using the open-source code GeoPDEs[7]. It is concluded that none of the two interface descriptions is preferable with respect to thequality of the integration. Thus, the choice of the interface type depends only on the studiedproblem and the available curve description, but not on the numerical aspects of the integration.1 INTRODUCTIONModelling the geometry and discretising the domain of a problem with a mesh is an importantstep in the process of solving any boundary value problem. Regarding the construction of ageometric representation, a common way to do it is to define a parametric shape and map thisgeometry into the physical domain. However, finding a suitable parametrization for complexgeometries is in general not an easy task. One way to mitigate this problem is to define thegeometry through a background mesh, where an interface specifies a region to be visible. Thisapproach allows us to analyse any geometry by simply selecting a part of the background meshto define the domain, as can be seen in the scheme presented in Figure 1. This concept is presentin methods such as the Immersed Boundary Method [8] or the Finite Cell Method [6].In the process of selecting the domain, some elements are cut by this interface. In consequence,among several challenges, the numerical integration of cut elements is an important point to beverified [4]. In particular, we would like to investigate, how the quadrature points are distributedover the cut elements, and we would like to observe how this is affected by the type of theG. H. TEIXEIRA, M. LOIBL, B. MARUSSIGFigure 1: Scheme of domain definition with an interfaceinterface used. Thereby, we compare parametric and implicit interface descriptions for cutelements. Depending on the type of interface an appropriate quadrature scheme is applied todefine the points and weights in each cut element. In particular, in this work, for parametriccurve interfaces, a reparametrisation of the element is performed and the quadrature is obtainedas in [2]. Whereas for implicit curve interfaces, an approximation of the curve by a Bézier curveis carried out and the corresponding points and weights are defined as in [5].2 INTERFACE REPRESENTATIONS2.1 Parametric curvesA parametric curve is defined by the map from the parameter domain to the physical domain,where each point in the parameter domain has its own image in the physical domain. Béziercurves, B-splines, and NURBS are examples of parametric curves. In the comparisons presentedhere, NURBS curves C(ξ) are used to create the interface, which can be defined by the linearcombination of NURBS basis functions Rj,p of degree p scaled by control points ciC(ξ) =n∑j=1Rj,p(ξ)ci. (1)At least one edge of a cut element is defined by this interface. And for such elements where theinterface is defined by a NURBS curve or portion of it, we employ the open-source code Ginkgo[2] to compute the quadrature points and weights. In this approach, the first step to obtainthe quadrature points is defined by the reparametrisation of the element cut by the interface.However, sometimes the reparametrisation of an element with only one B-spline surface is notpossible. So, the cut element is divided into several tiles, where each of them can be reconstructedusing B-spline surfaces. Once this is done, the quadrature points from a reference element aremapped to each one of these small parts of the reparametrised element.2.2 Implicit curvesAn implicit curve is defined by a locus of points for which an equation takes a specific value.The points with the same value form isolines and, for convenience, the zero set is usually chosenas the isoline that implicitly describes the curve. These points are defined to fulfil the conditionf(x, y) = 0. (2)The use of implicit curves to define the domain can be applied to solve topology optimisationproblems via the Level-set Method [1]. In order to compute the integral in implicitly defined2G. H. TEIXEIRA, M. LOIBL, B. MARUSSIG0 0.2 0.4 0.6 0.8 100.20.40.60.810 0.2 0.4 0.6 0.8 100.20.40.60.810 0.2 0.4 0.6 0.8 100.20.40.60.81Figure 2: Procedure for obtaining the quadrature pointsdomains, we apply the open-source code Algoim [5] to determine quadrature points and weights.In this approach, the implicit curve is approximated by a Bézier curve and then the algorithmis used to obtain the quadrature rule for domains implicit defined by multivariate polynomials.The quadrature points and weights are obtained by dividing the domain in pieces where in eachone of them the appropriate quadrature scheme is applied for one-dimensional integrals [5].3 METHODOLOGYThe process of defining the quadrature rule on the cut elements begins with the definition ofa background mesh and is followed by the insertion of an interface that creates the domain ofinterest. The difference between the approaches considered here is related to how these interfacesare defined for an identical geometry. Once the interface has been defined, the elements areidentified as active and non-active. The non-active elements are neglected in the computationand the active elements are divided into two subgroups cut and uncut elements. For the uncutelements, a Gaussian quadrature is defined in a reference element and then mapped to thephysical element, while for the cut elements the appropriate approach is applied depending onthe type of interface. Figure 2 describes the steps to obtain the quadrature rule. In a pre-existingbackground mesh, an interface defines the domain where, depending on the type of interface, thequadrature points are provided using the appropriate open-source code (Algoim [5] or Ginkgo[2]), which in sequence allows the numerical integration to be performed by simply evaluatingthe function at the quadrature points and scaling by the corresponding quadrature weights. Inorder to compare the different tools, first we consider the following verification settings:• Computation of the area enclosed by an interface for different refinement levels; and• computation of the area enclosed by an interface for a fixed refinement level while thegeometry moves.With these two settings we can verify not only the accuracy but also the robustness of thetwo approaches comparing the relative error with respect to the size of mesh and the relativeerror with respect to a given step. Once the behaviour of the tools is observed, we investigatehow the definition of the interface could affect the solution of 2D elasticity problems, where abackground mesh is defined using B-splines. For this, the open-source code for IsogeometricAnalysis GeoPDEs [7] is used. In this situation, two different scenarios of cut curves are set3G. H. TEIXEIRA, M. LOIBL, B. MARUSSIG0 0.2 0.4 0.6 0.8 100.20.40.60.81a)0 0.2 0.4 0.6 0.8 100.20.40.60.81b)Implicit curve Parametric curveFigure 3: Quadrature points: a) circle example; b) semicircle exampleup, one interface with no boundary conditions applied to it and one interface with boundaryconditions applied to it.4 NUMERICAL RESULTS4.1 Area computationOur first comparison study considers simply the computation of the area enclosed by aninterface. In this test, we verify the accuracy by computing the relative error for a fixed interfacewhile refining the background mesh. The first example is defined by a circle with radius R = 0.2and centered at (0.5, 0.5). In this example, the geometry is fixed using a single interface and themesh is refined by h-refinement. Figure 3a shows the shape and quadrature points for a mesh ofsize h = 0.25, which results in 36 quadrature points for both interface types. While the relative10−2 10−1 10010−1210−1010−810−610−410−2Size of elementRelativeErrorParametric curveImplicit curve9363610832461211162340936361082524689001836a)10−2 10−1 10010−1310−1010−710−410−1Size of elementRelativeErrorParametric curveImplicit curve4545361082976481350268299541622164689721944b)Figure 4: Convergence study: a) circle example; b) semicircle example4G. H. TEIXEIRA, M. LOIBL, B. MARUSSIG0 0.2 0.4 0.6 0.8 100.20.40.60.81a)0 0.2 0.4 0.6 0.8 100.20.40.60.81b)Implicit curve Parametric curveFigure 5: Quatrature points: a) Line ; b) Triangleerror with respect to the size of the elements is shown in Figure 4a. The second example is asemicircle, which is built using two interfaces. In particular, the interface using implicit curves isdefined by the intersection of a straight line y = −x+0.3536 with a circle with radius R = 0.25and centered at (0.1768,0.1768). The mesh is refined by h-refinement and Figure 3b shows thequadrature points obtained for the same element size as the circle example. Figure 4b showsthe relative error and the number of quadrature points needed for each size of mesh.4.2 Moving interfacesAnother aspect verified in this part of the work is the robustness of the approaches. For this,we look at different situations where it might be difficult to deal with the integration. This isdone by setting a mesh of size h = 0.25 and moving an interface around, trying to find positionsor situations where the codes are crashing. We investigate two situations: The first one is defined0.02 0.04 0.06 0.08 0.110−2010−1810−1610−1410−12StepL2ErrorParametric curveImplicit curvea)0.2 0.4 0.6 0.8 1 1.2 1.410−3010−2210−1410−6102StepL2ErrorParametric curveImplicit curveb)Figure 6: Robustness study: a) line; b) triangle5G. H. TEIXEIRA, M. LOIBL, B. MARUSSIGExact tractionSymmetrySymmetryTxRθrLLExacttractiona)0 0.2 0.4 0.6 0.8 100.20.40.60.81b)1.8e-010.30.40.50.60.70.80.91.01.1e+00∥ux∥c)Implicit curve Parametric curveFigure 7: Infinite plate with a hole: a) problem scheme; b) quadrature points; c) displacement magnitudeby a straight vertical line moving along the horizontal direction, as shown in Figure 5a. Whilethe second is given by a triangle in which the node (0, 0.5) rotates clockwise with respect to theorigin, while the nodes at (0, 0) and (0.5, 0) remain fixed, as shown in Figure 5b. Figure 6 showsthe relative error with respect to the given steps for both examples. For the line example, thestep represents the distance with respect to the first position of the line, which is at x = 0.5,and for the triangle example, the step represent the angle between the line that connects theorigin to the moving node and its original position, the line that connects the origin to the node(0, 0.5). As it is possible to see, none of the two approaches fails for any of the investigatedsteps.4.3 Elasticity benchmarksThe 2D elasticity problems are solved using the open-source code for Isogeometric AnalysisGeoPDEs [7] with quadrature points and weights provided by Algoim [5] and Ginkgo [2] .Thegeometries are defined by the intersection between a B-spline background mesh and one of thetwo types of an interface. And it is considered E = 1 and ν = 0.3.4.3.1 Infinite plate with a holeThe first example is the well-known benchmark of an infinite plate with a hole subjectedto a traction Tx = 10 at infinity, which can be modeled by a quarter of the geometry due tosymmetry. To do this, the background mesh Ω = [0, 1] × [0, 1] is cut by a circle with radiusR = 0.25, as shown in Figure 7a, and exact traction is applied at the top and right side of themesh, while symmetry boundary conditions are applied at the bottom and left side of the mesh.The analytical solution [3] is given asux =Tx ·Ri8µ[rRiA+ 2Rir(A+B)− 2R3ir3B](3)uy =Tx ·Ri8µ[rRiC + 2Rir(D + E)− 2R3ir3E], (4)6G. H. TEIXEIRA, M. LOIBL, B. MARUSSIGσyy = λ(εxx + εyy) + 2µεyyτxy = 2µεxyLHa)0 0.2 0.4 0.6 0.8 100.20.40.60.81b)5.6e-150.20.40.60.81.01.21.4e+00∥ux∥c)Implicit curve Parametric curveFigure 8: Plane square: a) problem scheme; b) quadrature points; c) displacement magnitudewhereA = (k + 1) cos(θ);B = cos(3θ);C = (k − 3) sin(θ)D = (1− k) sin(θ);E = sin(3θ); k = (3− 4ν);µ = E2(1 + ν).Figure 7b shows the position of the quadrature points on the cut element, obtained with the toolAlgoim [5], when the interface is defined by an implicit curve, and by the tool Ginkgo [2], whenthe interface is defined by a parametric curve. In Figure 7c, the magnitude of the component uxof the displacement vector u = [ux, uy]T is shown. The graph in Figure 9a shows the comparisonof the error with the size of the element, from which it is possible to see that solving the problemwith an element trimmed by an implicit curve gives as good results as when solving the problemdefining by a parametric interface. In both approaches we see almost the same accuracy.4.3.2 Square plateIn this example, we considered a square background mesh Ω = [0, 1] × [0, 1] trimmed atthe top, in the position y = 0.75, by a horizontal straight line, as shown in Figure 8a and 8b,which also shows the quadrature points obtained from both approaches and from the respectiveopen-source codes (Algoim [5] or Ginkgo [2]). The boundary conditions and the source term aremanufactured in a such way that the exact solution is given byux =sin(2πx) sin(2πy) (5)uy =sin(2πx) sin(2πy) (6)Figure 8c shows the magnitude of the component ux of the displacement field obtained fromthis problem. Figure 9b presents the graphs of convergence for different element sizes. Bothapproaches give similar accuracy, which is the same behavior identified on the previous exampleof an infinite plate with a hole. Hence, we can conclude that both approaches are well-suitedfor numerical approximation schemes.7G. H. TEIXEIRA, M. LOIBL, B. MARUSSIG10−110−610−510−410−310−2p+1Size of elementL2ErrorParametric curveImplicit curvea)10−110−510−410−310−2p+1Size of elementL2ErrorParametric curveImplicit curveb)Figure 9: Convergence study for p=2: a) infinite plate with a hole; b) square plate5 CONCLUSIONSIn this work, we compare two situations of cut elements that occur due to the intersection ofembedded interfaces in a background mesh. The two approaches differ in the interface definition,that is, either using parametric or implicit curves. We employ the open-source codes Algoim [5]and Ginkgo [2] to provide the quadrature points and weights over the cut elements. Numericalexamples of area computation are defined, from which we can observe similar behaviour in termsof convergence and the number of quadrature points used. Furthermore, using the open-sourcecode GeoPDEs, we apply both approaches to solve 2D elasticity problems with a backgroundmesh defined by B-splines. In this situation, we can compare the behaviour of these proceduresnot only when the integrand is equal to 1 (area computation problem) but also to observe theeffects of the definition of the interface in the solution of a different problem (2D elasticity).In all examples, we observe almost the same accuracy. Therefore, the examples show similaranswers for both approaches, which yield that the choice between the type of interface dependsmainly on the application of the problem, the advantage of using one rather than the other, andthe interface definition available but not on the numerical aspects of the integration.ACKNOWLEDGEMENTSThe work of Guilherme Henrique Teixeira is supported by the joint DFG/FWF CollaborativeResearch Centre CREATOR (CRC – TRR361/F90) at TU Darmstadt, TU Graz and JKULinz/Johann Radon Institute for Computational and Applied Mathematics, Linz.REFERENCES[1] G. Allaire, F. Jouve, and A. Toader. A level-set method for shape optimization. ComptesRendus. Mathématique, 334:1125–1130, 2002. doi: 10.1016/S1631-073X(02)02412-3.[2] P. Antolin, A. Buffa, and M. Martinelli. Isogeometric analysis on v-reps: First results.Computer Methods in Applied Mechanics and Engineering, 355:976–1002, 2019. doi: 10.8G. H. TEIXEIRA, M. LOIBL, B. MARUSSIG1016/j.cma.2019.07.015.[3] L. Coradello. Accurate isogeometric methods for trimmed shell structures. PhD thesis, Écolepolytechnique federale de Lausanne, 2021.[4] B. Marussig and T. J.R. Hughes. A review of trimming in isogeometric analysis: Challenges,data exchange and simulation aspects. Archives of Computational Methods in Engineering,25:1059–1127, 2018. doi: 10.1007/s11831-017-9220-9.[5] R. I. Saye. High-order quadrature on multi-component domains implicitly defined by multi-variate polynomials. Journal of Computational Physics, 448, 2022. doi: 10.1016/j.jcp.2021.110720.[6] D. Schillinger and M. Ruess. The finite cell method: A review in the context of higher-orderstructural analysis of cad and image-based geometric models. Archives of ComputationalMethods in Engineering, 22:391–455, 2015. doi: 10.1007/s11831-014-9115-y.[7] R. Vázquez. A new design for the implementation of isogeometric analysis in octave andmatlab: Geopdes 3.0. Computers and Mathematics with Applications, 72:523–554, 2016. doi:10.1016/j.camwa.2016.05.010.[8] X. Wei, B. Marussig, P. Antolin, and A. Buffa. Immersed boundary-conformal isogeometricmethod for linear elliptic problems. Computational Mechanics, 68:1385–1405, 2021. doi:10.1007/s00466-021-02074-6.9

Comparison of Integration Methods for Cut Elements

https://www.scipedia.com/wd/images/a/ad/Draft_Sanchez_Pinedo_702696023pap_732.pdf

Comparison of Integration Methods for Cut Elements

Abstract

Similar works

Full text

Available Versions

Scipedia