|
| 1 | +#include <crv.h> |
| 2 | +#include <crvBezier.h> |
| 3 | +#include <crvTables.h> |
| 4 | +#include <crvSnap.h> |
| 5 | +#include <crvMath.h> |
| 6 | +#include <crvBezierShapes.h> |
| 7 | +#include <gmi_analytic.h> |
| 8 | +#include <gmi_null.h> |
| 9 | +#include <apfMDS.h> |
| 10 | +#include <apfMesh2.h> |
| 11 | +#include <apf.h> |
| 12 | +#include <apfShape.h> |
| 13 | +#include <PCU.h> |
| 14 | +#include <lionPrint.h> |
| 15 | +#include <mth.h> |
| 16 | +#include <mth_def.h> |
| 17 | +#include <pcu_util.h> |
| 18 | +#include <ostream> |
| 19 | +/* This file contains miscellaneous tests relating to bezier shapes and |
| 20 | + * blended bezier shapes |
| 21 | + */ |
| 22 | + |
| 23 | +static apf::Mesh2* makeOneTriMesh(int order, apf::MeshEntity* &ent); |
| 24 | +static apf::Mesh2* makeOneTetMesh(int order, apf::MeshEntity* &ent); |
| 25 | + |
| 26 | + |
| 27 | +int main(int argc, char** argv) |
| 28 | +{ |
| 29 | + MPI_Init(&argc,&argv); |
| 30 | + PCU_Comm_Init(); |
| 31 | + lion_set_verbosity(1); |
| 32 | + if ( argc != 7 ) { |
| 33 | + if ( !PCU_Comm_Self() ) { |
| 34 | + printf("Usage: %s <ent_type> <order> <blend_order> <xi_0> <xi_1> <xi_2>>\n", argv[0]); |
| 35 | + printf("<ent_type> can only be 2 (for triangles) and 4 (for tets)\n"); |
| 36 | + printf("<order> is the order of bezier\n"); |
| 37 | + printf("<blend_order> can be -1, 0, 1, 2 (-1 means no blending)\n"); |
| 38 | + printf("<xi_0> <xi_1> <xi_2> inquiry point in the parent entity)\n"); |
| 39 | + } |
| 40 | + MPI_Finalize(); |
| 41 | + exit(EXIT_FAILURE); |
| 42 | + } |
| 43 | + |
| 44 | + int type = atoi(argv[1]); |
| 45 | + int order = atoi(argv[2]); |
| 46 | + int b = atoi(argv[3]); |
| 47 | + PCU_ALWAYS_ASSERT_VERBOSE(type == 2 || type == 4, |
| 48 | + "<ent_type> can only be 2 or 4!"); |
| 49 | + PCU_ALWAYS_ASSERT_VERBOSE(b > -2 && b < 3, |
| 50 | + "<blend_order> must be between -1 and 2!"); |
| 51 | + |
| 52 | + apf::MeshEntity* ent = 0; |
| 53 | + apf::Mesh2* m = type == 2 ? makeOneTriMesh(order,ent) : makeOneTetMesh(order,ent); |
| 54 | + |
| 55 | + double xi0 = atof(argv[4]); |
| 56 | + double xi1 = atof(argv[5]); |
| 57 | + double xi2 = atof(argv[6]); |
| 58 | + |
| 59 | + |
| 60 | + // set the order |
| 61 | + apf::FieldShape* bezierShape = crv::getBezier(order); |
| 62 | + int non = bezierShape->getEntityShape(type)->countNodes(); |
| 63 | + apf::Vector3 xi(xi0, xi1, xi2); |
| 64 | + apf::NewArray<double> vals(non); |
| 65 | + |
| 66 | + if (b == -1) // regular shape functions |
| 67 | + bezierShape->getEntityShape(type)->getValues(m,ent,xi,vals); |
| 68 | + else // blended shape functions |
| 69 | + { |
| 70 | + crv::setBlendingOrder(type, b); |
| 71 | + if (type == 2) |
| 72 | + crv::BlendedTriangleGetValues(m,ent,xi,vals); |
| 73 | + else |
| 74 | + crv::BlendedTetGetValues(m,ent,xi,vals); |
| 75 | + } |
| 76 | + |
| 77 | + |
| 78 | + printf("shape values are \n"); |
| 79 | + for (int i = 0; i < non; i++) { |
| 80 | + printf("%d, %E \n", i, vals[i]); |
| 81 | + } |
| 82 | + |
| 83 | + PCU_Comm_Free(); |
| 84 | + MPI_Finalize(); |
| 85 | +} |
| 86 | + |
| 87 | +static apf::Mesh2* makeOneTriMesh(int order, apf::MeshEntity* &ent) |
| 88 | +{ |
| 89 | + apf::Mesh2* mesh = apf::makeEmptyMdsMesh(0, 2, false); |
| 90 | + |
| 91 | + double vert_coords[3][6] = { |
| 92 | + {0.,0.,0., 0., 0., 0.}, |
| 93 | + {1.,0.,0., 0., 0., 0.}, |
| 94 | + {0.,1.,0., 0., 0., 0.}, |
| 95 | + }; |
| 96 | + |
| 97 | + // each edge is defined by the bounding verts |
| 98 | + int edge_info[3][2] = { |
| 99 | + {0,1}, |
| 100 | + {1,2}, |
| 101 | + {2,0} |
| 102 | + }; |
| 103 | + |
| 104 | + apf::MeshEntity* verts[3]; |
| 105 | + apf::MeshEntity* edges[3]; |
| 106 | + |
| 107 | + for (int i = 0; i < 3; i++) { |
| 108 | + apf::Vector3 coords(vert_coords[i][0], |
| 109 | + vert_coords[i][1], |
| 110 | + vert_coords[i][2]); |
| 111 | + apf::Vector3 params(vert_coords[i][3], |
| 112 | + vert_coords[i][4], |
| 113 | + vert_coords[i][5]); |
| 114 | + verts[i] = mesh->createVertex(0, coords, params); |
| 115 | + } |
| 116 | + for (int i = 0; i < 3; i++) { |
| 117 | + apf::MeshEntity* down_vs[2] = {verts[edge_info[i][0]], |
| 118 | + verts[edge_info[i][1]]}; |
| 119 | + edges[i] = mesh->createEntity(apf::Mesh::EDGE, 0, down_vs); |
| 120 | + } |
| 121 | + |
| 122 | + ent = mesh->createEntity(apf::Mesh::TRIANGLE, 0, edges); |
| 123 | + |
| 124 | + mesh->acceptChanges(); |
| 125 | + |
| 126 | + apf::changeMeshShape(mesh, crv::getBezier(order),true); |
| 127 | + |
| 128 | + printf ("Made tri with verts:\n"); |
| 129 | + for (int i = 0; i < 3; i++) { |
| 130 | + printf("v0: (%e,%e,%e)\n", vert_coords[i][0], vert_coords[i][1], vert_coords[i][2]); |
| 131 | + } |
| 132 | + |
| 133 | + printf ("all nodes of the tri:\n"); |
| 134 | + apf::Element* elem = apf::createElement(mesh->getCoordinateField(),ent); |
| 135 | + apf::NewArray<apf::Vector3> nodes; |
| 136 | + apf::getVectorNodes(elem,nodes); |
| 137 | + apf::destroyElement(elem); |
| 138 | + for (int i=0; i<(int)nodes.size(); i++) |
| 139 | + { |
| 140 | + printf("node %d: (%e,%e,%e)\n", i, nodes[i][0], nodes[i][1], nodes[i][2]); |
| 141 | + } |
| 142 | + return mesh; |
| 143 | +} |
| 144 | + |
| 145 | +static apf::Mesh2* makeOneTetMesh(int order, apf::MeshEntity* &ent) |
| 146 | +{ |
| 147 | + |
| 148 | + apf::Mesh2* mesh = apf::makeEmptyMdsMesh(0, 3, false); |
| 149 | + |
| 150 | + double vert_coords[4][6] = { |
| 151 | + {0.,0.,0., 0., 0., 0.}, |
| 152 | + {1.,0.,0., 0., 0., 0.}, |
| 153 | + {0.,1.,0., 0., 0., 0.}, |
| 154 | + {0.,0.,1., 0., 0., 0.}, |
| 155 | + }; |
| 156 | + |
| 157 | + // each edge is defined by the bounding verts |
| 158 | + int const (*edge_info)[2] = apf::tet_edge_verts; |
| 159 | + int face_info[4][3] = { |
| 160 | + {0,1,2}, |
| 161 | + {0,4,3}, |
| 162 | + {1,5,4}, |
| 163 | + {2,5,3} |
| 164 | + }; |
| 165 | + |
| 166 | + |
| 167 | + apf::MeshEntity* verts[4]; |
| 168 | + apf::MeshEntity* edges[6]; |
| 169 | + apf::MeshEntity* faces[4]; |
| 170 | + |
| 171 | + for (int i = 0; i < 4; i++) { |
| 172 | + apf::Vector3 coords(vert_coords[i][0], |
| 173 | + vert_coords[i][1], |
| 174 | + vert_coords[i][2]); |
| 175 | + apf::Vector3 params(vert_coords[i][3], |
| 176 | + vert_coords[i][4], |
| 177 | + vert_coords[i][5]); |
| 178 | + verts[i] = mesh->createVertex(0, coords, params); |
| 179 | + } |
| 180 | + for (int i = 0; i < 6; i++) { |
| 181 | + apf::MeshEntity* down_vs[2] = {verts[edge_info[i][0]], |
| 182 | + verts[edge_info[i][1]]}; |
| 183 | + edges[i] = mesh->createEntity(apf::Mesh::EDGE, 0, down_vs); |
| 184 | + } |
| 185 | + |
| 186 | + for (int i = 0; i < 4; i++) { |
| 187 | + apf::MeshEntity* down_es[3] = {edges[face_info[i][0]], |
| 188 | + edges[face_info[i][1]], |
| 189 | + edges[face_info[i][2]]}; |
| 190 | + faces[i] = mesh->createEntity(apf::Mesh::TRIANGLE, 0, down_es); |
| 191 | + } |
| 192 | + |
| 193 | + ent = mesh->createEntity(apf::Mesh::TET, 0, faces); |
| 194 | + |
| 195 | + mesh->acceptChanges(); |
| 196 | + |
| 197 | + apf::changeMeshShape(mesh, crv::getBezier(order),true); |
| 198 | + |
| 199 | + printf ("Made tet with verts:\n"); |
| 200 | + for (int i = 0; i < 4; i++) { |
| 201 | + printf("v0: (%e,%e,%e)\n", vert_coords[i][0], vert_coords[i][1], vert_coords[i][2]); |
| 202 | + } |
| 203 | + |
| 204 | + printf ("all nodes of the tet:\n"); |
| 205 | + apf::Element* elem = apf::createElement(mesh->getCoordinateField(),ent); |
| 206 | + apf::NewArray<apf::Vector3> nodes; |
| 207 | + apf::getVectorNodes(elem,nodes); |
| 208 | + apf::destroyElement(elem); |
| 209 | + for (int i=0; i<(int)nodes.size(); i++) |
| 210 | + { |
| 211 | + printf("node %d: (%e,%e,%e)\n", i, nodes[i][0], nodes[i][1], nodes[i][2]); |
| 212 | + } |
| 213 | + return mesh; |
| 214 | +} |
0 commit comments