High performance adaptive finite element modeling of complex CAD geometry

(1)

High performance adaptive finite

element modeling of complex CAD

geometry

S T E F A N I E S T R U N K

Master of Science Thesis Stockholm, Sweden 2013

(2)

(3)

High performance adaptive finite element

modeling of complex CAD geometry

S T E F A N I E S T R U N K

Master’s Thesis in Scientific Computing (30 ECTS credits) Master Programme in Scientific Computing 120 credits Royal Institute of Technology year 2013 Supervisor at KTH ware Johan Hoffman and S M Ashraful Kadir

Examiner was Michael Hanke TRITA-MAT-E 2013:40 ISRN-KTH/MAT/E--13/40--SE

Royal Institute of Technology School of Engineering Sciences

KTH SCI SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

(4)

(5)

Abstract

CAD (Computer Aided Design) and finite element analysis are of fundamental importance for numerical simulations. The general approach is to design a model using CAD software, create a mesh of a domain that includes this model and use finite element analysis to perform simulations on that mesh. When using more advanced simulation techniques, like adaptive finite element methods, it is more and more desired to use CAD information, not only for the creation of the initial mesh but also during the simulation.

In this thesis, an approach is presented how to use CAD data during adaptive mesh refinement in a finite element simulation. An error indicator is presented to find the elements in a mesh, which need to be improved for a better geometric approximation and it is shown how to integrate the different approaches into an existing high performance finite element solver.

(6)

(7)

Referat

Adaptiv finita-element-modellering av complex CAD-geometri CAD (Computer Aided Design) och finita-element-analys är grundläggande för numerisk simulering. Man konstruerar en modell med CAD-program, skapar ett beräkningsnät påen domän som innehåller modellen, och använder finita-element- analys för beräkningar på nätet. I mer avancerade simuleringar, som för adaptiva finita-element-metoder, är det önskvärt att använda CAD information inte bara för att skapa det första nätet utan under nätförfiningarna i adaptionen under simu- leringen. I detta arbete presenteras ett sätt att använda CAD-data för adaptiv nätförfining i en finita-element-simulering. En fel-indikator ges för att hitta de element som ska förfinas för att förbättra geometrisk approximation och vi beskriver hur de olika angreppssätten kan integreras i ett finita-element programpaket för högpresterande datorer.

(8)

(9)

Introduction

Finite element methods are a numerical tool for approximating the solution of differential equations. The general idea is to divide the continuous domain into a finite number of elements. On these finite elements, one defines basis functions and, after inserting these basis functions into the differential equation, one obtains a system of equations, which can be solved numerically.

To obtain a sufficiently accurate solution, this method requires a fine enough mesh. But using a mesh that is very fine in all regions of the domain is computa- tionally expensive. With the help of a posteriori error indicators, that are obtained by solving the primal and an associated dual problem [1] on the coarse mesh, it is possible to start with the minimum amount of initial elements and refine the mesh in regions that require a higher resolution. This optimized method is called

“adaptive finite element method” and has the great benefit that it minimizes the computational cost.

The meshes for such simulations are usually created using CAD (Computer Aided Design) models. However, when using a coarse initial mesh the approximation of the geometries of such models is also very coarse. Straightforward refinement techniques are not able to refine the mesh with respect to the geometry of the original CAD model. This means, also the adaptively refined meshes do not provide a better approximation of the geometry than the initial coarse mesh. Therefore, an integration of the geometric information into the refinement process is desirable.

In this thesis, an approach is presented to gradually refine the mesh according to the geometry of the CAD model. For that, three components were implemented and tested: (i) an initial projection to establish a connection between the initial mesh and the geometric model, (ii) a geometric error estimation, which allows to measure the deviation of the mesh from the CAD model and (iii) an improved refinement algorithm that creates new elements, such that the approximation of the geometric boundary is improved.

An example for such a refinement procedure is given in figure 1.1. It shows the two dimensional model of three circular arcs, which are initially approximated by an extremely coarse mesh. After multiple local mesh refinements on the boundary,

1

(12)

2 CHAPTER 1. INTRODUCTION

Figure 1.1. The idea of local mesh refinement on the boundary with a geometric correction.

one can see that the mesh in the lower right corner of figure 1.1 approximates the geometry almost exactly without the need to have a fine mesh resolution in all regions.

To obtain such a local refinement with respect to geometric data, two possible procedures are available: One can either do refinement with respect to geometrical boundaries, i. e., add elements at domain boundaries that give a continuously improved approximation of the geometry, or one can use advanced element types to build up the mesh. For example, the shape of such elements can be described with NURBS, as it was done by R. Sevilla and Huerta [2]. But such an approach requires a different finite element method that is able to account for such complex element shapes. Among others, this approach was introduced by R. Sevilla and Huerta [2]

or Hughes et al. [3]. An integration of NURBS based elements into a finite element solver for high performance computing was discussed by Cottrell et al. [4].

But since the goal of this project was to extend an already existing high performance finite element framework, which provides assembly and solution routines based on triangular or tetrahedral meshes, this was not a suitable approach.

Therefore, we decided to improve the geometric representation with each refinement step, but otherwise keep the finite element analysis unchanged. Ashraful Kadir introduced in his master thesis [5] a serial approach, which consists of funda-

(13)

3 mental algorithms from computer graphics. The most important algorithm is the point projection of a point in space onto a curve or surface. A detailed discussion of different point projection algorithms can be found in [5].

Since adaptive finite element methods rely on a posteriori error indicators to mark the necessary refinement areas, such an indicator for regions, where the mesh does not represent the geometry sufficiently accurate, is necessary. Much research has been done in using curvature informations to generate a geometric error estimate. Frey and Borouchaki [6] for example used an error analysis based on the Hessian of the surface or Kim et. al. [7] introduced an approach using the discrete curvature norm. In contrast, the approach presented in this thesis measures the area or volume (for 2D- or 3D-meshes, respectively) enclosed between mesh boundary and geometry. Naturally, a larger area or volume corresponds to a larger deviation from the geometry and hence, the geometric error is larger. A straightforward approach was used to approximate the area or volume between the mesh and the geometry, with the help of few sampling points.

The next chapter will introduce some basic concepts, to give the necessary background and fundamental mathematical methods for the used algorithms. In chapter 3 a short overview over the FEniCS project and the library DOLFIN HPC will be given, which was used for the implementations in this project. The geometric correction algorithms itself are part of chapter 4 and their optimization and outcomes are discussed in chapter 5. A summary and an outlook onto possible follow up developments are given in chapter 6.

(14)

(15)

Chapter 2

Basics Concepts

