Skip to content

Curved Meshing

Jump to bottom
Daniel Zaide edited this page Feb 22, 2016 · 19 revisions

Curved Meshing Support

This wiki will describe the current status of curved meshing in CRV, and the current and upcoming work.

Documentation

Doxygen exists for all header files in the crv directory, using make doc from the build directory. Mathematical details and other notes and thoughts are in bezier.tex which is in the SCOREC/docs repository.

Curved Shapes

Only simplices are supported at this time for 2D (x-y) and 3D meshes. There is no support for surface meshes at this time. 2D (x-y) support exists solely for debugging purposes.

The following curved entities are supported:

  • 2nd order Lagrange, all shapes
  • 2nd-6th order Full Bezier Simplices
  • 2nd-6th order Blended Bezier Simplices
  • 4th order Full G1 tetrahedra (G1 Triangles on the model boundaries, Bezier Tetrahedra elsewhere)
  • 4th order Blended G1 tetrahedra (G1 Triangles on the model boundaries, Blended Tetrahedra elsewhere)
  • 1st order Bezier Simplices (This exists largely to check the infrastructure surrounding bezier, and is identical to 1st order Lagrange)

In theory, up to 18th order Bezier Simplices are supported, however the integration points only go up to 6th order, so there would be an error when trying to measure areas or volumes.

Note: only one order is supported at a time, due to the implementation. Two meshes loaded with different orders will not perform correctly, nor will projecting between different order Bezier shapes.

Curved Meshing Algorithm

Given a linear mesh and a geometry (Parasolid File, analytic geometry definition), boundary edges and triangles can be snapped to the geometry at node locations in parametric space.

The current set of points used are from Chen and Babuska, 1995, and are better for interpolating than evenly spaced locations. Point locations can be found in crvBezierPoints.cc.

For Bezier entities, interpolating point locations are converted to Bezier control points. On the interior, the mesh remains linear, thus only meshes with convex domains can be safely curved at this time.

The algorithm is as follows for Bezier entities:

  1. For each entity on the boundary, snap the nodes to the CAD geometry, resulting in each high order node interpolating the geometry. This is done hierarchically, edges then faces.
  2. Synchronize all part boundaries across cores.
  3. Convert the interpolating points to Bezier control points using the transformation matrix going downward, with internal nodes on the tetrahedron, followed by internal nodes on the triangle and finally edge nodes converted to control points. Shape blending is used for entities on the boundary (see docs)
  4. Synchronize all part boundaries across cores.
  5. (optional) Attempt to fix all invalid entities.

Usage

The main usage is in crv.h, the only installed header at this time. Curving the mesh is done using a MeshCurver object. For example, to switch to second order lagrange,

apf::changeMeshShape(m,apf::getLagrange(2),true);
crv::InterpolatingCurver ic(m,2);
ic.run();

This will work for any interpolating shape, with the object simply snapping to the geometry (the ability to snap to geometry is required.)

For bezier shapes, changeMeshShape is done internally. It is as simple as

crv::BezierCurver bc(mesh,order,blending_order);
bc.run();

where order is the order of Bezier shape, and blending_order is the order of shape blending to be used (see bezier.tex in docs repo for more info), with blending_order = 0 resulting in a full bezier representation.

If blending is used, note that only orders of 1 and 2 exactly preserve linear elements, and that higher than 2, while supported, will take a high-order mesh of linear elements and give strange results.

Shape correction (fixing invalid elements) after curving can be done with

ma::Input* shapeFixer = crv::configureShapeCorrection(mesh);
crv::adapt(shapeFixer);

crv countNumberInvalidElements(mesh) can be used to check for invalid elements, or the shape correction can just be run. There is no current guarantee all invalid elements will be eliminated, and again, more details are in the tex file mentioned.

There is also a way to convert a mesh from a curved interpolating shape to a Bezier shape.

apf::changeMeshShape(m,apf::getLagrange(2),true);
crv::InterpolatingCurver ic(mesh,2);
ic.run();
crv::BezierCurver bc(mesh,order,0);
bc.run();

Meshes can also easily be changed in order this way (no guarantee the resulting mesh has all valid elements), with

crv::BezierCurver bc1(mesh,order,0);
bc1.run();
crv::BezierCurver bc2(mesh,new_order,0);
bc2.run();

Going to a higher order is much more accurate than going downward in order, which uses an approximation everywhere.

Visualization

There is also support for visualization of high-order fields, on linear or curved meshes, through writeCurvedVtuFiles, which uses subdivision to output into current paraview format. This is all in crvVtk.cc with function definitions in crv.h. There is currently no support for IP fields on curved meshes yet, as far as output is concerned.

This works using subdivision. Triangles are divided into more triangles, Tetrahedra are first divided into four hexes, which are then divided into more hexes.

They will also write two additional variables for each entity, detJacobian, the scaled Jacobian determinant, a quality measure, and minDetJacobian, the scaled Jacobian determinant, but the minimum of the entity. This allows for thresholding and filtering in paraview, as well as proper quality assessment.

ie, crv::writeCurvedVtuFiles(mesh,apf::Mesh::TET,4,"curvedMesh") writes a pvtu file and vtu files with each tet subdivided into 4 (each tet split into 4 hexes, each hex split into 2 hexes per direction).

crv::writeCurvedVtuFiles(mesh,apf::Mesh::TRIANGLE,4,"curvedMesh") writes a pvtu file and vtu files with each triangle subdivided into 16 triangles.

Viewing each entity type seperately is a useful debugging tool, since a 3D mesh with just faces subdivided is much smaller than the full subdivided tet mesh, and makes viewing the outer surface much easier.

There is crv::writeControlPointVtuFiles and crv::writeInterpolationPointVtuFiles, which write out the actual node values, with point data showing entity type, allowing for debugging point locations by looking at where the control points are, and what interpolating location their xi corresponds to.

A general approach to looking at invalid elements and determining causes of invalidity is to load the tet mesh with crv::writeCurvedVtuFiles(mesh,apf::Mesh::TET,order+2,"curvedMesh") and then in Paraview, Threshold by minDetJacobian to look at elements with minDetJacobian < 0. Then switch the data to look at to detJacobian to see the scaled Jacobian determinant in each invalid element.

Tests

There are currently 6 main test files, which contain examples of code usage:

  • bezierElevation - This tests the elevation functionality, increasing the order of Bezier shapes. Elevation is used mostly for validity and quality assessment, but can also be used to increase the order of a mesh.
  • bezierMesh - This tests the overall functionality of curving a mesh, contains usage and checks increasing and decreasing orders of a mesh
  • bezierMisc - This tests math and ordering used to support bezier functionality
  • bezierRefine - This tests the refinement operations and shape correction
  • bezierSubdivision - This tests subdivision functionality, which is used in refinement and validity and quality assessment
  • bezierValidity - This tests the validity assessment of curved elements, using a mesh that is invalid at certain orders

Nektar++

Nektar++ is currently modified to call APF for mesh queries. This is achieved using a Nektar++ input XML file with information about the APF mesh and Parasolid Geometry. All spatial location and Jacobian queries are made through APF, with information such as surface normals calculated from the Jacobian matrix.

  • Flow between two plates, domain is a cube.
  • Poisson's Equation on a sphere.
  • Hagen-Poiseuille Flow through a Cylinder
Clone this wiki locally