diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 256e58671..a5307fb7c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -65,6 +65,7 @@ endif() util_exe_func(measureAnisoStats measureAnisoStats.cc) util_exe_func(measureIsoStats measureIsoStats.cc) util_exe_func(visualizeAnisoSizes visualizeAnisoSizes.cc) +test_exe_func(bezierShapeEval bezierShapeEval.cc) if(ENABLE_SIMMETRIX) util_exe_func(runSimxAnisoAdapt runSimxAnisoAdapt.cc) endif() diff --git a/test/bezierShapeEval.cc b/test/bezierShapeEval.cc new file mode 100644 index 000000000..07abf4785 --- /dev/null +++ b/test/bezierShapeEval.cc @@ -0,0 +1,214 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +/* This file contains miscellaneous tests relating to bezier shapes and + * blended bezier shapes + */ + +static apf::Mesh2* makeOneTriMesh(int order, apf::MeshEntity* &ent); +static apf::Mesh2* makeOneTetMesh(int order, apf::MeshEntity* &ent); + + +int main(int argc, char** argv) +{ + MPI_Init(&argc,&argv); + PCU_Comm_Init(); + lion_set_verbosity(1); + if ( argc != 7 ) { + if ( !PCU_Comm_Self() ) { + printf("Usage: %s >\n", argv[0]); + printf(" can only be 2 (for triangles) and 4 (for tets)\n"); + printf(" is the order of bezier\n"); + printf(" can be -1, 0, 1, 2 (-1 means no blending)\n"); + printf(" inquiry point in the parent entity)\n"); + } + MPI_Finalize(); + exit(EXIT_FAILURE); + } + + int type = atoi(argv[1]); + int order = atoi(argv[2]); + int b = atoi(argv[3]); + PCU_ALWAYS_ASSERT_VERBOSE(type == 2 || type == 4, + " can only be 2 or 4!"); + PCU_ALWAYS_ASSERT_VERBOSE(b > -2 && b < 3, + " must be between -1 and 2!"); + + apf::MeshEntity* ent = 0; + apf::Mesh2* m = type == 2 ? makeOneTriMesh(order,ent) : makeOneTetMesh(order,ent); + + double xi0 = atof(argv[4]); + double xi1 = atof(argv[5]); + double xi2 = atof(argv[6]); + + + // set the order + apf::FieldShape* bezierShape = crv::getBezier(order); + int non = bezierShape->getEntityShape(type)->countNodes(); + apf::Vector3 xi(xi0, xi1, xi2); + apf::NewArray vals(non); + + if (b == -1) // regular shape functions + bezierShape->getEntityShape(type)->getValues(m,ent,xi,vals); + else // blended shape functions + { + crv::setBlendingOrder(type, b); + if (type == 2) + crv::BlendedTriangleGetValues(m,ent,xi,vals); + else + crv::BlendedTetGetValues(m,ent,xi,vals); + } + + + printf("shape values are \n"); + for (int i = 0; i < non; i++) { + printf("%d, %E \n", i, vals[i]); + } + + PCU_Comm_Free(); + MPI_Finalize(); +} + +static apf::Mesh2* makeOneTriMesh(int order, apf::MeshEntity* &ent) +{ + apf::Mesh2* mesh = apf::makeEmptyMdsMesh(0, 2, false); + + double vert_coords[3][6] = { + {0.,0.,0., 0., 0., 0.}, + {1.,0.,0., 0., 0., 0.}, + {0.,1.,0., 0., 0., 0.}, + }; + + // each edge is defined by the bounding verts + int edge_info[3][2] = { + {0,1}, + {1,2}, + {2,0} + }; + + apf::MeshEntity* verts[3]; + apf::MeshEntity* edges[3]; + + for (int i = 0; i < 3; i++) { + apf::Vector3 coords(vert_coords[i][0], + vert_coords[i][1], + vert_coords[i][2]); + apf::Vector3 params(vert_coords[i][3], + vert_coords[i][4], + vert_coords[i][5]); + verts[i] = mesh->createVertex(0, coords, params); + } + for (int i = 0; i < 3; i++) { + apf::MeshEntity* down_vs[2] = {verts[edge_info[i][0]], + verts[edge_info[i][1]]}; + edges[i] = mesh->createEntity(apf::Mesh::EDGE, 0, down_vs); + } + + ent = mesh->createEntity(apf::Mesh::TRIANGLE, 0, edges); + + mesh->acceptChanges(); + + apf::changeMeshShape(mesh, crv::getBezier(order),true); + + printf ("Made tri with verts:\n"); + for (int i = 0; i < 3; i++) { + printf("v0: (%e,%e,%e)\n", vert_coords[i][0], vert_coords[i][1], vert_coords[i][2]); + } + + printf ("all nodes of the tri:\n"); + apf::Element* elem = apf::createElement(mesh->getCoordinateField(),ent); + apf::NewArray nodes; + apf::getVectorNodes(elem,nodes); + apf::destroyElement(elem); + for (int i=0; i<(int)nodes.size(); i++) + { + printf("node %d: (%e,%e,%e)\n", i, nodes[i][0], nodes[i][1], nodes[i][2]); + } + return mesh; +} + +static apf::Mesh2* makeOneTetMesh(int order, apf::MeshEntity* &ent) +{ + + apf::Mesh2* mesh = apf::makeEmptyMdsMesh(0, 3, false); + + double vert_coords[4][6] = { + {0.,0.,0., 0., 0., 0.}, + {1.,0.,0., 0., 0., 0.}, + {0.,1.,0., 0., 0., 0.}, + {0.,0.,1., 0., 0., 0.}, + }; + + // each edge is defined by the bounding verts + int const (*edge_info)[2] = apf::tet_edge_verts; + int face_info[4][3] = { + {0,1,2}, + {0,4,3}, + {1,5,4}, + {2,5,3} + }; + + + apf::MeshEntity* verts[4]; + apf::MeshEntity* edges[6]; + apf::MeshEntity* faces[4]; + + for (int i = 0; i < 4; i++) { + apf::Vector3 coords(vert_coords[i][0], + vert_coords[i][1], + vert_coords[i][2]); + apf::Vector3 params(vert_coords[i][3], + vert_coords[i][4], + vert_coords[i][5]); + verts[i] = mesh->createVertex(0, coords, params); + } + for (int i = 0; i < 6; i++) { + apf::MeshEntity* down_vs[2] = {verts[edge_info[i][0]], + verts[edge_info[i][1]]}; + edges[i] = mesh->createEntity(apf::Mesh::EDGE, 0, down_vs); + } + + for (int i = 0; i < 4; i++) { + apf::MeshEntity* down_es[3] = {edges[face_info[i][0]], + edges[face_info[i][1]], + edges[face_info[i][2]]}; + faces[i] = mesh->createEntity(apf::Mesh::TRIANGLE, 0, down_es); + } + + ent = mesh->createEntity(apf::Mesh::TET, 0, faces); + + mesh->acceptChanges(); + + apf::changeMeshShape(mesh, crv::getBezier(order),true); + + printf ("Made tet with verts:\n"); + for (int i = 0; i < 4; i++) { + printf("v0: (%e,%e,%e)\n", vert_coords[i][0], vert_coords[i][1], vert_coords[i][2]); + } + + printf ("all nodes of the tet:\n"); + apf::Element* elem = apf::createElement(mesh->getCoordinateField(),ent); + apf::NewArray nodes; + apf::getVectorNodes(elem,nodes); + apf::destroyElement(elem); + for (int i=0; i<(int)nodes.size(); i++) + { + printf("node %d: (%e,%e,%e)\n", i, nodes[i][0], nodes[i][1], nodes[i][2]); + } + return mesh; +}