Geometries, like the shape of a car or an aircraft, can be represented in various ways. A widely used tool are NURBS (Non Uniform Rational B-Splines), which represent shapes in terms of weighted polynomial curves or surfaces. This way, the representation is both simple and powerful and hence, this was the method of choice. A short introduction into NURBS can be found in appendix A and the relevant literature, e. g., Piegl and Tiller [8] or Dimas and Briassoulis [9].

Three different components are considered throughout this thesis:

initial projection a preprocessing step that moves boundary vertices exactly onto the geometrical surfaces or curves and determines corresponding geometric parameters

geometric error a quantity to measure the deviation of the mesh from the geometry

extended refinement add new cells in locally marked regions and improve this way the approximation of the geometry

To obtain a basic understanding of a geometry and how it is represented, this first section gives a short introduction into geometrical modeling. The central mathematical method, used by all three components, is point projection and explained in the second part of this chapter.

2.1 Representation of geometries

A complex geometry consists of several geometric entities. The most important ones are curves and surfaces, which are in the following denoted as C(u) or S(u, v), respectively. How they are defined in the IGES standard, which was used in this project, is described in the specifications provided by the U.S. Product Data Asso- ciation [10]. A short introduction to the most important aspects of this standard is given in appendix B.

5

(16)

6 CHAPTER 2. BASICS CONCEPTS Simple geometries, like a circle or cylinder, can be represented by few or even just one geometric entity. However, relevant simulation models are a lot more complex than that. Often, up to hundreds or thousands of geometric entities are used to create models like cars or airplanes.

2.2 Point Projection and Point Inversion

Point inversion and point projection are both elementary concepts in geometric modeling. Point inversion determines the corresponding geometric parameters for a point already located on a geometry, i. e. u for a curve and (u, v) for a surface.

Point projection is more general applicable and allows to find the point on the geometry that is closest to an arbitrary point in space. Since this is what is needed throughout all previously mentioned components, this method is explained in the following.

To find the shortest distance between a point in space and a curve one needs to orthogonally project this point onto the geometry. Figure 2.1 shows the projection of the point P onto the curve C(u). C⁰(u) denotes the first derivative of the curve.

The point C(ui) is the point on the surface with the shortest distance to the point P.

b P

C^′(ui) C(ui)

b

Figure 2.1. Projection of the point P onto the curve C(u)

2.2.1 Formulation as a minimization problem

The point projection onto a NURBS curve can be formulated as a minimization problem, as it was done by Piegl and Tiller [8]:

Definition 2.1. The distance d between a given point P = (x, y, z) and a curve C(u) is minimum if f(u) = 0 with

d(P, C(u)) = inf{kC(ui) − P k |C(ui) ∈ C(u)} (2.1) and

f(u) = C⁰(u) · (C(u) − P ). (2.2)

(17)

2.2. POINT PROJECTION AND POINT INVERSION 7 In other words: If a point on the curve (C(ui)) is found, where the vector between the point C(ui) and P is orthogonal to the tangential at the curve point, one found the orthogonal projection of P onto C(u), hence the shortest distance.

Whithin this thesis, we assume that there exits always at least one point, such that the previous definition holds true.

2.2.2 Iterative Newton method and initial screening

In this project the approach suggested by Piegl and Tiller [8] was used to solve the minimization problem given in eq 2.2 using the iterative Newton method. The idea is to start from a sufficiently good initial guess for the sought-after geometric value and successively improve it, until a desired accuracy is reached.

One iteration of the Newton method that gives the improved value for the geometrical parameter can be formulated as

u_i+1= ui− f(ui)

f⁰(ui) = ui− C⁰(ui) · (C(ui) − P )

C⁰⁰(ui) · (C(ui) − P ) + |C⁰(ui)|². (2.3) In order to achieve a reliable convergence for the iterative scheme, one needs a suitable start value u0. To find this value, an initial screening is performed. A number of equidistant geometric parameters are tested for the distance of their corresponding points to P , as shown in figure 2.2. The parameter with smallest distance is consequently chosen as initial guess, for the example in the picture this is u3.

With this initial value it is possible to start above Newton iteration until the geometric parameter, associated with the point closest to P on the curve, is found.

P

b

bb b b b b b b

C(u3)

b

Figure 2.2. Initial screening to find a suitable start value

Convergence criteria

After each iteration step (equation 2.3) one has to check for convergence. The first check is the Euclidean distance between P and C(ui):

|C(ui) − P | ≤ ε1 (2.4)

Then one tests the cosine similarity (equation 2.5), which means that one calculates the cosine of the angle between the derivative C⁰(ui) and the vector between C(ui)

(18)

8 CHAPTER 2. BASICS CONCEPTS and P. The value is bound by -1 and 1, whereas a cosine of 0 means that the vectors are orthogonal to each other.

|C⁰(ui) · (C(ui) − P )|

|C⁰(ui)||C(ui) − P | ≤ ε₂ (2.5) If one of the equations 2.4 or 2.5 is satisfied, the iteration is stopped. Otherwise, a new value ui+1 is calculated with the iteration formula 2.3, for which the following has to be fulfilled:

First of all, one has to make sure that the value is still in the range, i. e., u_i+1∈[umin, u_max]. Special care is needed if the curve is closed, i. e., if origin and endpoint of the curve are the same point. This means, this point has two associated geometrical parameters. First the case of an open curve: If the newly calculated geometric parameter exceeds the bounds, one just has to set the value of ui+1 to u_min or umax, respectively. If the curve is closed and ui+1 is smaller than umin we can set ui+1= umax−(umin− u_i+1). Similar in the opposite case: if ui+1 is bigger than umax, one has to adjust it in the following way: ui+1= umin+ (ui+1− u_max).

This means, whenever the bounds are violated the difference (umax− u_min) is added or subtracted to bring the parameter back into the interval.

Additionally, after each iteration is tested whether the parameter did change significantly:

|(ui+1− u_i)C⁰(ui)| ≤ ε1 (2.6) If the variation is smaller than the threshold ε1, the geometric parameter is considered to be converged and the iteration is stopped.

(19)

Chapter 3

The FEniCS project

The implementation of this boundary smoothing method was done using the finite element framework FEniCS. The FEniCS project [11] consists of several open source components and has the aim to efficiently solve differential equations and automate the process. This project is developed in cooperation between several different research groups.

The core component of FEniCS is DOLFIN [12], a C++ library that provides functionalities for mesh representation, FEM assembly, interfaces to linear algebra backends. Moreover, it is the main user interface.

