Institute for
Numerical Simulation

INS

LOCAL ADAPTIVITY AND QUASI-OPTIMAL DISCRETIZATIONS
FOR MULTIVARIATE ISOGEOMETRIC ANALYSIS

WE PRESENT:

1. Drawback of traditional methods
2. Overview about Isogeometric Analysis
3. Local refinements and non-conformity
4. Quasi-optimal choice of the refinement type
5. Results for single patches
6. Spline Multigrid solver

1. Drawback of traditional methods:

Classical methods for simulation use as input a mesh which is triangular or quadrilateral in 2D and which is tetrahedral or hexahedral in 3D. They have several disadvantages which are related to geometry and hierarchy:

• Those traditional methods require dense meshes not only in order to attain a good accuracy for solving the PDE but also to obtain a good approximation of the physical geometry. Indeed, in order to obtain an accurate representation of the used geometries, an initial mesh is usually taken. For real world CAD objects, a large number of gridpoints are required in order to have a good approximation of the geometries. The resulting polygonal or polyhedral approximation, usually denoted by Ω, is considered as the initial mesh of simulation. Depending on the accuracy required for the PDE simulation, the domain Ω still needs to be further refined by using some a-posteriori error estimator. Hence, an unnecessarily large degree of freedom.
• To overcome the above problem, traditional methods take a very coarse mesh which is a polyhedral or polygonal approximation. That mesh is then repeatedly updated by shifting boundary nodes to the exact geometry. That is, nodes are translated toward the boundary when boundary refinements occur. Such a refinement technique can be easily applied to convex domains such as spheres or cylinders. When the domain has a concave boundary, shifting a boundary node to a curved boundary could cause an interference such as mesh folding because the new boundary node might well be inside some triangle. The resulting mesh interference has to be corrected. Thus, one has to apply a difficult geometric rectification on the fly beside solving the PDE.
• Many numerical analysts are fond of hierarchy on account of its efficiency. In the context of hierarchical methods like multigrid, multilevel or multiscale, one needs a sequence of hierarchical nested finite dimensional spaces such as

  

  Traditional methods suffer from difficult problems in generating such a hierarchical structure because the only available data are usually the finest space. Algorithms for mesh simplification (or geometry decimation) using multiresolution techniques do exist. However, they enable only the sequence of meshes to become coarser and coarser. The sequence of corresponding finite dimensional spaces are in general not nested

Return to content

2. Overview about Isogeometric Analysis:

The most widely used representations in CAD are Bézier, B-spline, Coons, Gordon, and NURBS patches. NURBS entities have become a CAD standard because they can describe the other representations exactly. In addition, with very few control points NURBS enables exact representations of algebraic entities like circular arcs, spheres, conic sections which are important components of CAD assemblies. Hence, NURBS entities become the most supported components in modern CAD exchange standards such as IGES or STEP. The IGA is featured by three good properties: CAD integration: it enables the coupling of simulation and modeling on the same CAD model. Small geometric degree of freedom: it does not increase the degree of freedom to capture geometric accuracy. It keeps the CAD domain of simulation unchanged from starting until finishing the computation. Hierarchical solver: it enables hierarchical solvers like multigrid even for complicated CAD objects. The sequence of nested spaces can be easily generated while projection and restriction operators are fast to compute.

Return to content

3. Local refinements and non-conformity:

The input is a small number of NURBS surfaces or solids M_i defined on parameter domains P_i which can be supposed to be the unit square or the unit cube. The treatment of the PDE carries the problem from the physical domain over to the parameter domain P. The entries of the Jacobian matrix of each NURBS mapping and the one for its inverse can be exactly computed even if the value of the inverse of the NURBS is unknown. The problem is completely solved on the parameter domain. For the spaces of approximation, the original IGA (see for e.g. works of Hughes, Basilevs, Buffa) uses NURBS functions as bases. But here we restrict ourselves to B-splines which are

For the multivariate case where the dimension is d=2 or d=3, we use the following bases

Hence, the space of approximation within a single B-spline domain is given by

Let us note that the main weakness of IGA so far is that it uses global refinements instead of local ones. Global refinements imply that an insertion of a knot entry in one direction spreads along the whole range of the other directions. Not only such a process increases the degree of freedom but the shape regularity condition of the spline segments could be violated. To circumvent such a problem, we allow our discretization to be non-conforming. In addition, the approximation space on one parameter domain P is given by

A typical simulation on a single 2D patch is illustrated in the next figure. On the left figure, the physical NURBS patch is displayed with the exact solution. The middle figure shows the exact solution on the parameter domain. On the right figure we see some non-conforming discretization obtained from a certain sequence of local refinements together with the corresponding computed solution. The main properties of our implementation are summarized as follows:

