Meshing for pVersion Finite Element Methods

1

• Meshing for p-Version Finite Element Methods
• M.S. Shephard, X.J. Luo and J.F. Remacle
• Scientific Computation Research Center
• Rensselaer Polytechnic Institute, Troy, NY
  12180-3590, USA
• R.M. OBara and M.W. Beall
• Simmetrix Inc., Clifton Park, NY 12065, USA
• In Collaboration With
• B.A. Szabo
• Washington University, St. Louis, MO 63130, USA
• R. Actis
• Engineering Software Research and Develop Inc

• Outline
• Curvilinear mesh requirements and representation
• Curvilinear mesh generation
• Influence of geometric approximationon p-version
  solution accuracy

2
p-Version Finite Element Method

• Capability of exponential rate of convergence
• Can produce a sequence of solution by increasing
  polynomial order
• Gaining popularity in industry
• p-version meshes
• Coarse mesh needed to maintain computational
  efficiency
• Mesh layout needs strong control and gradation
• Mesh must maintain appropriate level of geometric
  approximation
• Curvilinear Mesh Generation
• Approaches
• Directly generate curved meshes from curved
  geometry domains
• Start with straight-sized, planar faces meshes
  and curve mesh entities on curved boundaries and,
  as needed, mesh interior
• Key issues
• Shape of mesh entities classified on curved
  geometry model boundary
• Shape validity verification of curved mesh
  entities
• Abilities to construct and modify curved mesh
  entities as needed to obtain valid and acceptable
  curved meshes

3
Geometric Approximation Representations

• Mesh representation by Lagrange interpolants
• Lagrange interpolation to curve initially
  straight-sided meshes.
• Worked out for quadratic Lagrange elements
• Higher order that quadratic too messy - both
  geometric control and computations required
• results will show quadratic is not sufficient if
  p-value is raised
• Mesh representation by Bezier polynomial
• Effectively increase the geometric approximation
• order to any order desired
• Focus of current efforts
• Mesh representation by Spline methods
• More numerical stable
• Most geometry is piecewise (NURB or B-spline)

4
Application of Bezier Polynomial in Curved Meshing

• Advantageous properties of Bezier polynomials
• Can be as high a degree as desired
• Convex hull provides smoother and more
  controllable approximation
• Better properties to allow more efficient
  intersection checks
• Derivatives and products of Beziers are also
  Beziers
• Efficient algorithms for degree elevation and
  subdivision
• Technical issues to define mesh entity shape
• Interpolating and/or approximating model geometry
  
• Accounting for geometric modeling systems face
  parametric coordinates periodicity, degeneracy
  and distortion
• Chord length parameterization method for mesh
  edge on model face
• Cord length, d, for n1 interpolation points
  Qk, k0,,n defined as
• Edge parameters associated with locations are t0
  0, tn 1 and

5
Approximating Model Geometry - continuity

• Bezier mesh edges
• Three Bezier control points
  must be collinear to maintain continuity
  through the junction point
• Bezier mesh triangle faces
• All pairs of subtriangles through
• the common edge must be coplanar
• Each pair is an affine map of the
• two domain triangles
• Can be generalized up to continuity

6
Generic Hierarchic Structure of Mesh Entity Shape

• Mesh entity shape has been designed as hierarchic
  structure
• All of the inherited classes share the same
  interface as the base class
• Effective representation and definition of mixed
  order curved meshes
• Mesh entity needs to inflate shape based on
  entities they bound in order to support different
  polynomial orders

Inflating a face from quadratic to cubic
Need to inflate face closure based on bounding
edges
7
Mesh Region Shape Validity Determination

• Traditional validation methods test the Jacobian
  at integration points
• Goal - provide a general validation for Bezier
  Regions Relate Jacobian to the region control
  points and determine its minimum bound
• The Jacobian of a Bezier Region
• Partial derivatives of the region are themselves
  Bezier functions
• Jacobian determinant is defined by box-product of
  partial derivatives
• Since the product of 2 Bezier functions is also a
  Bezier function, the Jacobian determinant is also
  a Bezier function
• In the case of a tetrahedron the function is of
  order 3(p-1), where p is the order of the
  original shape
• Since a Bezier function is bounded by its convex
  hull, the Jacobian determinant function inside
  the region is bounded by the convex hull of its
  control points (which in this case are scalars)
• A region is valid globally if the minimum control
  point of the Jacobian determinant function is gt
  0

8
Testing a Cubic Bezier Tetrahedron

• In the cubic case, each partial derivative is
  composed of 10 control points which define a
  quadratic tetrahedron
• In the case of a fully curved tetrahedron, there
  are 1000 box product terms that make up 84 linear
  relationships
• If there are uncurved edges and faces then the
  number of calculations is reduced

• Vectors that make up the partial derivatives

Vectors that form the partial w/r to ?1
Vectors that form the partial w/r to ?2
Vectors that form the partial w/r to ?3
9
Invalid Region

• The box product terms that compose the Jacobian
  determinant function can be used to determine how
  a region should be corrected

Invalid Tet caused by moving P1 - note
that a (b x c) lt 0