A short introduction to DOLFIN will be given in the following section. Fur- ther information can be found in Logg et al. [11], which provides a fundamental background to the FEniCS project, including the mathematical background, and information on the software design.

3.1 DOLFIN and DOLFIN HPC

As already mentioned, DOLFIN is the core component of the FEniCS project and serves as the main user interface. It combines functionalities of several FEniCS components and other external packages and manages the communication between those components. Furthermore, it provides the data structures for simplex meshes as well as necessary algorithms for mesh refinement.

As the tasks to simulate get more and more complicated it is simply not possible anymore to simulate everything serially. To be able to use modern supercomputers, DOLFIN HPC [13] was developed as an optimized version of DOLFIN with respect to parallel computing. DOLFIN HPC is parallelized using a hybrid MPI/OpenMP approach.

A short explanation of the mesh data structures and the two different refinement algorithms used in DOLFIN HPC will be given, since the geometric correction will be directly incorporated into these algorithms.

9

(20)

10 CHAPTER 3. THE FENICS PROJECT

3.1.1 Meshes

A mesh in DOLFIN HPC is defined in terms of vertices and cells. All remaining entities, like edges, are derived from those. In parallel computations, each process stores a subset of the global mesh. All entities are owned exclusively by one process, but to be able to build up cells using vertices that are owned by a neighboring process, it is necessary to keep copies of those vertices as ghost entities. An additional object manages the global connectivity and assigns each local mesh entity a global index.

In DOLFIN HPC two refinement algorithms are available: uniform refinement, using a simple edge bisection scheme, and recursive longest edge bisection.

Simple edge bisection

The simple edge bisection scheme performs a uniform refinement of the whole mesh.

A new vertex is created at the midpoint of each edge and new cells are created by connecting all new and old vertices in each cell. This way, a triangle (2D) or tetrahedron (3D) is divided into 4 or 8 elements, respectively. Figure 3.1 shows a 2D-example.

Figure 3.1. Uniform refinement with simple edge bisection. Each cell of the initial (black) mesh is subdivided into four new cells (blue) by bisection of all edges.

Recursive longest edge bisection

To exploit the advantage of an adaptive finite element method, one needs a local mesh refinement. In DOLFIN HPC, a recursive longest edge bisection algorithm is implemented, which was first introduced by Rivara [14]. In [13] one can find details how this algorithm is implemented in DOLFIN HPC.

(21)

3.2. UNICORN HPC 11 The basic idea is to recursively bisect the longest edge until no hanging vertices remain. The algorithm, as it was developed by Jansson [15], is outlined in algorithm 3.1.

Algorithm 3.1 Parallel recursive longest edge bisection

Let R be the set of elements marked for refinement, T be the initial mesh, parti- tioned into k subsets Tⁱ

1: Bisect all elements from R ∩ Tⁱ

2: Propagate hanging nodes to neighboring processing elements

3: Mark elements with hanging nodes for bisection

4: if R= ∅ then

5: return

6: else

7: goto 1.

8: end if

Figure 3.2 shows an example how this refinement is done. In a first step all cells that need refinement are marked (gray cells in the picture). Then the longest edge of each marked cell (red) is determined and bisected. After that, it is possible that hanging vertices remain. In the picture, one of the bisected edges is on the boundary, hence the new vertex is not a hanging vertex. The other bisected edge is in the interior and the newly introduced vertex is now a hanging vertex (marked with the green rectangle). To resolve this, one needs to search for the opposite cell and recursively repeat the refinement, until all hanging vertices are resolved.

This recursive bisection is done locally on each process. But hanging vertices can also remain on process boundaries, hence these have to be propagated to the neighboring processing elements (PE). As soon as all processes are idle, global termination is achieved. To reduce the execution time of the refinement algorithm, load balancing is applied before to equally distribute the marked cells among all PEs.

3.2 Unicorn HPC

Unicorn HPC [16] is a unified continuum mechanics solver working on top of DOLFIN HPC. It uses the G2 (General Galerkin) method [1], a variant of FEM to discretize the Navier-Stokes equations, and uses an adaptive algorithm to obtain the desired accuracy. A short introduction into this method is given in the context of the flow simulation in chapter 5.

(22)

12 CHAPTER 3. THE FENICS PROJECT

Figure 3.2. Recursive longest edge bisection: two cells are marked to be refined. In the left upper corner is the mesh with the marked cells. The right upper corner shows the longest edge of each marked cell. The left lower image is after the refinement of the marked cells and before any hanging nodes are resolved. The right lower figure is the final refined mesh. The curve represents a part of the geometric boundary, which is supposed to be approximated by the mesh.

(23)

Chapter 4

Geometric correction algorithms and

their integration into FEniCS

After the introduction of DOLFIN HPC and Unicorn in the previous section, this chapter explains the three components of geometric correction and how they were incorporated into this existing library. This required to introduce a number of additional libraries and changes to DOLFIN HPC itself.

4.1 Additional Libraries

Two new libraries were necessary for the implementation of geometric boundary smoothing in DOLFIN HPC: NURBS++ and LIBGEOM. These will be shortly presented in the upcoming sections.

4.1.1 NURBS++

NURBS++ is an OpenSource library that provides data structures to represent surfaces and curves and basic algorithms to manipulate those. It was developed by Lavoie [17]. A detailed description of the library can be found in the official documentation [18].

4.1.2 LIBGEOM

The geometry library LIBGEOM connects NURBS++ and DOLFIN HPC. It was first developed by Ashraful Kadir [5] and extended during this project.

This library allows to parse CAD data, in particular geometries as defined in the IGES format (see appendix B ), store the different entities and provide both an interface to algorithms implemented in NURBS++ and own methods, like the geometric error calculation.

13

(24)

14 CHAPTER 4. GEOMETRIC CORRECTION ALGORITHMS

4.1.3 Connection

Figure 4.1 shows how the different modules are connected with DOLFIN HPC and Unicorn. LIBGEOM is the main user interface for geometric operations and hence used by both DOLFIN HPC and applications like Unicorn directly.

Figure 4.1. Connection between the different modules for boundary smoothing

4.2 Boundary smoothing

Boundary smoothing makes use of the three components that were briefly introduced in the beginning of chapter 2: (i) initial projection of boundary vertices onto the geometry, (ii) determination of the geometric error and (iii) mesh refinement with respect to geometrical boundaries. Each of these three components will be explained in the following.

4.2.1 Initial projection