• Full adaptivity using a-posteriori error estimator,
• Use of local refinements,
• Optimal element distribution,
• Hierarchical solver.

Return to content

4. Quasi-optimal choice of the refinemen type:

We construct an automatic self-adaptive discretization which starts from an initial mesh that is always supposed to be a very coarse tensor product mesh. We use an a-posteriori error estimator which permits to evaluate the numerical errors without knowing the exact solution. Only elements admitting too large errors need to be refined. That avoids the need of global uniform refinements. In order to evaluate the error within an element Q of the mesh on each patch, we use the following estimator for a d-dimensional problem where d=2,3.

where f is the RHS, M describe the NURBS patch parametrization, A is obtained from the Laplace operator, (n_i (Q),k_i(Q)) are the spline properties on the x _i-direction of the element Q and h_i(Q) is the length of the element Q on the x_i-direction. Afterwards, the elements Q having the largest value of the above estimator are refined in order to obtain a new mesh. When an element in the mesh needs to be refined, it is not known beforehand which kind of refinement ought to be applied. There are several possibilities for refinements. One can apply subdivisions on different directions. One can also insert new knot entries. The choice of the type of refinement should dynamically depend on both the solution of the PDE to be solved and the NURBS blocks. We introduced a method for choosing an optimal refinement which would reduce the error most in the next mesh. Optimality is not understood in the strict sense that there is no other discretization which yields a more accurate result. Optimality here is gauged with respect to some constant frames which depend only on the current discretization. Suppose that the current solution in u_h. Consider a larger space E where the solution is u_E. It is evident that u_E is at least as accurate as u_h. We search for the space E which has the largest error reduction. It is clear that a very small deviation of u_E from u_h could produce only a very small error reduction as in

Therefore, we should try to maximize that deviation so that the error is likely to be reduced most. Unfortunately, u_E is unknown and it makes no sense to compute it on a finer mesh for every possible mesh. Therefore, we need to find a way to gauge the deviation. Based on hierarchical spline space decompositions, we use the following estimates

Here, the constants c₁ and c₂ do not depend on the refinement type or position. In addition, the determination of the estimator r_Q amounts to solving a system which is linear, which is symmetric positive definite and which is very small.
Return to content

5. Results for single patches:

We present now a few refinement results on single patches. First, we consider the following exact solutions which correspond respectively to an internal layer and to an internal accumulation. The RHS of the Laplace equation are computed accordingly. Inside some tight positions of the domains, those functions admit unity values. Elsewhere, they decay exponentially to zero. The size of the internal features can be controlled by the parameter ω. The NURBS function which describes the patch is exactly the same as the one in Section 3. The position of the layer can be controlled by the parameter c. Similarly, the one for the internal accumulation by the parameters (a,b).

Repeated applications of the refinement process yield the following sequence of discretizations on the parameter domain. Further results including numerical comparison with global refinements and with uniform discretizatrions as well as 3D simulations can be found in our papers. The second refinement sequence plainly illustrates the geometric situation in Section 1 where there is no need to perform a refinement next to the boundary although the geometry has curvilinear boundaries.

Return to content

6. Spline multigrid linear solver:

The corresponding linear system is solved by using a modified multigrid algorithm. Below, we summarize the two-grid case which uses two nested spaces. The difference from the standard two-grid method is about the restriction and prolongation operators R and P. For the standard two-grid, some templates combining neighboring nodal values are used. Such a nodal combination cannot be used here because there are no nodes inside the elements. The multigrid process consists in applying the two-grid operation several times during the coarse grid correction.

• Use discrete B-splines for the prolongation operator P,
• Use non-uniform spline-wavelets for the restriction operator R,
• Use the same pre-smoothing and post-smoothing method as in standard multigrid,
• Use the same coarse grid correction as in standard multigrid.

Here, we use the Gauss-Seidel iteration for those smoothing steps. We numerically investigated that spline multigrid whose convergence history until the 21-st iteration is displayed in the next Table. We examine also the corresponding error at each iteration as well as the ratio between the errors of two consecutive iterations.

ITERATION	ERROR	RATIO BTW. CONSEC. ERRORS
0	2.558340e+000	nothing
3	2.932692e-005	10.433512
6	4.572813e-008	5.410555
9	7.159492e-010	3.048095
12	2.787349e-011	2.837645
15	1.242582e-012	2.800942
18	5.710444e-014	2.717146
20	7.999520e-015	2.386161
21	3.352464e-015	1.940152

Return to content

Last Update: on March 31st, 2010, in Bonn, by Maharavo Randrianarivony.