SCOREC
MeshAdapt Snapping with Quadratic Triangles
When running MeshAdapt with quadratic triangles, I get the following error:
mesh->getType(e) == apf::Mesh::TET failed at /users/creanj/lore/build/core_simmetrix/ma/maShapeHandler.cc + 63
signal (6): Aborted
while loading /lore/creanj/error_estimate2/tmp.jl, in expression starting on line 51
gsignal at /lib64/libc.so.6 (unknown line)
abort at /lib64/libc.so.6 (unknown line)
PCU_Assert_Fail at /users/creanj/lore/build/core_simmetrix/pcu/pcu_util.c:16
getQuality at /users/creanj/lore/build/core_simmetrix/ma/maShapeHandler.cc:63
getQuality at /users/creanj/lore/build/core_simmetrix/ma/maShapeHandler.cc:64
collectBadElements at /users/creanj/lore/build/core_simmetrix/ma/maSnapper.cc:86
trySnapping at /users/creanj/lore/build/core_simmetrix/ma/maSnapper.cc:114
run at /users/creanj/lore/build/core_simmetrix/ma/maSnapper.cc:260
apply at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:667
applyLocallyWithModification at /users/creanj/lore/build/core_simmetrix/apf/apfCavityOp.cc:44
applyToDimension at /users/creanj/lore/build/core_simmetrix/apf/apfCavityOp.cc:110
applyOperator at /users/creanj/lore/build/core_simmetrix/ma/maOperator.cc:49
snapAllVerts at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:685
snapOneRound at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:778
snapTaggedVerts at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:799
snap at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:815
adapt at /users/creanj/lore/build/core_simmetrix/ma/ma.cc:38
runMA at /users/creanj/.julia/v0.6/PumiInterface/src/funcs1.cc:1412
runMA at /users/creanj/.julia/v0.6/PumiInterface/src/apf.jl:1319 [inlined]
This case works correctly with a linear mesh.
I'm not sure if we support improving the geometric approximation ('snapping') with quadratic meshes. @mortezah
@JaredCrean2: Can you send me the call stack? Or better yet, can you send me your example code/ meshes/models?
There is minimal reproducible example in /users/creanj/lore/test_ma_quadratic. The line the triggers the error ma::runUniformRefinement(m) in test_init2.cc. If you run
./test_init2 circlemesh_linear.smb circle.smd`
a straight-sided mesh is uniformly refined (paraview file output_adapted). If you run
./test_init2 circlemesh_quadratic.smb circle.smd`
a curved mesh is used as input and the error occurs. Here is the backtrace:
#0 0x00007fffeba7b207 in raise () from /lib64/libc.so.6
#1 0x00007fffeba7c8f8 in abort () from /lib64/libc.so.6
#2 0x00007fffece283c8 in PCU_Assert_Fail (
msg=msg@entry=0x7fffffff1f20 "mesh->getType(e) == apf::Mesh::TET failed at /users/creanj/lore/build/core_simmetrix/ma/maShapeHandler.cc + 63 \n")
at /users/creanj/lore/build/core_simmetrix/pcu/pcu_util.c:16
#3 0x00007fffee925ab8 in ma::QuadraticHandler::getQuality (e=0x47b,
this=0x8bb6d0)
at /users/creanj/lore/build/core_simmetrix/ma/maShapeHandler.cc:63
#4 0x00007fffee9457f5 in ma::QuadraticHandler::getQuality (this=0x8bb6d0,
e=0x47b) at /users/creanj/lore/build/core_simmetrix/ma/maShapeHandler.cc:64
#5 0x00007fffee94bd23 in ma::collectBadElements (a=a@entry=0x8bb1b0, es=...,
normals=..., bad=...)
at /users/creanj/lore/build/core_simmetrix/ma/maSnapper.cc:86
#6 0x00007fffee94bfb6 in ma::trySnapping (adapter=0x8bb1b0, tag=0x8c32f0,
vert=0x299, badElements=...)
at /users/creanj/lore/build/core_simmetrix/ma/maSnapper.cc:114
#7 0x00007fffee94d737 in ma::Snapper::run (this=this@entry=0x7fffffff6018)
at /users/creanj/lore/build/core_simmetrix/ma/maSnapper.cc:260
#8 0x00007fffee93f66d in ma::SnapAll::apply (this=0x7fffffff5ff0)
at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:667
#9 0x00007fffed67a36e in apf::CavityOp::applyLocallyWithModification (
this=0x7fffffff5f70, d=0)
at /users/creanj/lore/build/core_simmetrix/apf/apfCavityOp.cc:44
#10 0x00007fffed67a935 in apf::CavityOp::applyToDimension (
this=this@entry=0x7fffffff5f70, d=0)
at /users/creanj/lore/build/core_simmetrix/apf/apfCavityOp.cc:110
#11 0x00007fffee93496f in ma::applyOperator (a=a@entry=0x8bb1b0,
o=o@entry=0x7fffffff5ff0)
at /users/creanj/lore/build/core_simmetrix/ma/maOperator.cc:49
#12 0x00007fffee93e728 in ma::snapAllVerts (a=a@entry=0x8bb1b0,
t=t@entry=0x8c32f0, isSimple=isSimple@entry=false,
successCount=@0x7fffffff61a8: 0)
at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:685
#13 0x00007fffee93e8d7 in ma::snapOneRound (a=a@entry=0x8bb1b0,
t=t@entry=0x8c32f0, isSimple=isSimple@entry=false,
successCount=@0x7fffffff61a8: 0)
at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:778
#14 0x00007fffee93e930 in ma::snapTaggedVerts (a=a@entry=0x8bb1b0,
tag=0x8c32f0) at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:799
#15 0x00007fffee93e9ab in ma::snap (a=a@entry=0x8bb1b0)
at /users/creanj/lore/build/core_simmetrix/ma/maSnap.cc:815
#16 0x00007fffee925dc7 in ma::adapt (in=0x8ed460)
at /users/creanj/lore/build/core_simmetrix/ma/ma.cc:38
#17 0x00000000004014d5 in main (argc=3, argv=0x7fffffff6568)
at /users/creanj/.julia/v0.6/PumiInterface/src/test_init2.cc:64
I peeked at the source code, and it looks like there is no quality metric for quadratic triangles, but there is for quadratic tets.
@JaredCrean2
I don't think we ever use ma::adapt for anything but linear meshes.
For higher order meshes one has to use crv::adapt. To do that the following steps are required:
- the mesh has to be converted to a Bezier mesh
- crv::adapt can be set up pretty much in the same way as ma::adapt is set up
- call to crv::adapt
- for visualization of the meshes to VTK different functions calls are needed (see https://github.com/SCOREC/core/wiki/Curved-Meshing)
I have updated your example with the above changes and it works fine (see below).
Let me know if you have any questions.
#include <iostream>
#include <fstream>
#include <string.h>
#include <apf.h>
#include <gmi.h>
#include <gmi_mesh.h>
#include <gmi_null.h>
#include <apfMDS.h>
#include <apfMesh2.h>
#include <PCU.h>
//#include <apfNumbering.h>
//#include <apfShape.h>
#include <ma.h>
#include <crv.h>
#include "mpi.h"
#include <gmi_sim.h>
#include <SimUtil.h>
// prototype for constructing minimal reproducable examples, using gmi_sim
int main (int argc, char** argv)
{
if (argc < 2 || argc > 3)
{
std::cerr << "Usage: test_init mesh_name.smb [geometry_file]" << std::endl;
return 1;
}
const char* meshname = argv[1];
char geoname[512];
if (argc == 3)
strcpy(geoname, argv[2]);
else
strcpy(geoname, ".null");
std::cout << "Entered init\n" << std::endl;
std::cout << "geoname = " << geoname << std::endl;
MPI_Init(0,NULL); // initilize MPI
PCU_Comm_Init(); // initilize PUMI's communication
// load mesh using null geometry
gmi_register_null();
gmi_register_mesh();
std::cout << "initializing geo_sim" << std::endl;
Sim_readLicenseFile(0);
gmi_sim_start();
gmi_register_sim();
std::cout << "loading geometric model" << std::endl;
gmi_model* g = gmi_load(geoname);
std::cout << "finished loading geometric model" << std::endl;
apf::Mesh2* m = apf::loadMdsMesh(g, meshname );
std::cout << "finished loading mesh" << std::endl;
apf::reorderMdsMesh(m);
// convert to second order Bezier
int order = 2;
crv::BezierCurver bc(m, order, 0);
bc.run();
// write the curve meshes before adapt
crv::writeCurvedVtuFiles(m, apf::Mesh::TRIANGLE, order + 2, "bezier_before_adapt");
crv::writeCurvedWireFrame(m, order + 10, "bezier_before_adapt");
// setup uniform refine and call crv::adapt
ma::Input* in = ma::configureUniformRefine(m);
crv::adapt(in);
// write the curve meshes after adapt
crv::writeCurvedVtuFiles(m, apf::Mesh::TRIANGLE, order + 2, "bezier_after_adapt");
crv::writeCurvedWireFrame(m, order + 10, "bezier_after_adapt");
return 0;
}
From ma::adapt(...), would it be possible to detect that the mesh is curved, print a message, and call abort?
I remember (I might be wrong though) in the past that was actually the case. I will look into it.
But, yes we should be able to detect the order of the mesh and abort ma::adpat if the order is greater than 1, and maybe tell the user to use crv::adapt instead.
I am finally getting back to this, and found a problem with the conversion to Bezier representation. The files are in /users/creanj/lore/test_ma_quadratic2. The domain is the unit circle, and I wrote some code to check the position of coordinate nodes on the boundary by printing out the distance from the origin. Before converting to Bezier all points are on the unit circle, but afterwards the mid-edge nodes are not.
Running the code as ./test_init3 circlemesh_quadratic.smb circle.smd, I get the output
ntered init
geoname = circle.smd
initializing geo_sim
loading geometric model
finished loading geometric model
finished loading mesh
before Bezier
examining vertices
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
examining edges
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
after Bezier
examining vertices
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
r = 1
examining edges
r = 1.00655
r = 1.00688
r = 1.00514
r = 1.00622
r = 1.0068
r = 1.00746
r = 1.00607
r = 1.00943
r = 1.00758
r = 1.00716
r = 1.00773
r = 1.00694
r = 1.00767
r = 1.00722
r = 1.00803
r = 1.00803
r = 1.0081
r = 1.00751
r = 1.00684
r = 1.00765
r = 1.0084
r = 1.00693
r = 1.00766
r = 1.00576
r = 1.00692
r = 1.01009
When the mesh is converted to Bezier the getPoint api will return the position of the control points. Control points are the same as interpolation points only for the vertices but they are not interpolating at other nodes, such as mid edge nodes. That is why you see values different than 1 for edges.
If you want to get the actual coordinates of a point at the mid edge you should do the following for Bezier meshes
apf::MeshElement* me = apf::createMeshElement(mesh, edge);
apf::Vector3 xi(0.,0.,0.);
apf::Vector3 x;
apf::mapLocalToGlobal(me, xi, x);
apf::destroyMeshElement(me);
after the above code x will hold the actual (physical) coordinates of the mid point.
Thanks for the explanation. If getPoint does not return the coordinates of a mesh node, I am going to have to review a lot of code for correctness. Is getPoint returning the control points going to be the final API, or is this an intermediate step towards the final API?
Also, how does setPoint work for Bezier meshes (if I have the xyz coordinates I want to move a mid-edge node to)?
I think getPoint always returned the position of control points for Bezier meshes. Note that if you have Quadratic Lagrange mesh, getPoint would actually return the (physical) position of the midpoints.
If your mesh is Bezier and you have the actual position of the mid-point, there is a way to get the corresponding control point and then set it using setPoint.
The other thing you could covert your Bezier mesh to an interpolation mesh (i.e. quadratic Lagrange). Make your modifications on that mesh and then covert it back to Bezier. (The only reason you want to convert back to Bezier is to be able to run crv::adapt)
You can see an example of using InterpolationCurver here https://github.com/SCOREC/core/blob/d9eeb2b028afd6182b4d3bc1063d3aae3a348760/test/curvetest.cc#L93
I did an experiment yesterday where I keep the mesh as Lagrange most of the time (so my r-adaptation code works), and convert back and forth to Bezier when doing h-adaptation for 2nd order coordinate fields, and it seems to work. I think I will stick with this unless you see any problems with it.
Also, the docs for setPoint and getPoint could be clarified to mention that Bezier meshes behave differently (currently they refer to the "spatial coordinates" of a node).