• Jacobian Control point, Jl for a tet region of
  degree d is equal to

10
Correcting an Invalid Region by Shape Manipulation

• Procedure for Correcting an Invalid Region
• For each Jl lt 0
• Identify the min Box Product term contributing to
  Jl
• Identify the region control points involved in
  the the box product vectors that can be moved
• Control Points may be constrained to
• be associated with mesh entities on the model
  boundary
• constrained to prevent other mesh regions from
  becoming invalid
• Identify the min angle to make the Box Product
  Term gt 0
• Determine which control points of region that
  defines the vector should be displaced in order
  to rotate the vector
• If this change would result in invalidating a
  neighboring region - modify that region in order
  to accommodate the shape change
• If the region is still invalid then perform one
  of the following
• Degree elevate the regions shape if possible in
  order to increase the degrees of freedom
  available
• Sub-divide the region in order to refine the mesh
  and introduce more degrees of freedom

11
Correcting an Invalid Entity by Local Modification

• Determine key mesh entities to prevent shape
  movement
• Identify each
• Find the dominating pair of partial vectors
  causing each
• Count the appearance number of each partial
  vectors
• Key entities - for the partial vector(s) which
  has the biggest appearance number, find the mesh
  entities which control points form the vector(s)
• Vector appears most frequently in the
  counting
• Mesh edge is the key entity

12
Correcting an Invalid Entity by Local Modification

• Analyze current situation to determine most
  effective operations
• Determine how much space needed for the key
  entity to fix the invalidity
• Small re-curving
• Large - others

13
Correcting an Invalid Entity by Local Modification

• Determine the existence of small edge length,
  face area or region volume in the neighboring of
  key mesh entity
• Exists apply edge, face, region collapse to
  produce more space
• Appropriate apply swap, split or compound
  operation

split
14
Quality of Curved Mesh for p-Version Method

• Quality of curved mesh is still an open issue
• Minimum determinant of Jacobian
• Determinant of Jacobian variation inside one
  element
• Normalized Minimum determinant of Jacobian
• Geometric approximation error
• Two main difficulties
• Lack of mathematical proof
• Hard to relate quality to finite element
• solution accuracy
• Example
• Compare two valid curved meshes
• based on the same geometric model

15
Quality of Curved Mesh for p-Version Method

• Quadratic mesh geometry
• Case (a) focusedon maximizingthe min. Jacobian
  -leads to strong interior meshentity curving

(a)
(b)
16
Quality of Curved Mesh for p-Version Method

• Meshes comparison
• Volume exact volume of the model 3.56241E-04
• Maximum distance deviation
• The maximum distance between the sample point of
  each mesh entity classified on curved model
  boundary and its corresponding closest point in
  the boundary
• Normalized with respect to the longest
  diagonal edge length of the model box
• The volume variation between linear mesh and
  curved mesh is small
• The geometric approximation has been improved by
  using curved mesh especially for case (b)
• Case (a) introduces excessive interior distortion

17
Curvilinear Mesh Refinement

• By maintaining the original curved shape, Bezier
  curved mesh supports generic refinement up to any
  order
• nth order Bezier mesh edge refinement at location
  t
• 
  
• nth order Bezier mesh triangle face refinement at
  location
• where
• are corresponding
• control points of triangle

18
Curvilinear mesh refinement

• nth Bezier tetrahedron region refinement at
  location
• where and
  are control points of original region
  
• Examples
• Edge split

19
Curvilinear mesh refinement

• Face split
• Region split
• Need a capability that produces similar effects
  for swaps

20
Singularity isolation

• Generalized advancing layers method to isolate
  vertex and edge singularities
• Generate long elements aligned with singular
  edges
• Strong geometric grading in the radial direction
  is employed
• Special function is applied to construct meshes
  for uncovered domain which may appears when
  singular edges and singular vertices meet with
  each other
• Examples
• Isolating the singularity around a crack

21
Singularity isolation

• Three edge singularities meeting at a vertex
  singularity

22
Examples

23
Examples
24
Examples
Mesh curved Using Quadratic Beziers
Mesh curved Using Cubic Beziers
25
Geometric Approximation in p-Version Method

• Investigate the influence of geometric
  approximation on the solution accuracy in
  p-version finite element method
• Model problem - Infinite plane with an elliptical
  hole
• Uniform tensile stress in vertical direction
• Only one quarter domain has been investigated
  double symmetry
• Traction (natural) boundary condition on edges BC
  and CD
• Symmetric essential boundary condition on edges
  AB and DE
• Polynomial order (p) varies from 1 to 10
• Geometric approximation order (q) varies from 1
  to 4

26
Error Norms

• Energy norm - must control energy norm to control
  pollution errors
• For this linear problem the easy to compute
  potential ( )and to also show it is the
  negative of strain energy
• Finite element potential energy ( ) is the
  product of load vector and finite element
  displacements over the loaded boundary
• Using this the relative error in energy norm
• norm - error in peak stress
• Exact value at vertex A is
• Direct computation of the finite element maximum
  stress
