An improvement of the classical finite element method is proposed in [1], the NURBS-Enhanced Finite Element Method (NEFEM). It is able to exactly represent the geometry by means of the usual CAD description of the boundary with Non-Uniform
Rational B-Splines (NURBS). For elements not intersecting the boundary, a standard
finite element interpolation and numerical integration is used. But elements intersecting the NURBS boundary need a specifically designed piecewise polynomial interpolation and numerical integration. This document presents preliminary work on the 3D extension of NEFEM

Fernández Méndez, Sonia

Huerta, Antonio

Sevilla Cárdenas, Rubén

UPCommons. Portal del coneixement obert de la UPC

7th Workshop on Numerical Methods in AppliedScience and Engineering (NMASE 08)Vall de Nu´ria, 9 a 11 de enero de 2008c©LaCa`N, www.lacan-upc.es3D NURBS-ENHANCED FINITE ELEMENT METHODR. Sevilla, S. Ferna´ndez-Me´ndez and A. HuertaLaboratori de Ca`lcul Nume`ric (LaCa`N)Departament de Matema`tica Aplicada IIIUniversitat Polite`cnica de Catalunyae-mail: {ruben.sevilla, sonia.fernandez, antonio.huerta}@upc.eduweb: http://www-lacan.upc.eduKey words: NURBS, Finite Elements, CAD, exact geometry representation, high-orderisoparametric finite elementsAbstract. An improvement of the classical finite element method is proposed in [1],the NURBS-Enhanced Finite Element Method (NEFEM). It is able to exactly representthe geometry by means of the usual CAD description of the boundary with Non-UniformRational B-Splines (NURBS). For elements not intersecting the boundary, a standardfinite element interpolation and numerical integration is used. But elements intersectingthe NURBS boundary need a specifically designed piecewise polynomial interpolation andnumerical integration. This document presents preliminary work on the 3D extension ofNEFEM.1 INTRODUCTIONThe importance of the geometrical model in finite element (FE) simulations has recentlybeen pointed out by several authors, see [1, 2, 3] to name a few. In some applications,such as compressible flow problems, if a Discontinuous Galerkin (DG) formulation isadopted, see [4], an important loss of accuracy is observed when a linear approximationof the boundary is used, see [2]. Bassi and Rebay [2] show that, in the presence ofcurved boundaries, a meaningful high-order accurate solution can only be obtained if thecorresponding high-order approximation of the geometry is employed (i.e. isoparametricFE).Maxwell equations are also very sensitive to the quality of the boundary representation.Reference [5] studies the error induced by the approximation of curvilinear geometries withisoparametric elements. The 3D Maxwell equations are solved in a sphere with isopara-metric FE and with an exact mapping of the geometry. The exact mapping provides moreaccurate results with errors differing by an order in magnitud. Thus, in some applica-tions, an isoparametric representation of the geometry is far from providing an optimalnumerical solution for a given spatial discretization.1R. Sevilla, S. Ferna´ndez-Me´ndez and A. HuertaNon-Uniform Rational B-Splines (NURBS, see [6]) are widely used for geometry de-scription in CAD (Computer Aided Design). This fact has motivated new numericalmethodologies considering an exact representation of the computational domain withNURBS, such as the isogeometric analysis [3] and the NURBS-Enhanced Finite ElementMethod (NEFEM) [1]. The isogeometric analysis considers the same NURBS basis func-tions for both the description of the entire geometry and for the approximation of thesolution. NEFEM also considers an exact representation of the domain but it differs fromthe isogeometric analysis in two main points: the geometry is given by the NURBS de-scription of the boundary (i.e. the information usually provided by CAD), and standardFE polynomial interpolation is considered for the approximation of the solution. Thus,in the large majority of the domain —for elements not intersecting the boundary— astandard FE interpolation and numerical integration is used, preserving the computa-tional efficiency of classical FE techniques. Specifically designed piecewise polynomialinterpolation and numerical integration is only required for those FE along the NURBSboundary.In [1] NEFEM is presented in two dimensional simulations. Poisson and Maxwellproblems demonstrate the applicability of the proposed method in both a continuous anddiscontinuous Galerkin frameworks. When the quantities of interest are defined on curvedboundaries NEFEM is at least one order of magnitude more precise than the correspondingisoparametric FE. Moreover, the exact representation of the boundary allows to mesh thedomain independently of the geometric complexity of the boundary whereas classicalisoparametric finite elements need h-refinement to properly capture the geometry.This paper presents preliminary work on the 3D extension of NEFEM. Section 2 intro-duces the basic concepts on NURBS surfaces. In section 3 fundamentals of NEFEM arepresented in three dimensions. Special attention is paid to the interpolation and numer-ical integration in those elements affected by the NURBS description on the boundary.Finally, in section 4 a Poisson example demonstrate the applicability of the proposedmethod.2 BASIC CONCEPTS ON NURBS SURFACESIn this section NURBS curves are briefly introduced in order to present NURBS sur-faces. For a more detailed presentation on NURBS see [6] .A qth-degree NURBS curve is a piecewise rational function defined in parametric formasC(λ) =q∑i=0νiBiCi,q(λ)q∑i=0νiCi,q(λ)0 ≤ λ ≤ 1, (1)where {Bi} are the coordinates of the control points (forming the control polygon), {νi}are the control weights, and the {Ci,q(λ)} are the normalized B-spline basis functions of2R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huertadegree q, which are defined recursively as:Ci,0(λ) ={1 if λ ∈ [λi, λi+1[,0 elsewhere,Ci,k(λ) =λ− λiλi+k − λiCi,k−1(λ) +λi+k+1 − λλi+k+1 − λi+1Ci+1,k−1(λ),for k = 1 . . . q, where λi, for i = 0, . . . ,m, are the knots or breakpoints, which are assumedordered 0 ≤ λi ≤ λi+1 ≤ 1.A NURBS surface of degree q in λ direction and degree l in κ direction is a piecewiserational function defined in parametric form asS(λ, κ) =q∑i=0l∑j=0νijBij Ci,q(λ)Cj,l(κ)q∑i=0l∑j=0νij Ci,q(λ)Cj,l(κ)0 ≤ λ, κ ≤ 1, (2)where {Bij} are the coordinates of the control points (forming the control net), {νij}are the control weights, and {Ci,q(λ)} and {Ci,l(κ)} are the normalized B-spline basisfunctions of degree q and l respectively, defined on the so-called knot vectorsΛq = {0, . . . , 0︸ ︷︷ ︸q+1, λq+1, . . . , λmq−q−1, 1, . . . , 1︸ ︷︷ ︸q+1},Λl = {0, . . . , 0︸ ︷︷ ︸l+1, κl+1, . . . , κml−l−1, 1, . . . , 1︸ ︷︷ ︸l+1}.Note that each 2D B-spline basis function is a cartesian product of two 1D B-splinebasis functions. Figure 1 shows two 2D B-spline basis functions associated with knotvectorsΛq = {0, 0, 0, 0, 0.4, 1, 1, 1, 1},Λl = {0, 0, 0, 0.2, 0.6, 0.6, 1, 1, 1}.This figure also represents 1D B-spline basis functions associated to Λq and Λl knotvectors.Note that NURBS curves change their definition at breakpoints. Then, NURBS sur-faces change their definition at knot lines, that is, at (λi, κ) for i = 1, . . . ,mq and (λ, κj)for j = 1, . . . ,ml. An example of a NURBS surface is represented in Figure 2 with thecorresponding control net. Knot lines are represented on the NURBS surface in order tostress the discontinuous nature of the parametrization.3R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huerta00.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.8100.20.40.60.81Figure 1: Two examples of 2D B-splines basis functions. Complete 1D B-spline basis are also represented:in red for knot vector Λq and, in blue for knot vector Λl.Figure 2: Control net (left) and NURBS surface with knot lines (right)3 NURBS-ENHANCED FINITE ELEMENT METHOD (NEFEM)Consider a physical domain Ω ⊂ R3 whose boundary, ∂Ω, or a portion of it is definedby NURBS surfaces. Every NURBS is assumed to be parametrized byS : [0, 1]2 −→ S([0, 1]2) ⊆ ∂Ω ⊂ R3.A regular partition of the domain Ω =⋃eΩe in tetrahedrons is assumed such that everyelement Ωe has at most one face on the NURBS boundary. Figure 3 shows a domain withpart of the boundary described by a NURBS surface corresponding to a cylinder, and avalid tetrahedral mesh for NEFEM.As usual in FE mesh generation codes, it is assumed that every curved boundary facebelongs to a unique NURBS, ΥFe ⊆ S([0, 1]2). That is, one element face can not bedefined by portions of two (or more) different NURBS surfaces. But on the contrary, itis important to note that knot lines, which characterize the piecewise nature of NURBS,are independent of the mesh discretization. Thus, NURBS parametrization can change4R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huerta(a) (b)Figure 3: (a) Physical domain with part of the boundary defined by a NURBS surface and (b) NEFEMtetrahedral meshits definition inside one face, that is breakpoints may belong to element faces and do notneed to coincide with FE nodes. This is another major advantage with respect to theisogeometric analysis [3].Every element of a 3D NEFEM tetrahedral mesh has 0, 1 or 3 edges on the NURBSboundary. An element with no edges defined by NURBS has planar faces and is de-fined and treated as a standard FE or DG element. Therefore, in the vast majority ofthe domain, interpolation and numerical integration are standard. Specifical numericalstrategies for the interpolation and the numerical integration are needed only for thoseelements affected by the NURBS description of the boundary.An element with one edge on the NURBS boundary has two curved faces defined fromthe NURBS edge and an interior tetrahedral vertex, see Figure 4 (b). Finally, for anelement with three edges on the NURBS, all his faces are curved: a boundary face definedby NURBS is denoted by ΥFe , and a curved interior face defined from a NURBS edge isdenoted by ΥEe , see Figure 4 (a).3.1 FE polynomial basisFor each element Ωe, a tetrahedral Te is defined using its vertices. A linear mappingΨ : I −→ Te is used, which goes from the reference tetrahedral I to the tetrahedralTe. The inverse of this linear transformation maps the tetrahedral Te into the referencetetrahedral I and, more important, also maps the actual element Ωe, which is in thephysical domain, into a curved element in local coordinates, namelyIe := Ψ−1(Ωe).Ie is called the local curved element for the actual element Ωe.Note that the reference tetrahedral I is the same for all elements Ωe. However, thelocal curved element Ie depends on the trimmed NURBS defining the curved face ΥFe of5R. Sevilla, S. Ferna´ndez-Me´ndez and A. HuertaΥEeΥFeQP1P2(a)ΥEeΥEeP1P2Q1Q2(b)Figure 4: (a) An element with a face on the NURBS boundary and (b) an element with an edge onthe NURBS boundary. Curved faces on the NURBS boundary are represented in orange, curved facesdefined from an edge on the NURBS boundary in blue and planar faces in greenΩe, and therefore it is different for every element Ωe intersecting the NURBS boundary.In order to work with standard FE polynomial approximations, Lagrange polynomials(that is, standard nodal interpolation) can be considered. In fact, they can be defined onthe curved tetrahedral, Ie, in the reference domain or equivalently, in the actual elementin the physical domain, Ωe. The use of a linear transformation from the local (reference)coordinates ξ = (ξ, η)T in Ie to the cartesian coordinates x = (x, y)T in Ωe, ensures thata complete polynomial interpolation of degree m in ξ leads to a polynomial interpolationwith the same degree in x. Thus, consistency and accuracy of the approximation isensured even for elements Ωe far from being a straight-sided element.In order to make the computation of Lagrange polynomials, more systematic, for anydegree and for any distribution of nodes, the implementation proposed in [7] is adopted,see [1] for more details.3.2 Numerical integrationThe weak form that must be solved requires both integrations along element faces andin the element interiors. All integrals in elements not having an edge along the NURBSboundary are computed using standard procedures. Elements Ωe with one or three edgeson the NURBS boundary require special attention. Two cases must be considered for bothsurface integrals (usually related to the implementation of natural boundary conditions orto flux evaluation along the face in a DG context) or volume integrals (standard integralsin the element Ωe). As discussed in the previous section, since NEFEM uses polynomialsto approximate the solution, the difficulties in numerical integration are only restrictedto the definition of a proper numerical quadrature in a curved element, Ωe, and a curved6R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huertaface.3.2.1 Surface integralsTwo different integrals are studied separately: integrals on a curved face defined byNURBS, ΥFe , and integrals on a curved face defined by an edge on the NURBS boundary,ΥEe , see Figure 4 (a) and (b).A surface integral to be computed on a curved boundary face ΥFe defined by a trimmedNURBS, can be written as∫ΥFef dx dy dz =∫Λef(S(λ, κ)) ‖JS(λ, κ)‖ dλ dκ,where f is a generic function (usually polynomial), Λe is a triangle (or quadrilateral) inthe parametric space of the NURBS such that S(Λe) = ΥFe , and ‖JS(λ, κ)‖ denotes thenorm of the differential of the NURBS parametrization S (which, in general, is not apolynomial).Evaluation of the previous integral requires to take into account the piecewise natureof the NURBS definition, because several changes of NURBS parametrization can occurinside Λe. A delaunay triangulation of Λe is considered to define a numerical quadratureon ΥFe . The triangulation is constructed from Λe vertexs, intersections of knot lines withΛe and intersections of knot lines inside Λe. Then, different numerical quadratures areconsidered at each subelement of the triangulation. In this work symmetric quadraturesfor triangles are considered, see [8].On the other hand, surface integrals on interior curved faces must be computed whena DG formulation is adopted. An interior curved face ΥEe can be parametrized using thedefinition of the NURBS edge and the interior node Q of the tetrahedral element byΠQ : [0, 1]2 −→ ΥEe(%, σ) 7−→ (1− σ)S(λ(%), κ(%))+ σQ,where λ(%) = λ1+%(λ2−λ1) and κ(%) = κ1+%(κ2−κ1), with λ1, λ2 ∈ Λq and κ1, κ2 ∈ Λlsuch that S(λi, κi) = Pi, for i = 1, 2, i.e. P 1 and P 2 are the vertexs of the element onthe NURBS boundary, see Figure 4 (b).Then, the integral on the interior curved face can be written as∫ΥEef dx dy dz =∫ 10∫ 10f(ΠQ(%, σ)) ‖JΠQ(%, σ)‖ d% dσ,where f is a generic function (usually polynomial), and ‖JΠQ(r, s)‖ denotes the norm ofthe differential of the mapping ΠQ (which, in general, is not a polynomial). Again, theevaluation of the previous integral requires to take into account the piecewise nature ofthe NURBS parametrization. Several changes of NURBS definition can occur inside theboundary edge and composite quadratures must be considered for r direction.7R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huerta3.2.2 Element integralsNEFEM also requires to compute integrals in an element Ωe with two or four curvedfaces, see Figure 4. Thus, a different numerical quadrature for every curved element Ωeis needed.For an element with one face on the NURBS boundary, one alternative is to define atransformation from the prism Λe × [0, 1] to the curved element Ωe, namely,ϕ : Λe × [0, 1] −→ Ωe(λ, κ, ϑ) 7−→ ϕ(λ, κ, ϑ) = (1− ϑ)S(λ, κ) + ϑQ. (3)The numerical quadrature on Λe × [0, 1] is defined from a cartesian product of thequadrature on Λe proposed in Section 3.2.2 and a Gauss-Legendre quadrature on the ϑdirection. Then, the element integral is computed as∫Ωef dx dy dz =∫Λe∫ 10f(ϕ(λ, κ, ϑ)) |Jϕ(λ, κ, ϑ)| dλ dκ dϑ,On the other hand, an element with a curved face not on the boundary, can beparametrized using the definition of the edge on the NURBS boundary and interior nodes,Q1 and Q2, of the element Ωeφ : [0, 1]3 −→ Ωe(%, σ, τ) 7−→ φ(%, σ, τ) = (1− τ)ΠQ1(%, σ) + τQ2.(4)Then, integral is evaluated using Gauss-Legendre quadratures in each direction. To ac-count for the rational definition of the NURBS edge composite quadratures are consideredin the parameter %. Then, the integral is computed as,∫Ωef dx dy dz =∫ 10∫ 10∫ 10f(φ(%, σ, τ)) |Jφ(%, σ, τ)| d% dσ dτ,Note that parametrizations (3) and (4) are linear in the third parameter (ϑ and τrespectively). Then, exact integration can be carried out in this direction with p + 1integration points.4 NUMERICAL EXAMPLEThe following model problem is solved in three dimensions−∆u = f in Ωu = ud on Γd∇u · n = gn on Γn8R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huertawhere Ω is an sphere of unit radius, Γd ∪ Γn = ∂Ω and n is the outward unit normalvector on ∂Ω. The source is given by f(x, y, z) = x cos(y) + y sin(z) + z cos(x), in such away that the analytical solution of the problem is known and smooth,u(x, y, z) = x cos(y) + y sin(z) + z cos(x).Dirichlet boundary conditions, corresponding to the analytical solution, and Neumannboundary conditions, also corresponding to the analytical normal flux, are imposed onΓd = ∂Ω ∩ {(x, y, z) ∈ R3 | z ≤ 0} and Γn = ∂Ω ∩ {(x, y, z) ∈ R3 | z > 0} respectively.FEM and NEFEM solutions in a mesh with only eight elements are represented inFigure 5 for quadratic and cubic approximation. It is important to note that not only thesolution is captured with lower accuracy with isoparametric FE but also geometric errorsare observable.(a) FEM p=2 (b) NEFEM p=2(c) FEM p=3 (d) NEFEM p=3Figure 5: Surface plot of FEM and NEFEM solutions5 CONCLUDING REMARKSThe 3D extension of NEFEM is studied in this work. Although the extension is concep-tually easy, special attention must be paid to geometrical aspects. Specifically designed9R. Sevilla, S. Ferna´ndez-Me´ndez and A. Huertapolynomial interpolation and numerical integration is needed for elements affected by theNURBS description of the boundary. Several strategies are presented to perform thenumerical integration on curved elements.Finally, a Poisson example is solved to demonstrate the applicability of the proposedmethod in a continuous Galerkin framework.REFERENCES[1] R. Sevilla, S. Ferna´ndez-Me´ndez, and A. Huerta, “NURBS-Enhanced Finite ElementMethod (NEFEM),” 2007. Submitted.[2] F. Bassi and S. Rebay, “High-order accurate discontinuous finite element solution ofthe 2D Euler equations,” J. Comput. Phys., vol. 138, no. 2, pp. 251–285, 1997.[3] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finiteelements, NURBS, exact geometry and mesh refinement,” Comput. Methods Appl.Mech. Eng., vol. 194, no. 39–41, pp. 4135–4195, 2005.[4] B. Cockburn, Discontinuous Galerkin methods for Computational Fluid Dynamics,vol. 3 (Fluids) of Encyclopedia of Computational Mechanics, ch. 4. John Wiley &Sons, 2004.[5] D. Xue and L. Demkowicz, “Control of geometry induced error in hp finite element(FE) simulations. I. Evaluation of FE error for curvilinear geometries,” Int. J. Numer.Anal. Model., vol. 2, no. 3, pp. 283–300, 2005.[6] L. Piegl and W. Tiller, The NURBS Book. London: Springer-Verlag, 1995.[7] J.S. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured grids.I. Time-domain solution of Maxwell’s equations,” J. Comput. Phys., vol. 181, no. 1,pp. 186–221, 2002.[8] S. Wandzura and H. Xiao, “Symmetric quadrature rules on a triangle,” Comput. Math.Appl., vol. 45, no. 12, pp. 1829–1840, 2003.Con el apoyo de Universitat Polite`cnica de Catalunya y E.T.S. d’Enginyers de Camins, Canals i Portsde Barcelona10

3D nurbs-enhanced finite element method

http://hispana.mcu.es/es/registros/registro.do?tipoRegistro=MTD&idBib=1094775

3D nurbs-enhanced finite element method

Abstract

Similar works

Full text

Available Versions

UPCommons. Portal del coneixement obert de la UPC