Branch data Line data Source code
1 : : #include "meshkit/AF2Algorithm.hpp"
2 : :
3 : : // C++
4 : : #include <cstddef>
5 : : #include <vector>
6 : : #include <map>
7 : :
8 : : // MeshKit
9 : : #include "meshkit/AF2DfltRuleAppVisitor.hpp"
10 : : #include "meshkit/AF2Edge3D.hpp"
11 : : #include "meshkit/AF2RuleApplication.hpp"
12 : : #include "meshkit/Error.hpp"
13 : :
14 : : // for debugging
15 : : #include <iostream>
16 : :
17 : 14 : AF2Algorithm::AF2Algorithm(const std::list<const AF2Rule*> & ruleListArg) :
18 : 14 : ruleList(ruleListArg)
19 : : {
20 : : // do nothing beyond copying the list of rules in initialization
21 : 14 : }
22 : :
23 : 28 : AF2Algorithm::~AF2Algorithm()
24 : : {
25 : : // do nothing, the standard deletion of the rule list is sufficient,
26 : : // since this does not take ownership of the rules
27 : 14 : }
28 : :
29 : 0 : AF2Algorithm::AF2Algorithm(const AF2Algorithm & toCopy)
30 : : {
31 : : // Note: the default implementation would work at this point, but
32 : : // might not in the future, and there shouldn't be much need to
33 : : // copy this.
34 [ # # ]: 0 : MeshKit::Error notImpl(MeshKit::MK_NOT_IMPLEMENTED);
35 [ # # ]: 0 : notImpl.set_string("AF2Algorithm copy construction is not supported.");
36 [ # # ]: 0 : throw notImpl;
37 : : }
38 : :
39 : 0 : AF2Algorithm& AF2Algorithm::operator=(const AF2Algorithm & rhs)
40 : : {
41 : : // Note: the default implementation would work at this point, but
42 : : // might not in the future, and there shouldn't be much need to
43 : : // assign this.
44 [ # # ]: 0 : MeshKit::Error notImpl(MeshKit::MK_NOT_IMPLEMENTED);
45 [ # # ]: 0 : notImpl.set_string("AF2Algorithm assignment operator is not supported.");
46 [ # # ]: 0 : throw notImpl;
47 : : }
48 : :
49 : 35 : AF2AlgorithmResult* AF2Algorithm::execute(
50 : : const AF2LocalTransformMaker* const & transformMaker,
51 : : const double* coords, unsigned int numPoints,
52 : : const unsigned int* edges, unsigned int numEdges,
53 : : const moab::EntityHandle* vertexHandles, int debug) const
54 : : {
55 : : typedef std::list<const AF2Rule*>::const_iterator RuleConstItr;
56 : :
57 : : // these are used for debugging
58 : : static int face=0;
59 : 35 : face++;
60 : :
61 : 35 : unsigned long int nextPointId = 0ul;
62 [ + - ]: 35 : std::list<AF2Point3D*> allPoints;
63 [ + - ]: 70 : std::list<const AF2Polygon3D*> allFaces;
64 [ + - ]: 70 : AF2Front front;
65 : :
66 : : // initialize the front
67 : : initFront(front, allPoints, nextPointId,
68 [ + - ]: 35 : coords, numPoints, edges, numEdges, vertexHandles);
69 [ - + ]: 35 : if (debug>1)
70 : : {
71 [ # # ][ # # ]: 0 : std::cout<<" On surface " << face << " initial front " << allPoints.size() << " points and " << numEdges << " edges.\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
72 : : }
73 : 35 : int step = 0;
74 : : // while the front is not empty and there is still progress
75 [ + - ][ + + ]: 10016 : while (!front.isEmpty() && front.getMaximumQuality() < 50u)
[ + - ][ + + ]
[ + + ]
76 : : {
77 : : // select a neighborhood on the advancing front
78 [ + - ]: 9981 : AF2Neighborhood* ngbhd = front.selectNeighborhood(transformMaker);
79 [ + - ]: 9981 : AF2Edge3D* baselineEdge = ngbhd->getBaselineEdge3D();
80 : :
81 : : // attempt to apply each of the rules to the neighborhood
82 [ + - ]: 9981 : AF2DfltRuleAppVisitor ruleAppVisitor;
83 [ + - ][ + - ]: 127413 : for (RuleConstItr itr = ruleList.begin(); itr != ruleList.end(); ++itr)
[ + + ]
84 : : {
85 [ + - ]: 117432 : (*itr)->applyRule(*ngbhd,
86 [ + - ][ + - ]: 234864 : baselineEdge->getQualityLevel(), ruleAppVisitor);
87 : : }
88 : :
89 : : // process the results of attempting to apply rules
90 : : const AF2RuleApplication* bestRuleApp =
91 [ + - ]: 9981 : ruleAppVisitor.getBestRuleApplication();
92 [ + + ]: 9981 : if (bestRuleApp == NULL)
93 : : {
94 : : // there were no successful rule applications, so decrease the quality
95 [ + - ]: 3438 : baselineEdge->decreaseQuality();
96 : : }
97 : : else
98 : : {
99 : : // there was a successful rule application, so process it and advance
100 : :
101 : : // add any new points added by the rule application
102 [ + - ]: 6543 : std::map<const AF2Point2D*, AF2Point3D*> newPointsMap;
103 [ + - ][ + + ]: 10594 : for (unsigned int npi = 0; npi < bestRuleApp->getNumNewPoints(); ++npi)
104 : : {
105 : : processNewPoint(bestRuleApp->getNewPoint(npi), nextPointId,
106 [ + - ][ + - ]: 4051 : ngbhd, newPointsMap, allPoints, front);
107 : : }
108 : :
109 : : // add the faces added by the rule application
110 [ + - ][ + + ]: 15807 : for (unsigned int nfi = 0; nfi < bestRuleApp->getNumNewFaces(); ++nfi)
111 : : {
112 : : processNewFace(bestRuleApp->getNewFace(nfi),
113 [ + - ][ + - ]: 9264 : ngbhd, newPointsMap, allFaces, front);
114 : : }
115 : : // at a successful application, dump the content if we want to
116 : 6543 : step++;
117 [ - + ]: 6543 : if (debug>2)
118 : : {
119 [ # # ]: 0 : output_intermediate_result(allPoints, allFaces, face, step);
120 : 6543 : }
121 : : }
122 : 9981 : }
123 : :
124 [ + - ][ + + ]: 35 : if (front.isEmpty())
125 : : {
126 : : // the advancing front algorithm successfully completed
127 [ - + ]: 34 : if (debug>0)
128 : : {
129 [ # # ][ # # ]: 0 : std::cout<<" On surface " << face << " generated " << allPoints.size() << " points and " << allFaces.size() << " triangles\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
130 : : }
131 [ + - ][ + - ]: 34 : return new AF2AlgorithmResult(allPoints, allFaces);
132 : : }
133 : :
134 : : // the advancing front algorithm failed
135 [ + - ]: 1 : release(allPoints, allFaces);
136 [ + - ][ + - ]: 36 : return new AF2AlgorithmResult();
137 : : }
138 : :
139 : 35 : void AF2Algorithm::initFront(AF2Front & front, std::list<AF2Point3D*> & pntList,
140 : : unsigned long & pntId,
141 : : const double* coords, unsigned int numPoints,
142 : : const unsigned int* edges, unsigned int numEdges,
143 : : const moab::EntityHandle* vertexHandles) const
144 : : {
145 : : // make the point objects and add them to the front
146 [ + - ]: 35 : std::vector<AF2Point3D*> pntVector;
147 [ + - ]: 35 : pntVector.reserve(numPoints);
148 [ + + ]: 1243 : for (unsigned int pi = 0; pi < numPoints; ++pi)
149 : : {
150 : : AF2Point3D* point = new AF2Point3D(pntId,
151 [ + - ][ + - ]: 1208 : coords[3*pi], coords[3*pi + 1], coords[3*pi + 2]);
152 : 1208 : ++pntId;
153 [ + + ]: 1208 : if (vertexHandles != NULL)
154 : : {
155 [ + - ]: 1142 : point->setCommittedHandle(vertexHandles[pi]);
156 : : }
157 [ + - ]: 1208 : pntVector.push_back(point);
158 [ + - ]: 1208 : pntList.push_back(point);
159 [ + - ]: 1208 : front.addPoint(point);
160 : : }
161 : :
162 : : // make the edge objects and gather them in a list of edges
163 [ + - ]: 70 : std::list<AF2Edge3D*> edgeList;
164 [ + + ]: 1243 : for (unsigned int ei = 0; ei < numEdges; ++ei)
165 : : {
166 [ + - ][ - + ]: 1208 : if (edges[2*ei] >= numPoints || edges[2*ei + 1] >= numPoints)
167 : : {
168 [ # # ]: 0 : MeshKit::Error badArg(MeshKit::MK_BAD_INPUT);
169 [ # # ]: 0 : badArg.set_string("An edge index exceeds the number of points.");
170 [ # # ]: 0 : throw badArg;
171 : : }
172 [ + - ]: 1208 : AF2Point3D* edgeStart = pntVector[edges[2*ei]];
173 [ + - ]: 1208 : AF2Point3D* edgeEnd = pntVector[edges[2*ei + 1]];
174 : : // mark the endpoints of the edges as on the initial boundary
175 : : // Note: This is not done for isolated points, since the front
176 : : // does not start to advance from them.
177 [ + - ]: 1208 : edgeStart->limitDistanceToBoundary(0);
178 [ + - ]: 1208 : edgeEnd->limitDistanceToBoundary(0);
179 [ + - ][ + - ]: 1208 : AF2Edge3D* edge = new AF2Edge3D(edgeStart, edgeEnd);
180 [ + - ]: 1208 : edgeList.push_back(edge);
181 : : }
182 : :
183 : : // initialize the front with the edges
184 [ + - ][ + - ]: 70 : front.advanceFront(edgeList);
185 : 35 : }
186 : :
187 : 9264 : void AF2Algorithm::processNewFace(const AF2Polygon2D* newFace2D,
188 : : AF2Neighborhood* & ngbhd,
189 : : std::map<const AF2Point2D*, AF2Point3D*> & newPointsMap,
190 : : std::list<const AF2Polygon3D*> & allFaces, AF2Front & front) const
191 : : {
192 [ + - ]: 9264 : std::list<const AF2Point3D*> facePoints3D;
193 [ + - ]: 18528 : std::list<AF2Edge3D*> edgeList;
194 : :
195 : 9264 : AF2Point3D* firstVertex = NULL;
196 : 9264 : AF2Point3D* prevVertex = NULL;
197 : 9264 : AF2Point3D* curVertex = NULL;
198 [ + - ][ + + ]: 37056 : for (unsigned int fvi = 0; fvi < newFace2D->getNumVertices(); ++fvi)
199 : : {
200 : : // get the 3-D point corresponding to the 2-D point . . .
201 [ + - ]: 27792 : const AF2Point2D* faceVertex2D = newFace2D->getVertex(fvi);
202 : : // . . . either from the neighborhood . . .
203 [ + - ]: 27792 : curVertex = ngbhd->getCorrespondingPoint(faceVertex2D);
204 : : // . . . or from the map of new points that were added.
205 [ + + ]: 27792 : if (curVertex == NULL)
206 : : {
207 [ + - ]: 6772 : curVertex = newPointsMap[faceVertex2D];
208 : : }
209 : :
210 : : // add it to the list of points for future creation of the 3-D face
211 [ + - ]: 27792 : facePoints3D.push_back(curVertex);
212 : :
213 : : // process the edge
214 [ + + ]: 27792 : if (firstVertex == NULL)
215 : : {
216 : : // at the first vertex just store the first vertex and what will be
217 : : // the previous vertex
218 : 9264 : firstVertex = curVertex;
219 : 9264 : prevVertex = curVertex;
220 : : }
221 : : else
222 : : {
223 : : // create and store an edge from the current to the previous vertex
224 [ + - ][ + - ]: 18528 : edgeList.push_back(new AF2Edge3D(curVertex, prevVertex));
[ + - ]
225 : : // update what will be the previous vertex
226 : 18528 : prevVertex = curVertex;
227 : : }
228 : : }
229 : : // create and store an edge from the first vertex to the last vertex
230 [ + - ][ + - ]: 9264 : edgeList.push_back(new AF2Edge3D(firstVertex, curVertex));
[ + - ]
231 : :
232 : : // create the three-dimensional face and add it to the list of faces
233 [ + - ][ + - ]: 9264 : allFaces.push_back(new AF2Polygon3D(facePoints3D));
[ + - ]
234 : :
235 : : // advance the front with the edges
236 [ + - ][ + - ]: 18528 : front.advanceFront(edgeList);
237 : 9264 : }
238 : :
239 : 4051 : void AF2Algorithm::processNewPoint(const AF2Point2D* newPoint2D,
240 : : unsigned long & pntId,
241 : : AF2Neighborhood* & ngbhd,
242 : : std::map<const AF2Point2D*, AF2Point3D*> & newPointsMap,
243 : : std::list<AF2Point3D*> & allPoints, AF2Front & front) const
244 : : {
245 [ + - ]: 4051 : AF2Point3D* newPoint3D = ngbhd->transformPoint(newPoint2D, pntId);
246 : 4051 : ++pntId;
247 [ + - ]: 4051 : newPointsMap[newPoint2D] = newPoint3D;
248 [ + - ]: 4051 : allPoints.push_back(newPoint3D);
249 [ + - ]: 4051 : front.addPoint(newPoint3D);
250 : 4051 : }
251 : :
252 : 1 : void AF2Algorithm::release(std::list<AF2Point3D*> & allPoints,
253 : : std::list<const AF2Polygon3D*> & allFaces) const
254 : : {
255 : : typedef std::list<const AF2Polygon3D*>::iterator FaceItr;
256 : : typedef std::list<AF2Point3D*>::iterator PointItr;
257 : :
258 [ + - ][ + - ]: 3 : for (FaceItr itr = allFaces.begin(); itr != allFaces.end(); ++itr)
[ + + ]
259 : : {
260 [ + - ][ + - ]: 2 : delete (*itr);
261 : : }
262 : :
263 [ + - ][ + - ]: 7 : for (PointItr itr = allPoints.begin(); itr != allPoints.end(); ++itr)
[ + + ]
264 : : {
265 [ + - ]: 6 : delete (*itr);
266 : : }
267 : 1 : }
268 : :
269 : : #include "moab/Core.hpp"
270 : : #include "moab/ReadUtilIface.hpp"
271 : 0 : void AF2Algorithm::output_intermediate_result (std::list<AF2Point3D*> & allPoints,
272 : : std::list<const AF2Polygon3D*> & allFaces,int face, int step) const
273 : : {
274 [ # # ]: 0 : moab::Core mb;
275 : 0 : int num_nodes = (int) allPoints.size();
276 [ # # ]: 0 : std::vector<double> newPointCoords;
277 : : typedef std::list<AF2Point3D*>::const_iterator ConstPointItr;
278 [ # # ][ # # ]: 0 : for (ConstPointItr pItr = allPoints.begin();
[ # # ][ # # ]
279 [ # # ]: 0 : pItr != allPoints.end(); ++pItr)
280 : : {
281 [ # # ][ # # ]: 0 : newPointCoords.push_back((*pItr)->getX());
[ # # ]
282 [ # # ][ # # ]: 0 : newPointCoords.push_back((*pItr)->getY());
[ # # ]
283 [ # # ][ # # ]: 0 : newPointCoords.push_back((*pItr)->getZ());
[ # # ]
284 : : }
285 : : // Commit the new points to MOAB
286 [ # # ]: 0 : moab::Range newPointsRange;
287 : : mb.create_vertices(
288 [ # # ][ # # ]: 0 : &newPointCoords[0], num_nodes, newPointsRange);
289 : : // Set the map between local ID and moab handle
290 [ # # ]: 0 : std::map<long, moab::EntityHandle> idToHandle;
291 [ # # ]: 0 : ConstPointItr np3dItr = allPoints.begin();
292 [ # # ][ # # ]: 0 : for (moab::Range::const_iterator nphItr = newPointsRange.begin();
[ # # ][ # # ]
293 [ # # ]: 0 : nphItr != newPointsRange.end(); ++nphItr)
294 : : {
295 : :
296 [ # # ][ # # ]: 0 : long localID = (*np3dItr)->getLocalId();
297 [ # # ][ # # ]: 0 : idToHandle[localID]=*nphItr;
298 [ # # ]: 0 : ++np3dItr;
299 : : }
300 : :
301 : 0 : int numTriangles = allFaces.size();
302 : :
303 : : // pre-allocate connectivity memory to store the triangles
304 : : moab::ReadUtilIface* readInterface;
305 [ # # ]: 0 : mb.query_interface(readInterface);
306 : :
307 : : moab::EntityHandle firstHandle;
308 : : moab::EntityHandle* triConnectAry;
309 : : readInterface->get_element_connect(
310 [ # # ]: 0 : numTriangles, 3, moab::MBTRI, 0, firstHandle, triConnectAry);
311 : :
312 : 0 : unsigned int caIndex = 0u;
313 : : typedef std::list<const AF2Polygon3D*>::const_iterator ConstTriangleItr;
314 [ # # ][ # # ]: 0 : for (ConstTriangleItr tItr = allFaces.begin();
[ # # ][ # # ]
315 [ # # ]: 0 : tItr != allFaces.end(); ++tItr)
316 : : {
317 [ # # ]: 0 : for (unsigned int vi = 0; vi < 3; ++vi)
318 : : {
319 [ # # ][ # # ]: 0 : triConnectAry[caIndex + vi] = idToHandle[(*tItr)->getVertex(vi)->getLocalId()];
[ # # ][ # # ]
320 : : }
321 : 0 : caIndex += 3;
322 : : }
323 : :
324 [ # # ][ # # ]: 0 : std::stringstream tempStep;
325 [ # # ][ # # ]: 0 : tempStep<<"Face_"<< face << "_Step" <<step<< ".vtk";
[ # # ][ # # ]
[ # # ]
326 [ # # ][ # # ]: 0 : mb.write_file(tempStep.str().c_str());
327 : :
328 : :
329 : 0 : return;
330 [ + - ][ + - ]: 156 : }
331 : :
|
