Branch data Line data Source code
1 : :
2 : :
3 : : #include "meshkit/MKCore.hpp"
4 : : #include "meshkit/SurfaceFacetMeshReader.hpp"
5 : : #include "meshkit/CurveFacetMeshReader.hpp"
6 : : #include "meshkit/ModelEnt.hpp"
7 : : #include <iostream>
8 : : #include <iGeom.h>
9 : : #include <cmath>
10 : :
11 : : #include "moab/Core.hpp"
12 : :
13 : :
14 : : namespace MeshKit
15 : : {
16 : :
17 : : // Output mesh types for this class
18 : : moab::EntityType SurfaceFacetMeshReader_types[] = {moab::MBVERTEX, moab::MBEDGE, moab::MBTRI, moab::MBMAXTYPE};
19 : 42 : const moab::EntityType* SurfaceFacetMeshReader::output_types()
20 : 42 : { return SurfaceFacetMeshReader_types; }
21 : :
22 : 2 : SurfaceFacetMeshReader::SurfaceFacetMeshReader(MKCore *mk_core, const MEntVector &ments)
23 : 2 : : MeshScheme(mk_core, ments)
24 : : {
25 : 2 : mk = mk_core;
26 : 2 : }
27 : :
28 : 0 : SurfaceFacetMeshReader::~SurfaceFacetMeshReader()
29 : : {
30 [ # # ]: 0 : }
31 : :
32 : 2 : void SurfaceFacetMeshReader::setup_this()
33 : : {
34 : : //set the faceting parameters to default values if they are not already set
35 : :
36 : : //set the faceting tolerance
37 [ - + ]: 2 : if(!facet_tol) facet_tol = 1e-4;
38 : :
39 : : //set the geom_resabs value
40 [ - + ]: 2 : if(!geom_res) geom_res = 1e-6;
41 : :
42 : : //create a solid curve mesher
43 [ + - ][ + - ]: 2 : CurveFacetMeshReader *cfmr = (CurveFacetMeshReader*) mk_core()->construct_meshop("CurveFacetMeshReader");
[ + - ]
44 : :
45 : : //set the parameters of the curvemesher to match those of the surface mesher
46 : 2 : cfmr->set_mesh_params(facet_tol,geom_res);
47 : :
48 [ + - ][ + - ]: 18 : for(MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++)
[ + + ]
49 : : {
50 [ + - ]: 16 : ModelEnt *me = mit->first;
51 : : //do an initial check that the model entity is of the correct dimension
52 [ + - ][ - + ]: 16 : if(me->dimension() != 2)
53 : : {
54 [ # # ][ # # ]: 0 : std::cout << "Found a model entity with an incorrect dimension for this mesher" << std::endl;
55 [ # # ]: 0 : mentSelection.erase(mit);
56 : 0 : continue;
57 : : }
58 : :
59 : : //make sure that the model entity's children have been meshed
60 [ + - ]: 16 : MEntVector children;
61 [ + - ]: 16 : me->get_adjacencies(1, children);
62 : : //std::cout << "Number of children found: " << children.size() << std::endl;
63 : :
64 [ + - ][ + - ]: 76 : for(MEntVector::iterator j = children.begin(); j != children.end(); j++)
[ + + ]
65 : : {
66 [ + - ]: 60 : ModelEnt *child = *j;
67 [ + - ][ + + ]: 60 : if(child->is_meshops_list_empty()) cfmr-> add_modelent(child);
[ + - ]
68 : : }
69 : :
70 : 16 : }
71 : :
72 : : //have the children be meshed first
73 : 2 : mk_core()->insert_node(cfmr, this, mk_core()->root_node());
74 : :
75 : 2 : }
76 : :
77 : 2 : void SurfaceFacetMeshReader::execute_this()
78 : : {
79 : :
80 [ + - ][ + - ]: 18 : for(MEntSelection::iterator mit = mentSelection.begin(); mit !=mentSelection.end(); mit++)
[ + + ]
81 : : {
82 [ + - ]: 16 : ModelEnt *me = mit->first;
83 : :
84 : : //---------------------------------//
85 : : // FACET //
86 : : //---------------------------------//
87 [ + - ]: 16 : facet(me);
88 : :
89 [ + - ]: 16 : me->set_meshed_state(COMPLETE_MESH);
90 : : }
91 : 2 : }
92 : :
93 : 16 : void SurfaceFacetMeshReader::facet(ModelEnt *surf)
94 : : {
95 : :
96 [ + - ]: 16 : iGeom::EntityHandle h = surf->geom_handle();
97 [ + - ]: 16 : iMesh::EntitySetHandle sh = IBSH(surf->mesh_handle());
98 : :
99 : : //get the facets for this surface/ ref_face
100 [ + - ]: 16 : std::vector<double> pnts;
101 [ + - ]: 32 : std::vector<int> conn;
102 [ + - ][ + - ]: 16 : iGeom::Error errval = mk->igeom_instance()->getFacets(h,facet_tol,pnts,conn);
103 [ - + ][ # # ]: 16 : if (errval != iBase_SUCCESS) throw Error(MK_FAILURE, "Unable get the facets of the surface.");
104 : :
105 : : //create vector for keeping track of the vertices
106 [ + - ]: 32 : std::vector<iBase_EntityHandle> verts;
107 : :
108 : : // loop over the facet points
109 [ + + ]: 4280 : for(unsigned int j=0; j<pnts.size(); j+=3)
110 : : {
111 : : //create a new vertex handle
112 : : iBase_EntityHandle v;
113 [ + - ][ + - ]: 4264 : mk->imesh_instance()->createVtx(pnts[j],pnts[j+1],pnts[j+2],v);
[ + - ][ + - ]
[ + - ]
114 [ + - ]: 4264 : verts.push_back(v);
115 : : }
116 : :
117 : :
118 : : // get the geometric vertices of the surface and merge them into the mesh by
119 : : // a proximity check
120 : :
121 : : //get the geometric vertices
122 [ + - ]: 32 : MEntVector geom_verts;
123 [ + - ]: 16 : surf->get_adjacencies(0,geom_verts);
124 : 16 : unsigned int matches = 0;
125 : :
126 : : //loop over all the vertices we've created
127 [ + + ]: 4280 : for(unsigned int j=0; j<verts.size() ; j++)
128 : : {
129 : : //loop over all geometric vertices adjacent to this surface
130 [ + + ]: 18512 : for(unsigned int k=0; k<geom_verts.size(); k++)
131 : : {
132 : : // if thhe current vert is within the geom_res proximity of a geometric vert
133 [ + - ][ + - ]: 14248 : if(vtx2vtx_dist(geom_verts[k]->geom_handle(),verts[j]) < geom_res)
[ + - ][ + - ]
[ + + ]
134 : : {
135 : : //replace the vertex with the geometric vertex
136 : :
137 : : //capture vert to remove from the list
138 [ + - ]: 60 : iMesh::EntityHandle vert_to_del = verts[j];
139 : :
140 : : //replace the vertex with the geom vertex in the vector
141 [ + - ]: 60 : std::vector<moab::EntityHandle> dum_handle;
142 [ + - ][ + - ]: 60 : geom_verts[k]->get_mesh(0, dum_handle);
143 [ + - ][ + - ]: 60 : verts[j] = IBEH(dum_handle[0]);
144 : :
145 : : //delete the former vertex
146 [ + - ][ + - ]: 60 : mk_core()->imesh_instance()->deleteEnt(vert_to_del);
[ + - ]
147 : :
148 : 60 : matches++;
149 : : }
150 : : }
151 : : }
152 : :
153 : : //now check that we have matched the correct number of vertices
154 [ - + ]: 16 : if (matches != geom_verts.size())
155 : : {
156 [ # # ][ # # ]: 0 : if(matches > geom_verts.size()) std::cout << "Warning: coincident vertices in surface " << surf->id() << std::endl;
[ # # ][ # # ]
[ # # ]
157 [ # # ][ # # ]: 0 : if(matches < geom_verts.size()) std::cout << "Warning: one or more geom vertices could not be matched in surface " << surf->id() << std::endl;
[ # # ][ # # ]
[ # # ]
158 : : }
159 : :
160 : : //vector for keeping track of the triangles
161 [ + - ]: 32 : std::vector<iBase_EntityHandle> tris;
162 : :
163 : : //loop over the connectivity
164 [ + + ]: 4248 : for(int j = 0; j < (int)conn.size()-2; j+=3)
165 : : {
166 : : //get the appropriate points for a triangle and add them to a vector
167 [ + - ]: 4232 : std::vector<iBase_EntityHandle> tri_verts;
168 [ + - ][ + - ]: 4232 : tri_verts.push_back(verts[conn[j]]);
[ + - ]
169 [ + - ][ + - ]: 4232 : tri_verts.push_back(verts[conn[j+1]]);
[ + - ]
170 [ + - ][ + - ]: 4232 : tri_verts.push_back(verts[conn[j+2]]);
[ + - ]
171 : :
172 : : //create a new triangle handle
173 : : iBase_EntityHandle t;
174 [ + - ][ + - ]: 4232 : mk->imesh_instance()->createEnt(iMesh_TRIANGLE, &tri_verts[0], 3, t);
[ + - ]
175 : 4232 : tri_verts.clear();
176 [ + - ]: 4232 : tris.push_back(t);
177 : 4232 : }
178 : : //std::cout << "Created " << tris.size() << " triangles" << std::endl;
179 : :
180 : : //add verticess and edges to the entity set
181 [ + - ][ + - ]: 16 : mk->imesh_instance()->addEntArrToSet(&verts[0], verts.size(), sh);
[ + - ]
182 [ + - ][ + - ]: 32 : mk->imesh_instance()->addEntArrToSet(&tris[0], tris.size(), sh);
[ + - ]
183 : :
184 : 16 : }
185 : :
186 : 2 : void SurfaceFacetMeshReader::set_mesh_params(double faceting_tolerance, double geom_resabs)
187 : : {
188 : : //Assign the faceting values if they are passed into the function
189 [ + - ]: 2 : if(faceting_tolerance) facet_tol = faceting_tolerance;
190 [ + - ]: 2 : if(geom_resabs) geom_res = geom_resabs;
191 : 2 : }
192 : :
193 : 14248 : double SurfaceFacetMeshReader::vtx2vtx_dist(iGeom::EntityHandle vtx1, iMesh::EntityHandle vtx2)
194 : : {
195 : :
196 : : double x1,y1,z1;
197 : : double x2,y2,z2;
198 : :
199 [ + - ][ + - ]: 14248 : mk_core()->igeom_instance()->getVtxCoord(vtx1, x1, y1, z1);
[ + - ]
200 [ + - ][ + - ]: 14248 : mk_core()->imesh_instance()->getVtxCoord(vtx2, x2, y2, z2);
[ + - ]
201 : :
202 [ + - ][ + - ]: 14248 : double dist = pow((x1-x2),2) + pow((y1-y2),2) + pow((z1-z2),2);
[ + - ]
203 : :
204 : 14248 : return sqrt(dist);
205 : : }
206 : :
207 [ + - ][ + - ]: 156 : }
|