• Search for is conducted over the Gauss
  Quadrature points and the vertices of each
  element
• Definition of relative error in maximum stress

27
Test Problems

• Two models with different parameter are
  selected
• where m is a third parameter to describe the
  inner ellipse. m 0 corresponds to a circle and
  m 1 is a sharp crack
• Isotropic material properties
• Youngs Modulus 1.0
• Poisson ratio 0.3
• Plain strain

28
Meshes and Geometric Approximation Shape

• 
  m 0.25
• Mesh edge geometric approximation shape
• interpolant interpolating points are
  equally or unequally spaced in parametric space
• continuity enforced through vertices A and
  E (can get for even p2 due to symmetries)

29
Geometric Approximation Shape around Vertex A

• m 0.25
• 
  q 2
• q 3
  q 4

area shown in plots
30
Error in Energy Norm for Model with m 0.25

• 
  
• (a) shape
  (b) shape
• When , the error in mapping begins to
  dominate the solution error
• Finite element error approaches a limit when p
  increases which is essentially the geometric
  discretization error
• The geometric discretization error is less when
  geometric approximation order (q) increases

31
Error in Maximum Stress for Model with m 0.25

• (a) Interpolant shape
  (b) slope continuous shape
• q 1, the computed maximum stress overestimated
  the exact value at p 10 by relative error 122
• Expected behavior
• Sharp corner exists at point A
• Stress theoretically goes to infinity

32
Shape Result for Model with m 0.9

• (a) Energy norm
  (b) norm
• The differences between the results of and
  shapes are small
• Errors in energy norm decrease when q increases
  but still unacceptable
• The computed stress is far below the exact stress
  value and relative errors are very high at p
  10.
• (85 -75 ) for q 1 to 4

33
Geometric Approximation Shape around Vertex A

• m0.9 q2
• q3
  q4

34
Graded Mesh and Approximation Shape for m0.9

• Use two mesh edges to approximate the ellipse
  with gradation 0.15
• 
  q 2
• q3
  q4

35
Relative Error in Energy Norm for Graded mesh

• 
  
• (a) shape
  (b) shape
• Error in energy norm decrease comparing to
  ungraded mesh, but the difference is not dramatic
• shape produces a slightly better result
  than shape

36
Relative Error in Maximum Stress for Graded Mesh

• (a) shape
  (b) shape
• q 1, the computed maximum stress is unbounded
  as expected

37
Relative Error in Maximum Stress for Graded Mesh

• shapes underestimate the exact stress value
  for lower p and overestimate the exact value for
  higher p.
• Relative error at p 10
• 56.691 for q 2
• 21.325 for q 3
• 7.557 for q 4
• shapes always underestimate exact stress
  value
• Relative error at p 10
• -30.515 for q 2
• -17.934 for q 3
• -1.517 for q 4
• Error in maximum stress has been greatly improved
  by using curved graded mesh

38
Model Curvature Driven Parameterization

• Equal spaced parameterization for models with m
  0.9
• Geometric approximation shapes bad for
• Geometric approximation errors are still big even
  increasing q
• Results are unacceptable
• Curvature driven unequal spaced parameterization
  method
• Curvature variation around vertex A are high when
• Reduce the geometric approximation error at
  interesting vicinity A

39
Model Curvature Driven Parameterization

• Model edge is
  in parametric space
• Curvature
• Determination of n-1 interpolation points of
  nth order shape
• Unequal spaced interpolation points for model
  with m 0.9

40
Meshes and Shapes for Unequal Spaced Interpolant

• m0.9
  q2
• q3
  q4

41
Results for Unequal Spaced Interpolant

• (a) Energy norm
  (b) norm
• Errors in energy norm decrease comparing to equal
  spaced shape
• The computed stress still always underestimates
  the exact value by relative error 40 for
  quadratic shape
• The computed stress of cubic and quartic shapes
  highly overestimate the exact stress at p 3
  then decrease to converge to 39 when
  p 10
• 5.234 for q 3 2.641 for q 4

42
Error in Maximum Stress for Model with m 0.25

• q 2 and q 3
• shapes underestimate the exact stress value
  for lower p and overestimate the exact value for
  higher p.
• Relative error at p 10 45 for q 2 , 7.7
  for q 3
• shapes always underestimate exact stress
  value
• Relative error at p 10 16 for q 2,
  -5.0 for q 3
• q 4
• Both and shapes have similar behavior as
  shape of q 2, 3
• Substantial smaller relative error when p 10
• 2.8 for shape
• 0.29 for shape
• Consistent with the results obtained by using
  blending function method

43
Closing Remarks

• p-version finite element method requires careful
  construction of appropriate meshes to fully
  realize its advantageous exponential convergence
  rate properties.
• The application of Bezier polynomial in
  curvilinear mesh generation for p-version finite
  element progresses nicely
• The solution quality of p-version finite element
  method is strongly affected by the geometric
  approximation of curved domain
• Conventional assumption of quadratic geometric
  approximation shape is not adequate in p-version
  finite element method
• Choice of geometric approximation order depends
  on the solution accuracy requirement