Usually the coarse mesh will be created with the help of special CAD software. It is possible that the boundary vertices of the mesh don’t represent the geometric boundary exactly. Since in our case it is important to have a good initial pre- sentation of the geometric boundary, we can move the mesh boundary points, in a preprocessor step, exactly onto the geometry. A bounding box is chosen, that defines a mesh region in which this initial projection should be performed. This box can also have the size of the entire domain. All boundary vertices of the mesh, which are inside this bounding box, are projected onto the geometry using the point projection algorithm explained in section 2.2 and are moved to this projected point.

The outline of the initial projection algorithm can be found in algorithm 4.1.

Algorithm 4.1 Initial projection - preprocessing step

Let B be the boundary of the mesh M and S a subset of M with all the vertices included in the bounding box.

for all v ∈ B ∩ S do

find closest point on the geometry

store the corresponding parameter u, v and the patch-id save the new vertex position in the mesh

end for

(25)

4.2. BOUNDARY SMOOTHING 15 After performing this preprocessing step, geometric parameters are known for all boundary vertices and can be used consequently during the mesh refinement and geometric error computation.

4.2.2 Geometric error

The geometric error is a quantity that measures the deviation of the mesh from the geometric boundary. To estimate a geometric error we approximate the area (2D) or volume (3D) between the geometry and the mesh. The size of the area or volume is then the quantity for the error.

As before, for the sake of simplicity this explanation will be restricted to the 2D case. The 3D case is analog and its explanation can be found in appendix C.2.

The idea is to use the midpoint of a geometric boundary edge and determine the according projection on the geometry. With the two known edge vertices and the projected point, one can calculate the area of the triangle spanned by those three.

Applying this step recursively, one can calculate the area more precisely.

The algorithm is outlined in algoritm 4.2. Let e be an edge of the mesh, with its vertices v1, v2 ∈ G where G is the set of vertices, which are located on the geometric boundary. Figure 4.2 shows a two dimensional example with 2 recursive steps.

Algorithm 4.2 Geometric error calculation in 2D function recursiveError(e, recursion steps)

if recursion steps = 0 then return 0

end if

calculate midpoint

project midpoint onto the geometry

error = area(triangle(v1, v2, projected point)) e1 = edge(v1, projected point)

e₂ = edge(projected point, v2)

. Recursive calls improve the error estimate error += recursiveError(e1, recursion steps - 1)

error += recursiveError(e2, recursion steps - 1) return error

end function

4.2.3 Improvement of the boundary approximation using mesh refinement

The previous sections explained methods that allowed to position boundary vertices more exactly on the geometric boundary and to measure the deviation of the mesh from the geometric shape. But the main goal is to reduce this geometric error, i. e., improve the approximation of the geometric boundary through the mesh.

(26)

16 CHAPTER 4. GEOMETRIC CORRECTION ALGORITHMS

Figure 4.2. Geometrical error: The ellipsoid is initially represented by 5 vertices (red) forming 4 cells. In the first step the green triangle areas are calculated. In the next step two blue triangle areas are added to each green area. The triangles are created by bisecting the boundary edge and projecting the midpoint onto the geometry.

DOLFIN HPC provided already mesh refinement algorithms (see section 3.1.1), hence, it was natural to extend those such that they include boundary smoothing techniques. The basic idea is very simple: whenever a boundary edge is bisected, the newly introduced vertex is positioned exactly on the geometric boundary instead of the midpoint of the edge. In this sense, it does the same as the initial projection, which also moved boundary vertices onto the geometric boundary. However, in this case it is done immediately when the vertex is created.

In the following, the set G denotes all vertices which were projected in the preprocessing step. Whenever both edge endpoints of the edge, which needs to be bisected, are contained in G, the newly created point has to be projected and moved onto the geometric boundary. Algorithm 4.3 outlines the modification of the refinement algorithm.

Algorithm 4.3 Refinement with boundary smoothing determine the edge which needs to be bisected Let v1, v2 be the vertices of the edge

Let vm be the midpoint between v1 and v2

if v₁∈ G and v2∈ G then

vm = projection of the midpoint on the geometry store geometric parameters of the new point end if

add vm to the mesh

Figures 4.3 and 4.4 illustrate again how a geometric correction would impove the geometric mesh approximation. On the left sides of figures 4.3 and 4.4 one can see the result of the original refinement. In comparison on the right sides, one can see

(27)

4.2. BOUNDARY SMOOTHING 17 that boundary smoothing with respect to the geometry did improve the geometric approximation.

Figure 4.3. Simple edge bisection with boundary smoothing (right) and without (left)

Figure 4.4. Recursive longest edge bisection: on the left without geometric correc- tion and on the right with geometric correction

4.2.4 Parallelization approaches

The described algorithms are serial, but several approaches are imaginable to adapt these to distributed memory architectures.

The point projection algorithm introduced in section 2.2 is designed to search the shortest distance to one entity. But, as explained in section 2.1, geometric models often consist of several geometric entities. This means, in case of a complex

(28)

18 CHAPTER 4. GEOMETRIC CORRECTION ALGORITHMS geometry one has to search for each vertex through all entities to find the correct curve/surface.

Therefore, if a large number of vertices has to be projected, this can be extremely expensive, since for each vertex all geometric entities have to be considered. Hence, an even distribution of this work among all processes is crucial. More details about this load balancing step are given in the next chapter.

An approach to speedup the search for the correct geometric entity would be to distribute the geometry by assigning a subset of entities to each process. This would allow all processes to search for the correct entity at the same time, consequently decrease the time to find the projection of one point. Unfortunately, this approach has an enormous communication overhead, since each point has to be sent to all processes and the data needs to be recollected after the search.

A different approach would be to store the entire geometry on every process and use the already distributed mesh in DOLFIN HPC. This means, the projection of a point is executed locally within one process, but multiple points can be projected at the same time. Since the existing refinement methods were already fully parallel with respect to mesh entities, it was natural to choose this approach and let every process store the entire geometry instead.

4.3 Requirements for the initial mesh

In this section the importance of the quality of the initial mesh is discussed. Maybe the most important aspect is that we require at least a certain degree of correct approximation of the geometric boundary in the initial mesh. This means that the geometric boundary vertices should lie at least close to the geometry. An example was constructed to show two problems that can occur if the vertices are not close to the geometry. This was especially constructed to show the problems, no useful results were expected. Figure 4.5 shows on the left side the coarse mesh before the projection. The boundary mesh is colored in red and the geometry is outlined in green. The triangles, which will cause problems after the projection, are marked in yellow and orange and pink.

In the orange case, one point is already close to the geometry and due to the coarse boundary approximation, the two other points of this triangle are also supposed to be on the boundary. After the projection of the two vertices the triangle will collapse into a line. The second problem happens by projecting the boundary vertices around the yellow triangle. As you can see in figure 4.6 the cell will be inverted. This will result in an invalid mesh.

