Branch data Line data Source code
1 : : #include "meshkit/AF2DfltTriangleMeshOp.hpp"
2 : :
3 : : // C++
4 : : #include <cstddef>
5 : : #include <list>
6 : : #include <vector>
7 : :
8 : : // MOAB
9 : : #include "moab/Interface.hpp"
10 : : #include "moab/Range.hpp"
11 : : #include "moab/ReadUtilIface.hpp"
12 : : #include "moab/Types.hpp"
13 : :
14 : : // MeshKit
15 : : #include "meshkit/AF2Algorithm.hpp"
16 : : #include "meshkit/AF2AlgorithmResult.hpp"
17 : : #include "meshkit/AF2DfltPlaneProjMaker.hpp"
18 : : #include "meshkit/AF2DfltTriangleRules.hpp"
19 : : #include "meshkit/AF2Point3D.hpp"
20 : : #include "meshkit/Error.hpp"
21 : : #include "meshkit/MeshScheme.hpp"
22 : : #include "meshkit/ModelEnt.hpp"
23 : : #include "meshkit/RegisterMeshOp.hpp"
24 : :
25 : : namespace MeshKit
26 : : {
27 : :
28 : : moab::EntityType AF2DfltTriangleMeshOp::meshTypes[] =
29 : : {moab::MBVERTEX, moab::MBTRI, moab::MBMAXTYPE};
30 : :
31 : 160 : bool AF2DfltTriangleMeshOp::can_mesh(iBase_EntityType dimension)
32 : : {
33 : 160 : return dimension == iBase_FACE;
34 : : }
35 : :
36 : 0 : bool AF2DfltTriangleMeshOp::can_mesh(ModelEnt *me)
37 : : {
38 : 0 : return canmesh_face(me);
39 : : }
40 : :
41 : 1056 : const char* AF2DfltTriangleMeshOp::name()
42 : : {
43 : 1056 : return "AF2DfltTriangleMeshOp";
44 : : }
45 : :
46 : 42 : const moab::EntityType* AF2DfltTriangleMeshOp::output_types()
47 : : {
48 : 42 : return meshTypes;
49 : : }
50 : :
51 : 9 : AF2DfltTriangleMeshOp::AF2DfltTriangleMeshOp(
52 : : MKCore *meshkitCore, const MEntVector &meshEntVec) :
53 : 9 : MeshScheme(meshkitCore, meshEntVec)
54 : : {
55 : : // don't do anything beyond constructing superclass with correct arguments
56 : 9 : }
57 : :
58 : 27 : AF2DfltTriangleMeshOp::~AF2DfltTriangleMeshOp()
59 : : {
60 : : // don't do anything; superclass destructor will be called
61 [ - + ]: 18 : }
62 : :
63 : 0 : AF2DfltTriangleMeshOp::AF2DfltTriangleMeshOp(
64 : 0 : const AF2DfltTriangleMeshOp& toCopy) : MeshScheme(toCopy)
65 : : {
66 : : // don't do anything beyond constructing superclass with correct arguments
67 : 0 : }
68 : :
69 : 0 : AF2DfltTriangleMeshOp& AF2DfltTriangleMeshOp::operator=(
70 : : const AF2DfltTriangleMeshOp& rhs)
71 : : {
72 [ # # ]: 0 : Error notImpl(MeshKit::MK_NOT_IMPLEMENTED);
73 : : notImpl.set_string(
74 [ # # ]: 0 : "AF2DfltTriangleMeshOp assignment operator is not supported.");
75 [ # # ]: 0 : throw notImpl;
76 : : }
77 : :
78 : 9 : void AF2DfltTriangleMeshOp::execute_this()
79 : : {
80 [ + - ]: 9 : AF2DfltTriangleRules defaultRules;
81 [ + - ][ + - ]: 18 : AF2Algorithm algorithm(defaultRules.getRules());
82 : :
83 [ + - ]: 9 : int debug = get_debug_verbosity(); // will call MeshOp to get the debug verbosity level
84 : :
85 : : // if some debug level, dump out the mesh so far
86 [ - + ]: 9 : if (debug>0)
87 : : {
88 : : // get the moab instance for convenience
89 [ # # ][ # # ]: 0 : moab::Interface* moabInstance = mk_core()->moab_instance();
90 [ # # ]: 0 : moabInstance->write_file("meshBeforeAF2.vtk");
91 : : }
92 [ + - + - ]: 78 : for (MEntSelection::iterator selIt = mentSelection.begin();
[ + + ]
93 : 39 : selIt != mentSelection.end(); ++selIt)
94 : : {
95 : : // extract the model entity from the map iterator
96 [ + - ]: 30 : ModelEnt* modelEnt = selIt->first;
97 : :
98 [ + - ][ - + ]: 30 : if (modelEnt->get_meshed_state() >= COMPLETE_MESH)
99 : : {
100 : 0 : continue;
101 : : }
102 : :
103 : : // get the moab instance for convenience
104 [ + - ][ + - ]: 30 : moab::Interface* moabInstance = mk_core()->moab_instance();
105 : :
106 : : // query the boundary to find mesh edges
107 : : // Note: the senses returned by this method should be either FORWARD
108 : : // or REVERSE, not BOTH . . .
109 : : // but I think there is a chance that mesh edges that are on a
110 : : // geometric edge that has a sense of BOTH might not agree with
111 : : // nearby edges and thus might not form a wire when placed end
112 : : // to end. The AF2Algorithm doesn't require that the input be in
113 : : // wire order, and should be okay as long as the edge actually
114 : : // does appear with both senses.
115 [ + - ]: 30 : std::vector<moab::EntityHandle> bndryEdges;
116 [ + - ]: 60 : std::vector<int> bEdgeSenses;
117 [ + - ]: 30 : modelEnt->boundary(1, bndryEdges, &bEdgeSenses);
118 : :
119 : : // get lists of the vertex coordinates and handles
120 : : // along with an index from the edge handles into those lists
121 [ + - ]: 60 : std::vector<int> bndryIds;
122 [ + - ]: 60 : std::vector<double> coords;
123 [ + - ]: 60 : moab::Range bndryVertsRange;
124 : : modelEnt->get_indexed_connect_coords(
125 [ + - ]: 30 : bndryEdges, &bEdgeSenses, NULL, bndryIds, coords, &bndryVertsRange);
126 : : std::vector<moab::EntityHandle> bndryVertsVec(
127 [ + - ][ + - ]: 60 : bndryVertsRange.begin(), bndryVertsRange.end());
[ + - ]
128 : :
129 : : // convert the edges into the unsigned int format that is needed for
130 : : // input into the AF2Algorithm
131 : : // Note: The tags used inside get_indexed_connect_coords are actually
132 : : // unsigned int (as of May 2016)
133 [ + - ]: 60 : std::vector<unsigned int> inputEdges;
134 [ + - ]: 4628 : for (std::vector<int>::const_iterator bidItr = bndryIds.begin();
[ + - + - ]
[ + + ]
135 : 2314 : bidItr != bndryIds.end(); ++bidItr)
136 : : {
137 [ + - ][ + - ]: 2284 : inputEdges.push_back(*bidItr);
138 : : }
139 : :
140 : : // determine the number of boundary edges and vertices
141 : 30 : unsigned int numVertices = bndryVertsVec.size();
142 : 30 : unsigned int numEdges = inputEdges.size() / 2;
143 : :
144 : : // check that there is reasonable data to initialize the front
145 [ - + ]: 30 : if (numEdges <= 0)
146 : : {
147 [ # # ]: 0 : Error failErr(MeshKit::MK_FAILURE);
148 [ # # ]: 0 : failErr.set_string("There are no boundary mesh edges to use to initialize the advancing front in AF2DfltTriangleMeshOp.");
149 [ # # ]: 0 : throw failErr;
150 : : }
151 : : // This should never happen if there are edges, but as a sanity check
152 : : // since the coordinates of the first vertex may be used to check sizing
153 [ - + ]: 30 : if (numVertices <= 0)
154 : : {
155 [ # # ]: 0 : Error failErr(MeshKit::MK_FAILURE);
156 : : failErr.set_string(
157 [ # # ]: 0 : "No boundary vertices found in AF2DfltTriangleMeshOp.");
158 [ # # ]: 0 : throw failErr;
159 : : }
160 : :
161 : : // configure the sizing function, if any, that will be used
162 : 30 : SizingFunction* sizing = NULL;
163 [ + - ][ + + ]: 30 : if (modelEnt->sizing_function_index() > -1)
164 : : {
165 [ + - ]: 29 : int sfIndex = modelEnt->sizing_function_index();
166 [ + - ][ + - ]: 29 : SizingFunction* tempSizing = mk_core()->sizing_function(sfIndex);
167 : : // check the desired size at the coordinates of the first vertex
168 [ + - ][ + - ]: 29 : if (tempSizing->size(&coords[0]) > 0)
[ + - ]
169 : : {
170 : 29 : sizing = tempSizing;
171 : : }
172 : : }
173 : :
174 : : // construct the local transform maker and run the algorithm
175 : : AF2DfltPlaneProjMaker localTransformMaker(
176 [ + - ][ + - ]: 60 : modelEnt->igeom_instance(), modelEnt->geom_handle(), sizing);
[ + - ]
177 : :
178 : : AF2AlgorithmResult* meshResult = algorithm.execute(&localTransformMaker,
179 [ + - ][ + - ]: 30 : &coords[0], numVertices, &inputEdges[0], numEdges, &bndryVertsVec[0], debug);
[ + - ][ + - ]
180 : :
181 : : // throw failure if the algorithm did not succeed
182 [ + - ][ - + ]: 30 : if (!meshResult->isSuccessful())
183 : : {
184 [ # # ]: 0 : delete meshResult;
185 [ # # ]: 0 : Error failErr(MeshKit::MK_FAILURE);
186 [ # # ]: 0 : failErr.set_string("AF2DfltTriangleMeshOp failed.");
187 [ # # ]: 0 : throw failErr;
188 : : }
189 : :
190 : : // Collect the new points in a vector
191 [ + - ]: 30 : const std::list<AF2Point3D*>* pointList = meshResult->getPoints();
192 [ + - ]: 60 : std::vector<double> newPointCoords;
193 : : typedef std::list<AF2Point3D*>::const_iterator ConstPointItr;
194 : 30 : int numNewPoints = 0;
195 [ + - + - ]: 10266 : for (ConstPointItr pItr = pointList->begin();
[ + + ]
196 : 5133 : pItr != pointList->end(); ++pItr)
197 : : {
198 [ + - ][ + - ]: 5103 : if (!(*pItr)->isCommitted())
[ + + ]
199 : : {
200 : 3961 : ++numNewPoints;
201 [ + - ][ + - ]: 3961 : newPointCoords.push_back((*pItr)->getX());
[ + - ]
202 [ + - ][ + - ]: 3961 : newPointCoords.push_back((*pItr)->getY());
[ + - ]
203 [ + - ][ + - ]: 3961 : newPointCoords.push_back((*pItr)->getZ());
[ + - ]
204 : : }
205 : : }
206 : : // Commit the new points to MOAB
207 [ + - ]: 60 : moab::Range newPointsRange;
208 : : moabInstance->create_vertices(
209 [ + - ][ + - ]: 30 : &newPointCoords[0], numNewPoints, newPointsRange);
210 : : // Set the MOAB handles for the now-committed vertices
211 : 30 : ConstPointItr np3dItr = pointList->begin();
212 [ + - ][ + - ]: 7982 : for (moab::Range::const_iterator nphItr = newPointsRange.begin();
[ + - ][ + + ]
213 [ + - ]: 3991 : nphItr != newPointsRange.end(); ++nphItr)
214 : : {
215 [ + - ][ + - ]: 5103 : while ((*np3dItr)->isCommitted())
[ + + ]
216 : : {
217 [ + - ]: 1142 : ++np3dItr;
218 : : }
219 [ + - ][ + - ]: 3961 : (*np3dItr)->setCommittedHandle(*nphItr);
[ + - ]
220 [ + - ]: 3961 : ++np3dItr;
221 : : }
222 : :
223 : : // get a pointer to the list of new triangles and check how many there are
224 [ + - ]: 30 : const std::list<const AF2Polygon3D*>* faceList = meshResult->getFaces();
225 : 30 : int numTriangles = faceList->size();
226 : :
227 : : // pre-allocate connectivity memory to store the triangles
228 : : moab::ReadUtilIface* readInterface;
229 : : moab::ErrorCode moabRet;
230 [ + - ]: 30 : moabRet = moabInstance->query_interface(readInterface);
231 [ - + ][ # # ]: 30 : MBERRCHK(moabRet, moabInstance);
[ # # ][ # # ]
232 : : moab::EntityHandle firstHandle;
233 : : moab::EntityHandle* triConnectAry;
234 : : moabRet = readInterface->get_element_connect(
235 [ + - ]: 30 : numTriangles, 3, moab::MBTRI, 0, firstHandle, triConnectAry);
236 [ - + ][ # # ]: 30 : MBERRCHK(moabRet, moabInstance);
[ # # ][ # # ]
237 : :
238 : : // fill the connectivity array
239 : 30 : unsigned int caIndex = 0u;
240 : : typedef std::list<const AF2Polygon3D*>::const_iterator ConstTriangleItr;
241 [ + - + - ]: 18116 : for (ConstTriangleItr tItr = faceList->begin();
[ + + ]
242 : 9058 : tItr != faceList->end(); ++tItr)
243 : : {
244 [ + + ]: 36112 : for (unsigned int vi = 0; vi < 3; ++vi)
245 : : {
246 [ + - ][ + - ]: 27084 : triConnectAry[caIndex + vi] = (*tItr)->getVertex(vi)->getVertexHandle();
[ + - ]
247 : : }
248 : 9028 : caIndex += 3;
249 : : }
250 : :
251 [ + - ]: 30 : delete meshResult;
252 : :
253 : : // Insert the new entities into the entity range
254 [ + - ][ + - ]: 30 : selIt->second.insert(firstHandle, firstHandle + numTriangles - 1);
255 : :
256 : : // Commit the mesh
257 [ + - ][ + - ]: 30 : modelEnt->commit_mesh(selIt->second, COMPLETE_MESH);
258 : 39 : }
259 : 9 : }
260 : :
261 : 0 : const moab::EntityType* AF2DfltTriangleMeshOp::mesh_types_arr() const
262 : : {
263 : 0 : return output_types();
264 : : }
265 : :
266 : 9 : void AF2DfltTriangleMeshOp::setup_this()
267 : : {
268 : : // Check that the dimension is correct for each model entity and pass
269 : : // down the sizing functions that are set on the selected surfaces
270 : : // to their boundary edges if those edges don't have sizing functions set
271 [ + - + - ]: 78 : for (MEntSelection::iterator selIt = mentSelection.begin();
[ + + ]
272 : 39 : selIt != mentSelection.end(); ++selIt)
273 : : {
274 : : // extract the model entity from the map iterator
275 [ + - ]: 30 : ModelEnt* modelEnt = selIt->first;
276 : :
277 : : // check that the dimension of the selected model entity is two
278 [ + - ][ - + ]: 30 : if (modelEnt->dimension() != 2)
279 : : {
280 [ # # ]: 0 : Error dimErr(MeshKit::MK_WRONG_DIMENSION);
281 [ # # ]: 0 : dimErr.set_string("Found a selected entity of dimension not equal to 2 in AF2DfltTriangleMeshOp.");
282 [ # # ]: 0 : throw dimErr;
283 : : }
284 : :
285 : : // Try to check whether there is a valid sizing function set on the
286 : : // model entity.
287 : 30 : int sFIndex = -1;
288 [ + - ][ + + ]: 30 : if (modelEnt->sizing_function_index() > -1)
289 : : {
290 [ + - ]: 29 : int tempSFIndex = modelEnt->sizing_function_index();
291 [ + - ][ + - ]: 29 : SizingFunction* tempSizing = mk_core()->sizing_function(tempSFIndex);
292 : : // check the desired size at the coordinates of the first vertex
293 [ + - ][ + - ]: 29 : if (tempSizing->size() > 0)
294 : : {
295 : 29 : sFIndex = tempSFIndex;
296 : : }
297 : : }
298 : :
299 : : // If the sizing function appears to be valid, pass it down to adjacent
300 : : // children that do not have a sizing function set, because it is best
301 : : // to have the size on the boundaries about right.
302 [ + + ]: 30 : if (sFIndex > -1)
303 : : {
304 [ + - ]: 29 : MEntVector bndryEdges;
305 [ + - ]: 29 : modelEnt->get_adjacencies(1, bndryEdges);
306 [ + - + - ]: 358 : for (MEntVector::iterator beItr = bndryEdges.begin();
[ + + ]
307 : 179 : beItr != bndryEdges.end(); ++beItr)
308 : : {
309 [ + - ][ + - ]: 150 : int beSFIndex = (*beItr)->sizing_function_index();
310 [ + + ]: 150 : if (beSFIndex == -1)
311 : : {
312 [ + - ][ + - ]: 2 : (*beItr)->sizing_function_index(sFIndex);
313 : : }
314 : 29 : }
315 : : }
316 : : }
317 : :
318 : : // Set up the default MeshOp on any children that need one
319 : 9 : setup_boundary();
320 : :
321 : : // Make sure that this MeshOp depends on all MeshOps set on its facets
322 : 9 : ensure_facet_dependencies(false);
323 : 9 : }
324 : :
325 [ + - ][ + - ]: 156 : } // namespace MeshKit
|
