-
Notifications
You must be signed in to change notification settings - Fork 63
Curved Meshing
This wiki will describe the current status of curved meshing in CRV, and the current and upcoming work.
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. You may want to use pdflatex twice to build the table of contents in the document.
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.
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:
- 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.
- Synchronize all part boundaries across cores.
- 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)
- Synchronize all part boundaries across cores.
- (optional) Attempt to fix all invalid entities.
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 thus only up to that order is currently supported.
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.
To fix invalid elements or adapt to sizing fields, extensions of standard operations are implemented: refinement, coarsening, swapping, and repositioning. Repositioning only works for second order meshes.
Shape correction (fixing invalid elements) after curving can be done with
ma::Input* shapeFixer = crv::configureShapeCorrection(mesh);
crv::adapt(shapeFixer);
Of note is the input parameters shouldSnap and shouldTransferParametric which should be both true if the geometry file (CAD model) is present.
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 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.
For example, assume that we have a second-order Tet Mesh. Using apf::writeVtkFiles we can write a vtk-file which in addition to showing the vertex nodes also shows the mid-edge nodes. So, instead of straight sided Tets, we will have Tets with broken edges, e.g. the following image shows a straight-sided Tet along with a 10-point Tet in Paraview.
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
- Matrix inversion is one cost of inaccuracy at this time. Rather than use an external tool for inverting matrices, several functions in
crvMath.hare used. These lose accuracy for large matrix inversions, for example, the transformation matrix for an 8th order bezier tetrahedron is 165x165, and the inversion loses about 8 digits of accuracy. This is relevant for transforming Jacobian determinants, where the order of the Jacobian determinant shape is 3*(order-1), and the transformation matrix can get very large. - Only one order is supported at a time. Since shape functions are static in implementation, this approach was chosen. When variable order fields are supported in apf, this will need to be changed.
- Checking validity/measuring quality is the single biggest computational cost at this time. Operations such as coarsening and swapping, where the best quality configuration is chosen (or no configuration is deemed valid, and the operation is not performed) are significantly more expensive in the curved case
- Improve shape correction so ALL invalid elements can be fixed. There is room in coarsening/swapping to tweak these operations to improve their success rate. Repositioning can be examined as well. An optimizer may be necessary at some point.
- Curved Adaptation is NOT supported for G1 Meshes. This needs to be looked at, but the problem largely stems from that variable order meshes are not supported, and as such G1 triangles on boundary results in the G1 point structure on the interior, and there is no current support for handling the three additional points (probably not hard to implement).
- Currently, APF Fields do not have their shape written out into the file. As a result, when reading in a quadratic mesh in our format, we cannot distinguish between Lagrange and Bezier shapes.