A third problem is: if the initial mesh already has cells with a poor mesh quality, we cannot expect an improvement. The pink cells in figure 4.5 illustrate this problem. This third problem could be avoided by applying a suitable mesh smoother after the initial projection.

(29)

4.3. REQUIREMENTS FOR THE INITIAL MESH 19

Figure 4.5. On the left the initial mesh of the airfoil. the boundary mesh is in red and the geometry is marked in green. The yellow, orange and pink cells demonstrate different problems during the initial projection. On the right is the projected mesh

Figure 4.6. Close up of the leading edge before the projection and after the projec- tion. The problematic cells are marked in yellow and red.

(30)

(31)

Chapter 5

Optimizations and Results

After the explanation of the implemented algorithms in the previous sections, this chapter describes the tested optimizations to these methods and presents them along with the respective results.

5.1 Used models

For the performance analysis mainly two models were used: The front of an airplane and a wing section built from a NACA 0012 airfoil.

5.1.1 Front of an Airplane

The complete airplane (figure 5.1) geometry consists of 962 surface patches. For the tests, only the front of the airplane was used, which are only 52 surfaces. The corresponding initial airplane mesh has 694236 vertices. To limit the region to the cock- pit, a bounding box was created (x ∈ [0.05; 0.2], y ∈ [−0.25; 0.25], z ∈ [−0.2; 0.2]).

This bounding box contained 1278 boundary vertices. This restriction to the front was done because the rest of the airplane was already very good approximated by the initial mesh and the effects of boundary smoothing would not be well visible.

5.1.2 NACA Airfoil

The second model is a wing section built from a NACA 0012 airfoil [19]. The initial mesh was an extremely coarse approximation of this model and a finer approximation was created with three uniform refinement steps. Figure 5.2 shows the area around the airfoil in the initial mesh as well as after the first two refinement steps.

In Figure 5.3 and figure 5.4 one can compare pictures of the complete initial mesh and the mesh after three uniform refinements.

For the performance test we used the three times uniformly refined mesh. The mesh has 21904 vertices out of which 416 are on the geometric boundary around the airfoil. The geometry itself consists of 38 surfaces.

21

(32)

22 CHAPTER 5. OPTIMIZATIONS AND RESULTS

Figure 5.1. The rendered image of the airplane (top). Below one can see a close up of the front

Figure 5.2. Refine the mesh around a NACA wing section. Starting with an extremely coarse mesh, performing three uniform refinement steps. In this figure you can see the initial coarse mesh left, in the middle after one refinement and on the right after the second uniform refinement step.

5.2 Performance optimizations

All performance tests in this section were performed on a single machine with an Intel Core i7-3770 CPU (3.40 GHz, 4 cores/8 threads) and 16 GBytes of RAM. The measured times are given as wall clock times and the results are in seconds.

(33)

5.2. PERFORMANCE OPTIMIZATIONS 23

Figure 5.3. Initial coarse mesh.

Figure 5.4. Mesh after 3 uniform refinement steps with geometric correction

5.2.1 Load balancing in the initial projection

As briefly mentioned when explaining the possible parallelization approaches, an even distribution of the workload is important. This problem arises in the initial projection of the boundary vertices since the mesh is equally distributed over the processes without any consideration of the geometric boundaries. This means, some processes might only have interior vertices while others have a large amount of vertices on the boundary. When executing the initial projection with this partitioning of the mesh, processes without any boundary vertices would remain completely idle

(34)

24 CHAPTER 5. OPTIMIZATIONS AND RESULTS while others have to project a large amount of points. Hence, a proper redistribution prior to the execution of the initial projection is crucial.

An example for such a situation is given in figure 5.5, which shows the original distribution of the mesh of a rectangular domain that contains a circle. In the original distribution all the geometric boundary points belong to process one, since a projection onto the model is only necessary for points at the circle. After applying load balancing, the second process owns the same amount of geometric boundary vertices, as can be seen in figure 5.6.

This load balancing was performed using the already implemented load balancing methods of DOLFIN HPC, which are explained in Jansson [20]. The load balancer redistributes the mesh according to weights assigned to each cell, hence, assigning a larger weight to all cells that contain boundary vertices gave the desired mesh distribution.

Figure 5.5. The mesh of a circle in a box distributed over 2 processes. All triangles in red belong to process one and respectively all in green belong to process two.

Figure 5.6. Without load balancing (left) the geometric boundary belongs com- pletely to process one. After load balancing (right) both processes own the same amount of geometric boundary points.

Performance

To measure the runtime improvement when using load balancing, the airplane front was used. As already mentioned, this mesh has 1278 vertices in the bounding box to be projected. In addition to the load balancing, also an initial screening was applied with 8 sample points on each patch in each direction, as it will be explained in the next section. Figure 5.7 shows the measured run times of the preprocessing step with and without load balancing on different numbers of processes. It is obvious that

(35)

5.2. PERFORMANCE OPTIMIZATIONS 25 in the initial distribution of the mesh one process owned all the marked geometric boundary vertices. The distribution of the workload onto multiple processes reduced the run time significantly.

1 2 3 4

400 600 800 1,000 1,200

Processes

timeinsec

with load balancing without load balancing

Figure 5.7. Initial Projection of the airplane front with and without load balancing.

5.2.2 Initial screening to determine the geometric entity

A complex geometry consists of several surfaces and curves. Hence, it would be extremely expensive to perform the point projection on each geometric entity. In- stead, an initial search has to be performed to determine the closest entity to the point first, before the actual point projection can be executed.

In section 2.2 the method to determine a suitable start value for the Newton iteration was explained. The same approach can also be used to determine the closest entity. A small number of sampling points is used on each entity, to determine the entity with the closest distance, before a more exact search on this entity is performed.

A two dimensional example with three curves can be seen in figure 5.8. In this example 4 sampling points per curve are used.

Performance

The times in figure 5.9 were produced by initially projecting all boundary points of the airplane front onto the geometry. Both set ups, with and without initial screening, ran on 4 processes.

5.2.3 Use a-priori knowledge about entities during refinement

When performing the point projection during the refinement of boundary cells, a-priori knowledge about suitable geometric entities can be used. Often, both end-

(36)

b

Figure 5.8. Geometry with three curve entities. Do determine the closest entity one performs an initial screening with only a few sampling points on each entity. Here the black curve is the closest.

sec

1 10 100 1000 10000

100000 ≈54000

without initial

screening

432.107

with initial

screening

Figure 5.9. Initial screening of the airplane front on 4 processes with and without initial screening.

(37)

5.3. GEOMETRIC ERROR 27 points of the bisected edge are located on the same geometric entity, consequently also the new vertex that is created between them will belong to this entity. Hence, the search for the correct entity is omitted and the point projection can directly be executed onto this geometric entity. Only in cases where the two endpoints of the edge are associated with different entities, one has to search for the suitable patch.

Figure 5.10 shows the time benefit if one uses this optimization. In the set up the NACA wing section was refined on 4 nodes. The cells, which were marked for refinement, were determined using the geometric error indicator, as explained in the next section. Furthermore, one can see in figure 5.16 how many points were projected in each refinement step.

It might sound promising to limit the search to the geometric entities associated with the endpoints of the edge, if not both vertices are located on the same entity.

But experiments showed that in such cases too often an intermediate entity would have been the correct one and therefore this approach was discarded.

2 4 6 8

20 40 60 80 100

number of refinements

timeinsec

with optimization without optimization

Figure 5.10. Model: NACA wing section. The times were measured after each refinement around the geometry.

5.3 Geometric error

The main goal of boundary smoothing is to reduce the geometric error and therefore get a better approximation of the simulated model. To verify, that the presented and implemented mesh refinement techniques are effective in that sense, the evolution of the geometric error with successive refinement steps was investigated in different test cases.

In each of the test runs, the initial geometric error was computed and compared to the geometric error after each refinement step. For the computation of the geometric error, one recursion step was used.

(38)

28 CHAPTER 5. OPTIMIZATIONS AND RESULTS Two quantities were measured during those test cases: local and global geometric error. The local error is the geometric error per cell, as it was defined in the previous chapter. The global error is the accumulated error of all cells on all processes.

5.3.1 Refinement of the airplane front

In this test case the front geometry of the airplane was used, as described before.

Figure 5.11 shows the relevant boundary cells inside the chosen bounding box. Since the original mesh provided already a good approximation of the geometry, a volume of 1e − 9 was chosen as threshold for the geometric error. All cells with a geometric error larger than this threshold were marked for refinement. The plots in figure 5.12 reveal a clearly decreasing geometric error with each refinement step, hence, prove the refinement effective.

This is confirmed by the global error, shown in figure 5.13, which decreases with each refinement step.

Figure 5.11. Model: Airplane - Bounding box region marked on surface mesh

5.3.2 NACA wing section

The second setup used the mesh around the NACA wing, which was gained from uniform refinement, as it was shown in figure 5.4. The initial mesh had 21904 vertices. In each iteration, the geometric error around the geometry was calculated and cells marked for refinement, when the local error is larger as the threshold, which was chosen to be 10⁻⁷. In figure 5.15, one can see how the global error is again reduced after each refinement step. A curious situation occured during the first refinement: Since the refinement method always chooses the longest edge for bisection, it is possible that no additional boundary vertices are introduced,

(39)

5.3. GEOMETRIC ERROR 29

Figure 5.12. Cockpit view: The boundary cells colored according to the local error.

The upper left picture is the error function of the initial mesh. Upper right is the error after the first local refinement. The lower row is after the second and third refinement.

0 1 2 3

1.5 2 2.5 3 3.5

·10⁻⁶

amount of local refinements

globalError

Figure 5.13. The global geometric error of the airplane front after each refinement.

Calculated over the sum of all local errors. Where the local error is the approximated volume between a boundary face and the geometry.

(40)

30 CHAPTER 5. OPTIMIZATIONS AND RESULTS although the cells at the boundary are marked for refinement. This means, all of these cells had their longest edge inside the mesh and not on the boundary and the global error did not decrease, consequently.

As before for the airplane, the geometric error around the wing section is shown in figure 5.14. In these plots, only cells where the volume between the mesh and the geometry is larger than 10⁻⁷ are shown, i. e., cells that were also selected for refinement.

Regarding the computation time of the geometric error: the computation time increased with each refinement step. This is expected behavior, since the number of elements increased with each refinement step and hence the error had to be evaluated more often.

When comparing the time needed for a refinement step to the amount of projected points in each step, as it is done in figure 5.16, it is clearly visible that a larger amount of projected points is responsible for a longer execution time. Hence, the projection is the dominant part of the refinement. Figure 5.17 shows, that a larger number of processors successfully reduces the refinement time, due to the reduced workload per process.

5.4 Simulation around a NACA wing section with

UNICORN and DOLFIN

To use the geometric boundary correction in an application, the flow around the wing section (see section 5.1) was simulated using Unicorn HPC. Unicorn uses the General Galerkin (G2) method, a variant of the finite element method, which will be briefly explained in the next section, before simulation results are discussed.

5.4.1 Finite element method and G2

The finite element method is a popular tool for the discretization of partial differential equations (PDE). For that, the PDE is formulated in weak form and, using integration by parts, transformed into the variational form. The domain is divided into elements, which correspond to the cells of the computational mesh, and basis functions are defined with respect to each cell. These are inserted into the variational form, allowing to assemble a linear system of equations which can subsequently be solved using an iterative solver, e. g., Conjugated Gradients or GMRES.

General Galerkin is a specific variant of this method that uses piecewise linear basis functions in both space and time. Hence, it is particularly useful for the discretization of time dependent systems of differential equations like the Navier- Stokes equations.

(41)

5.4. SIMULATION AROUND A NACA WING SECTION 31

Figure 5.14. Geometric error around the NACA wing. Only cells are displayed where the volume between the face of the triangle and the geometry is larger than 1e − 07. On top is the error of the initial mesh. The other plots show the error after each refinement.

(42)

0 2 4 6

1 2 3 4 5 ·10⁻⁴

globalError

0 2 4 6

0 100 200 300 400 500

timeinsec

np2np3 np4

Figure 5.15. On the left side one can see how the global error is reduced with every local refinement around the geometry. The plot on the right side shows the time needed to calculate this error on 2, 3, and 4 nodes.

0 1 2 3 4 5 6 7 8

20 40 60 80 100 120

timeperrefinementinsec

refinement time

0 1 2 3 4 5 6 7 8

100 200 300 400 500 600 700

projectedpoints

projected points

Figure 5.16. This plot shows the correlation between the necessary time per re- finement on two processes and the amount of projected points in this refinement step.

(43)

5.4. SIMULATION AROUND A NACA WING SECTION 33

2 3 4

40 50 60 70 80 90 100 110 120

number of processes

timeinsec

2nd refinement step: 486 projected points 4th refinement step: 688 projected points

Figure 5.17. Times needed per refinement. The calculations were done on two, three or four processes

5.4.2 Simulation set up

The incompressible Navier-Stokes equations were considered over the domain Ω and within a time interval [0,T]:

˙u + (u · ∇)u + ∇p = g + 2ν∇ · ε(u) in Q = Ω × [0, T ] (5.1)

∇ · u= 0 in Q (5.2)

ˆu(·, 0) = ˆu⁰ in Ω (5.3)

u is the velocity, p is the pressure and ν is the kinematic viscosity.

The discretization of the Navier-Stokes equations using G2 was described by Hoffman and Johnson [1].

Here, (·, ·) denotes the inner product. As boundary conditions no slip conditions were applied, i. e., u = 0 on the boundary.

Figure 5.18 shows the results after the first and the fourth iteration of the adaptive solver. In the picture of the fourth iteration, it can be seen that the representation of the geometry was improved during the refinements. Especially around the leading edge the impovement is visible.

(44)

Figure 5.18. Flow simulation around a wing section using geometric correction during the refinement. Left: result after the first iteration right: result after the fourth iteration

(45)

Chapter 6

Summary and future work

During this project, the three main components, necessary for geometric boundary smoothing in high performance automated finite element software, were developed and implemented: (i) initial projection, (ii) geometric error estimation and (iii) extension of the refinement procedure. For that, two libraries were introduced, to handle the CAD data and provide an interface for important algorithms from geometric modeling.

The basic idea of initial projection was to project boundary points of the mesh onto the geometric boundary of the CAD model and correct the position of these vertices such that they are located exactly on the surface of the model. The execution time of this process was improved by a load balancing step that distributed the workload evenly among the available processing elements. For geometries that are built up from multiple geometric entities, it was necessary to find the correct entity first before the projection could be performed. This initial search was accelerated by using a smaller number of sampling points.

The deviation of the mesh from the geometric boundary was measured by an own quantity for the geometric error. It considers the surface or volume enclosed between geometry and mesh boundary and is computed with a recursive algorithm. The accuracy of the determined error can be improved by a larger number of recursion steps, at the cost of an increased computation time.

To improve the approximation of the geometric boundary, the existing refinement procedure was adapted. Newly created vertices that are on the domain boundaries are now positioned with respect to the geometric model.

However, some ideas remain how to improve the developed and implemented procedures: Up to now, the selection of cells for refinement in the fluid solver Uni- corn is solely based on the error indicators obtained from primal and dual solution.

Consequently, an additional refinement step was necessary, based on the geometric error indicators. Incorporating those two techniques into a single procedure would be the next step.

NURBS can be used in many different ways to describe the geometrical shape of a model. So far, the abilities of LIBGEOM are limited to the description of

35

(46)

36 CHAPTER 6. SUMMARY AND FUTURE WORK geometric boundaries by curves or surfaces. But boundaries can also be described by trimming holes into surfaces using curves or other surfaces. While a rudimentary support for intersecting surfaces is already given today, a representation of very complex geometries, which make extensive use of such modeling techniques, is not possible, yet. Hence, to make the developed library usable in a wider field, support for such complex entities has to be added.

Also, many more and composed geometric entities are included in the complete IGES standard, for example a straight line, for which no support was added, yet.

Such types could be easily transformed into a NURBS representation, but that would require an extensive parser that provides the conversion for all possible geometry types of the standard into those basic entities. Prior to the usage of the geometric data it was therefore necessary to perform a conversion of the geometric model into the NASA-IGES-NURBS-Only standard (see appendix B), which limits the available geometric entity types to a minimum. Including such a parser would make the applicability of the developed software more general and comfortable.

(47)

Bibliography

[1] J. Hoffman and C. Johnson. Computational Turbulent Incompressible Flow:

Applied Mathematics: Body and Soul 4. Applied Mathematics, Body and Soul.

Springer, 2007. ISBN 9783540465317.

[2] S. Fernández-Méndez R. Sevilla and A. Huerta. Nurbs-enhanced finite ele- ment method (nefem). a seamless bridge between cad and fem. Archives of Computational Methods in Engineering (invited), 18(Issue 4):441–484, 2011.

[3] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(39–41):4135 – 4195, 2005. ISSN 0045-7825. doi: 10.1016/j.cma.2004.10.008. URL http://www.

sciencedirect.com/science/article/pii/S0045782504005171.

[4] J.A. Cottrell, T.J.R. Hughes, and Y. Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. Wiley, 2009. ISBN 9780470749098. URL http:

//books.google.se/books?id=9Q9y0Xtz5E4C.

[5] Ashraful Kadir. Geometric Shape Preservation in Adaptive Refinement of Fi- nite Element Meshes. Master thesis, KTH - Royal Institute of Technology, Stockholm, Sweden., 2010.

[6] P. J. Frey and H. Borouchaki. Surface meshing using a geometric error estimate.

International Journal for Numerical Methods in Engineering, 58(2):227–245, 2003. ISSN 1097-0207. doi: 10.1002/nme.766. URL http://dx.doi.org/10.

1002/nme.766.

[7] David Levin Sun-Jeong Kim, Chang-Hun Kim. Surface simplification using a discrete curvature norm. Computers and Graphics, 26(5):657 – 663, 2002.

ISSN 0097-8493. doi: 10.1016/S0097-8493(02)00121-8. URL http://www.

sciencedirect.com/science/article/pii/S0097849302001218.

[8] L. Piegl and W. Tiller. The NURBs Book. Monographs in visual communica- tions. Springer Berlin Heidelberg, 1997. ISBN 9783540615453.

[9] E. Dimas and D. Briassoulis. 3d geometric modelling based on nurbs: a review.

Advances in Engineering Software, 30(9–11):741 – 751, 1999. ISSN 0965-9978.

37

(48)

38 BIBLIOGRAPHY

doi: 10.1016/S0965-9978(98)00110-0. URL http://www.sciencedirect.com/

science/article/pii/S0965997898001100.

[10] US Product Data Association. Iges 5.3, 1996. URL http://www.uspro.org/

documents/IGES5-3_forDownload.pdf.

[11] Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012. ISBN 978-3-642-23098-1. doi: 10.1007/978-3-642-23099-8.

[12] Anders Logg and Garth N. Wells. Dolfin: Automated finite element com- puting. ACM Transactions on Mathematical Software, 37(2), 2010. doi:

10.1145/1731022.1731030.

[13] J. Hoffman N. Jansson and J. Jansson. Parallel adaptive fem cfd. Report num- ber 8, Computation Technology Laboratory, KTH - Royal Institute of Tech- nology, Stockholm, Sweden., 2010. URL http://www.publ.kth.se/trita/

ctl-4/008.

[14] M. Rivara. Mesh refinement processes based on the generalized bisection of simplices. SIAM Journal on Numerical Analysis, 21(3):604–613, 1984. doi:

10.1137/0721042. URL http://epubs.siam.org/doi/abs/10.1137/0721042.

[15] Niclas Jansson. High performance adaptive finite element methods for turbulent fluid flow. Licentiate thesis, KTH - Royal Institute of Technology, Stockholm, Sweden., 2011.

[16] Johan Hoffman, Johan Jansson, Rodrigo V. de Abreu, Cem Degirmenci, Niclas Jansson, Kaspar Mueller, Murtazo Nazarov, and Jeanette H. Spuehler. Uni- corn: Parallel adaptive finite element simulation of turbulent flow and fluid- structure interaction for deforming domains and complex geometry. Computer and Fluids, in press, 2013. doi: 10.1016/j.compfluid.2012.02.003.

[17] Philippe Lavoie. URL http://sourceforge.net/projects/libnurbs.

[18] Philippe Lavoie. An introduction to nurbs++, 1999. URL http://libnurbs.

sourceforge.net/old/user.pdf.

[19] Robert M. Pinkerton Eastman N. Jacobs, Kenneth E. Ward. The characteris- tics of 78 related airfoil sections from test in the variable-density wind tunnel.

(460), 1993.

[20] Niclas Jansson. Adaptive Mesh Refinement for Large Scale Parallel Computing with DOLFIN. Master thesis, KTH - Royal Institute of Technology, Stockholm, Sweden., 2008.

[21] L. Piegl. On nurbs: a survey. Computer Graphics and Applications, IEEE, 11 (1):55–71, 1991. ISSN 0272-1716. doi: 10.1109/38.67702.

(49)

BIBLIOGRAPHY 39

[22] David P. Miller Austin L. Evans. Nasa iges and nasa-iges nurbs-only standard.

Handbook of Grid Generations, page Chapter 31.

[23] J.F. Thompson, B.K. Soni, and N.P. Weatherill. Handbook of Grid Generations.

CRC Press, 1999. ISBN 9780849326875. URL http://books.google.se/

books?id=ImaDT6ijKq4C.

(50)

(51)

Appendix A

Curves and Surfaces

A short overview over the possible representations of a model, using NURBS curves or surfaces, is given in this appendix.

Definition A.1. A curve C is a one dimensional connected point set in a two or three dimensional space with at least piecewise smoothness (see figure A.1).

Definition A.2. A surface S is a two dimensional point set in a three dimensional space, with at least piecewise smoothness (see figure A.2).

−1

0 1 −1

0 0 1

10 20

−20 0 20 40

−40

−20 0 20 40

Figure A.1. A curve in three dimensional space (left) and a curve in two dimensional space

This introduction to NURBS is based Piegl [21] and Dimas and Briassoulis [9].

A.1 Mathematical description

A NURBS curve can be described as a vector-valued, piecewise rational, polynomial function and is defined as:

C(u) = Pn

i=1w_iP_iN_i,p(u) Pn

i=1wiNi,p(u) (A.1)

41

(52)

42 APPENDIX A. CURVES AND SURFACES

Figure A.2. A surface in a three dimensional space

where wi are weights, Pi are control points and Ni,p are normalized B-Spline basis functions of degree p.

A surface is analogously defined in the following way:

S(u, v) = Pn

i=0

Pm

j=0wi,jPi,jNi,p(u)Nj,q(v) Pn

i=0

Pm

j=0wi,jNi,p(u)Nj,q(v) (A.2) A.1.1 Normalized B-Spline basis functions

The B-Spline basis functions are spanned over a knot vector U with r + 1 elements, where r = n + p + 1. As already mentioned, n is the number of control points and p is the degree of the basis functions.

The endpoints of the knot vector are repeated and the multiplicity is dependent on the degree, i.e., for non-uniform and non-periodic B-splines this is p + 1.

The basic functions are then defined over the whole knot vector [u0; ur], but in general it is transformed to the interval from 0 to 1. The recursive definition of the basis functions is:

Ni,0(u) =

(1 if ui≤ u ≤ u_i+1

0 otherwise (A.3)

N_i,p(u) = u − u_i

ui+p− u_iN_i,p−1(u) + u_i+p+1− u

ui+p+1− u_i+1N_i+1,p−1(u) (A.4) For the representation of surfaces, a second knot vector V is used, which has the size s+1, where s = m+q +2. V represents the second dimension in the parameter space. Therefor: m is the number of control points in this direction and q the degree of the basis functions.

A.1.2 NURBS basis functions for curves If we define the NURBS basis functions as:

Ri,p(u) = wiNi,p(u) Pn

j=0wjNj,p(u) (A.5)

(53)

A.1. MATHEMATICAL DESCRIPTION 43 we can reformulate (A.1) to a simpler representation, that only depends on the control points and the defined NURBS basis functions:

C(u) =

n

X

i=0

P_iR_i,p(u) (A.6)

A.1.3 Bivariate NURBS basis functions for surfaces We can rewrite equation A.2 for the representation of surfaces as:

S(u, v) =^Xⁿ

i=0 n

X

j=0

Pi,jRi,p;j,q(u, v) (A.7)

if we define the bivariate NURBS basis functions in the following way:

R_i,p;j,q(u, v) = w_i,jN_i,p(u)Nj,q(v) Pn

r=0

Pm

s=0wr,sNr,p(u)Ns,q(v) (A.8) A comprehensive and well written summary of the properties of the basis functions and NURBS in general was given in [9]. Information on how to compute the derivatives of NURBS can be found in the book of Piegl and Tiller [8].

High performance adaptive finite element modeling of complex CAD geometry

High performance adaptive finite

element modeling of complex CAD

geometry

High performance adaptive finite element

modeling of complex CAD geometry

Abstract

Referat

Contents

Chapter 1

Introduction

Chapter 2

Basics Concepts

2.1 Representation of geometries

2.2 Point Projection and Point Inversion

Chapter 3

The FEniCS project

3.1 DOLFIN and DOLFIN HPC

3.2 Unicorn HPC

Chapter 4

Geometric correction algorithms and

their integration into FEniCS

4.1 Additional Libraries

4.2 Boundary smoothing

4.3 Requirements for the initial mesh

Chapter 5

Optimizations and Results

5.1 Used models

5.2 Performance optimizations

5.3 Geometric error

5.4 Simulation around a NACA wing section with

UNICORN and DOLFIN

Chapter 6

Summary and future work

Bibliography

Appendix A

Curves and Surfaces

A.1 Mathematical description