Branch data Line data Source code
1 : : /*
2 : : * Intx2MeshInPlane.cpp
3 : : *
4 : : * Created on: Oct 24, 2012
5 : : * Author: iulian
6 : : */
7 : :
8 : : #include "moab/IntxMesh/Intx2MeshInPlane.hpp"
9 : : #include "moab/GeomUtil.hpp"
10 : : #include "moab/IntxMesh/IntxUtils.hpp"
11 : :
12 : : namespace moab
13 : : {
14 : :
15 : 2 : Intx2MeshInPlane::Intx2MeshInPlane( Interface* mbimpl ) : Intx2Mesh( mbimpl ) {}
16 : :
17 [ - + ]: 2 : Intx2MeshInPlane::~Intx2MeshInPlane() {}
18 : :
19 : 96 : double Intx2MeshInPlane::setup_tgt_cell( EntityHandle tgt, int& nsTgt )
20 : : {
21 : : // the points will be at most ?; they will describe a convex patch, after the points will be
22 : : // ordered and collapsed (eliminate doubles) the area is not really required get coordinates of
23 : : // the tgt quad
24 : 96 : double cellArea = 0;
25 : : int num_nodes;
26 [ + - ]: 96 : ErrorCode rval = mb->get_connectivity( tgt, tgtConn, num_nodes );
27 [ - + ]: 96 : if( MB_SUCCESS != rval ) return 1.; // it should be an error
28 : :
29 : 96 : nsTgt = num_nodes;
30 : :
31 [ + - ][ + - ]: 96 : rval = mb->get_coords( tgtConn, num_nodes, &( tgtCoords[0][0] ) );
32 [ - + ]: 96 : if( MB_SUCCESS != rval ) return 1.; // it should be an error
33 : :
34 [ + + ]: 384 : for( int j = 0; j < nsTgt; j++ )
35 : : {
36 : : // populate coords in the plane for intersection
37 : : // they should be oriented correctly, positively
38 [ + - ]: 288 : tgtCoords2D[2 * j] = tgtCoords[j][0]; // x coordinate,
39 [ + - ]: 288 : tgtCoords2D[2 * j + 1] = tgtCoords[j][1]; // y coordinate
40 : : }
41 [ + + ]: 192 : for( int j = 1; j < nsTgt - 1; j++ )
42 [ + - ]: 96 : cellArea += IntxUtils::area2D( &tgtCoords2D[0], &tgtCoords2D[2 * j], &tgtCoords2D[2 * j + 2] );
43 : :
44 : 96 : return cellArea;
45 : : }
46 : :
47 : 250 : ErrorCode Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc( EntityHandle tgt, EntityHandle src, double* P, int& nP,
48 : : double& area, int markb[MAXEDGES], int markr[MAXEDGES],
49 : : int& nsSrc, int& nsTgt, bool check_boxes_first )
50 : : {
51 : :
52 : 250 : int num_nodes = 0;
53 [ + - ][ - + ]: 250 : ErrorCode rval = mb->get_connectivity( src, srcConn, num_nodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
54 : :
55 : 250 : nsSrc = num_nodes;
56 [ + - ][ + - ]: 250 : rval = mb->get_coords( srcConn, num_nodes, &( srcCoords[0][0] ) );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
57 : :
58 : 250 : area = 0.;
59 : 250 : nP = 0; // number of intersection points we are marking the boundary of src!
60 [ + + ]: 250 : if( check_boxes_first )
61 : : {
62 [ + - ]: 1 : setup_tgt_cell( tgt, nsTgt ); // we do not need area here
63 : : // look at the boxes formed with vertices; if they are far away, return false early
64 [ + - ][ - + ]: 1 : if( !GeomUtil::bounding_boxes_overlap( tgtCoords, nsTgt, srcCoords, nsSrc, box_error ) )
65 : 0 : return MB_SUCCESS; // no error, but no intersection, decide early to get out
66 : : }
67 : : #ifdef ENABLE_DEBUG
68 : : if( dbg_1 )
69 : : {
70 : : std::cout << "tgt " << mb->id_from_handle( tgt ) << "\n";
71 : : for( int j = 0; j < nsTgt; j++ )
72 : : {
73 : : std::cout << tgtCoords[j] << "\n";
74 : : }
75 : : std::cout << "src " << mb->id_from_handle( src ) << "\n";
76 : : for( int j = 0; j < nsSrc; j++ )
77 : : {
78 : : std::cout << srcCoords[j] << "\n";
79 : : }
80 : : mb->list_entities( &tgt, 1 );
81 : : mb->list_entities( &src, 1 );
82 : : }
83 : : #endif
84 : :
85 [ + + ]: 1000 : for( int j = 0; j < nsSrc; j++ )
86 : : {
87 [ + - ]: 750 : srcCoords2D[2 * j] = srcCoords[j][0]; // x coordinate,
88 [ + - ]: 750 : srcCoords2D[2 * j + 1] = srcCoords[j][1]; // y coordinate
89 : : }
90 : : #ifdef ENABLE_DEBUG
91 : : if( dbg_1 )
92 : : {
93 : : // std::cout << "gnomonic plane: " << plane << "\n";
94 : : std::cout << " tgt \n";
95 : : for( int j = 0; j < nsTgt; j++ )
96 : : {
97 : : std::cout << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << "\n ";
98 : : }
99 : : std::cout << " src\n";
100 : : for( int j = 0; j < nsSrc; j++ )
101 : : {
102 : : std::cout << srcCoords2D[2 * j] << " " << srcCoords2D[2 * j + 1] << "\n";
103 : : }
104 : : }
105 : : #endif
106 : :
107 [ + - ][ - + ]: 250 : rval = IntxUtils::EdgeIntersections2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, markb, markr, P, nP );MB_CHK_ERR( rval );
[ # # ][ # # ]
108 : : #ifdef ENABLE_DEBUG
109 : : if( dbg_1 )
110 : : {
111 : : for( int k = 0; k < 3; k++ )
112 : : {
113 : : std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
114 : : }
115 : : }
116 : : #endif
117 : :
118 : 250 : int side[MAXEDGES] = { 0 }; // this refers to what side? src or tgt?
119 : : int extraPoints =
120 [ + - ]: 250 : IntxUtils::borderPointsOfXinY2( srcCoords2D, nsSrc, tgtCoords2D, nsTgt, &( P[2 * nP] ), side, epsilon_area );
121 [ + + ]: 250 : if( extraPoints >= 1 )
122 : : {
123 [ + + ]: 360 : for( int k = 0; k < nsSrc; k++ )
124 : : {
125 [ + + ]: 270 : if( side[k] )
126 : : {
127 : : // this means that vertex k of src is inside convex tgt; mark edges k-1 and k in
128 : : // src,
129 : : // as being "intersected" by tgt; (even though they might not be intersected by
130 : : // other edges, the fact that their apex is inside, is good enough)
131 : 92 : markb[k] = 1;
132 : 92 : markb[( k + nsSrc - 1 ) % nsSrc] =
133 : 92 : 1; // it is the previous edge, actually, but instead of doing -1, it is
134 : : // better to do modulo +3 (modulo 4)
135 : : // null side b for next call
136 : 92 : side[k] = 0;
137 : : }
138 : : }
139 : : }
140 : : #ifdef ENABLE_DEBUG
141 : : if( dbg_1 )
142 : : {
143 : : for( int k = 0; k < 3; k++ )
144 : : {
145 : : std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
146 : : }
147 : : }
148 : : #endif
149 : 250 : nP += extraPoints;
150 : :
151 : : extraPoints =
152 [ + - ]: 250 : IntxUtils::borderPointsOfXinY2( tgtCoords2D, nsTgt, srcCoords2D, nsSrc, &( P[2 * nP] ), side, epsilon_area );
153 [ + + ]: 250 : if( extraPoints >= 1 )
154 : : {
155 [ + + ]: 600 : for( int k = 0; k < nsTgt; k++ )
156 : : {
157 [ + + ]: 450 : if( side[k] )
158 : : {
159 : : // this is to mark that tgt edges k-1 and k are intersecting src
160 : 175 : markr[k] = 1;
161 : 175 : markr[( k + nsTgt - 1 ) % nsTgt] =
162 : 175 : 1; // it is the previous edge, actually, but instead of doing -1, it is
163 : : // better to do modulo +3 (modulo 4)
164 : : // null side b for next call
165 : : }
166 : : }
167 : : }
168 : : #ifdef ENABLE_DEBUG
169 : : if( dbg_1 )
170 : : {
171 : : for( int k = 0; k < 3; k++ )
172 : : {
173 : : std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
174 : : }
175 : : }
176 : : #endif
177 : 250 : nP += extraPoints;
178 : :
179 : : // now sort and orient the points in P, such that they are forming a convex polygon
180 : : // this will be the foundation of our new mesh
181 : : // this works if the polygons are convex
182 [ + - ]: 250 : IntxUtils::SortAndRemoveDoubles2( P, nP, epsilon_1 ); // nP should be at most 8 in the end ?
183 : : // if there are more than 3 points, some area will be positive
184 : :
185 [ + + ]: 250 : if( nP >= 3 )
186 : : {
187 [ + + ]: 663 : for( int k = 1; k < nP - 1; k++ )
188 [ + - ]: 415 : area += IntxUtils::area2D( P, &P[2 * k], &P[2 * k + 2] );
189 : : }
190 : :
191 : 250 : return MB_SUCCESS; // no error
192 : : }
193 : :
194 : : // this method will also construct the triangles/polygons in the new mesh
195 : : // if we accept planar polygons, we just save them
196 : : // also, we could just create new vertices every time, and merge only in the end;
197 : : // could be too expensive, and the tolerance for merging could be an
198 : : // interesting topic
199 : 200 : ErrorCode Intx2MeshInPlane::findNodes( EntityHandle tgt, int nsTgt, EntityHandle src, int nsSrc, double* iP, int nP )
200 : : {
201 : : // except for gnomonic projection, everything is the same as spherical intx
202 : : // start copy
203 : : // first of all, check against tgt and src vertices
204 : : //
205 : : #ifdef ENABLE_DEBUG
206 : : if( dbg_1 )
207 : : {
208 : : std::cout << "tgt, src, nP, P " << mb->id_from_handle( tgt ) << " " << mb->id_from_handle( src ) << " " << nP
209 : : << "\n";
210 : : for( int n = 0; n < nP; n++ )
211 : : std::cout << " \t" << iP[2 * n] << "\t" << iP[2 * n + 1] << "\n";
212 : : }
213 : : #endif
214 : :
215 : : // get the edges for the tgt triangle; the extra points will be on those edges, saved as
216 : : // lists (unordetgt)
217 : :
218 : : // first get the list of edges adjacent to the tgt cell
219 : : // use the neighTgtEdgeTag
220 : : EntityHandle adjTgtEdges[MAXEDGES];
221 [ + - ][ - + ]: 200 : ErrorCode rval = mb->tag_get_data( neighTgtEdgeTag, &tgt, 1, &( adjTgtEdges[0] ) );MB_CHK_SET_ERR( rval, "can't get edge tgt tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
222 : : // we know that we have only nsTgt edges here; [nsTgt, MAXEDGES) are ignored, but it is small
223 : : // potatoes
224 : :
225 : : // these will be in the new mesh, mbOut
226 : : // some of them will be handles to the initial vertices from src or tgt meshes (lagr or euler)
227 : :
228 [ + - ][ + - ]: 200 : EntityHandle* foundIds = new EntityHandle[nP];
229 [ + + ]: 940 : for( int i = 0; i < nP; i++ )
230 : : {
231 : 740 : double* pp = &iP[2 * i]; // iP+2*i
232 : : //
233 [ + - ]: 740 : CartVect pos( pp[0], pp[1], 0. );
234 : 740 : int found = 0;
235 : : // first, are they on vertices from tgt or src?
236 : : // priority is the tgt mesh (mb2?)
237 : 740 : int j = 0;
238 : 740 : EntityHandle outNode = (EntityHandle)0;
239 [ + + ][ + + ]: 2834 : for( j = 0; j < nsTgt && !found; j++ )
240 : : {
241 : : // int node = tgtTri.v[j];
242 [ + - ]: 2094 : double d2 = IntxUtils::dist2( pp, &tgtCoords2D[2 * j] );
243 [ + + ]: 2094 : if( d2 < epsilon_1 )
244 : : {
245 : :
246 : 125 : foundIds[i] = tgtConn[j]; // no new node
247 : 125 : found = 1;
248 : : #ifdef ENABLE_DEBUG
249 : : if( dbg_1 )
250 : : std::cout << " tgt node j:" << j << " id:" << mb->id_from_handle( tgtConn[j] )
251 : : << " 2d coords:" << tgtCoords2D[2 * j] << " " << tgtCoords2D[2 * j + 1] << " d2: " << d2
252 : : << " \n";
253 : : #endif
254 : : }
255 : : }
256 : :
257 [ + + ][ + + ]: 2520 : for( j = 0; j < nsSrc && !found; j++ )
258 : : {
259 : : // int node = srcTri.v[j];
260 [ + - ]: 1780 : double d2 = IntxUtils::dist2( pp, &srcCoords2D[2 * j] );
261 [ + + ]: 1780 : if( d2 < epsilon_1 )
262 : : {
263 : : // suspect is srcConn[j] corresponding in mbOut
264 : :
265 : 69 : foundIds[i] = srcConn[j]; // no new node
266 : 69 : found = 1;
267 : : #ifdef ENABLE_DEBUG
268 : : if( dbg_1 )
269 : : std::cout << " src node " << j << " " << mb->id_from_handle( srcConn[j] ) << " d2:" << d2 << " \n";
270 : : #endif
271 : : }
272 : : }
273 [ + + ]: 740 : if( !found )
274 : : {
275 : : // find the edge it belongs, first, on the tgt element
276 : : //
277 [ + + ]: 2184 : for( j = 0; j < nsTgt; j++ )
278 : : {
279 : 1638 : int j1 = ( j + 1 ) % nsTgt;
280 [ + - ]: 1638 : double area = IntxUtils::area2D( &tgtCoords2D[2 * j], &tgtCoords2D[2 * j1], pp );
281 : : #ifdef ENABLE_DEBUG
282 : : if( dbg_1 )
283 : : std::cout << " edge " << j << ": " << mb->id_from_handle( adjTgtEdges[j] ) << " " << tgtConn[j]
284 : : << " " << tgtConn[j1] << " area : " << area << "\n";
285 : : #endif
286 [ + + ]: 1638 : if( fabs( area ) < epsilon_1 / 2 )
287 : : {
288 : : // found the edge; now find if there is a point in the list here
289 : : // std::vector<EntityHandle> * expts = extraNodesMap[tgtEdges[j]];
290 [ + - ]: 546 : int indx = TgtEdges.index( adjTgtEdges[j] );
291 [ - + ]: 546 : if( indx < 0 ) // CID 181166 (#1 of 1): Argument cannot be negative (NEGATIVE_RETURNS)
292 : : {
293 [ # # ][ # # ]: 0 : std::cerr << " error in adjacent tgt edge: " << mb->id_from_handle( adjTgtEdges[j] ) << "\n";
[ # # ][ # # ]
294 [ # # ]: 0 : delete[] foundIds;
295 : 0 : return MB_FAILURE;
296 : : }
297 [ + - ]: 546 : std::vector< EntityHandle >* expts = extraNodesVec[indx];
298 : : // if the points pp is between extra points, then just give that id
299 : : // if not, create a new point, (check the id)
300 : : // get the coordinates of the extra points so far
301 : 546 : int nbExtraNodesSoFar = expts->size();
302 [ + + ]: 546 : if( nbExtraNodesSoFar > 0 )
303 : : {
304 [ + - ][ + - ]: 1714 : CartVect* coords1 = new CartVect[nbExtraNodesSoFar];
[ + - ][ + + ]
305 [ + - ][ + - ]: 481 : mb->get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
[ + - ]
306 : : // std::list<int>::iterator it;
307 [ + + ][ + + ]: 1382 : for( int k = 0; k < nbExtraNodesSoFar && !found; k++ )
308 : : {
309 : : // int pnt = *it;
310 [ + - ][ + - ]: 901 : double d3 = ( pos - coords1[k] ).length_squared();
311 [ + + ]: 901 : if( d3 < epsilon_1 )
312 : : {
313 : 393 : found = 1;
314 [ + - ]: 393 : foundIds[i] = ( *expts )[k];
315 : : #ifdef ENABLE_DEBUG
316 : : if( dbg_1 ) std::cout << " found node:" << foundIds[i] << std::endl;
317 : : #endif
318 : : }
319 : : }
320 [ + - ]: 481 : delete[] coords1;
321 : : }
322 : :
323 [ + + ]: 546 : if( !found )
324 : : {
325 : : // create a new point in 2d (at the intersection)
326 : : // foundIds[i] = m_num2dPoints;
327 : : // expts.push_back(m_num2dPoints);
328 : : // need to create a new node in mbOut
329 : : // this will be on the edge, and it will be added to the local list
330 [ + - ][ + - ]: 153 : mb->create_vertex( pos.array(), outNode );
331 [ + - ]: 153 : ( *expts ).push_back( outNode );
332 : 153 : foundIds[i] = outNode;
333 : 546 : found = 1;
334 : : #ifdef ENABLE_DEBUG
335 : : if( dbg_1 ) std::cout << " new node: " << outNode << std::endl;
336 : : #endif
337 : : }
338 : : }
339 : : }
340 : : }
341 [ - + ]: 740 : if( !found )
342 : : {
343 [ # # ]: 0 : std::cout << " tgt polygon: ";
344 [ # # ]: 0 : for( int j1 = 0; j1 < nsTgt; j1++ )
345 : : {
346 [ # # ][ # # ]: 0 : std::cout << tgtCoords2D[2 * j1] << " " << tgtCoords2D[2 * j1 + 1] << "\n";
[ # # ][ # # ]
347 : : }
348 [ # # ][ # # ]: 0 : std::cout << " a point pp is not on a tgt polygon " << *pp << " " << pp[1] << " tgt polygon "
[ # # ][ # # ]
[ # # ]
349 [ # # ][ # # ]: 0 : << mb->id_from_handle( tgt ) << " \n";
[ # # ]
350 [ # # ]: 0 : delete[] foundIds;
351 : 0 : return MB_FAILURE;
352 : : }
353 : : }
354 : : #ifdef ENABLE_DEBUG
355 : : if( dbg_1 )
356 : : {
357 : : std::cout << " candidate polygon: nP " << nP << "\n";
358 : : for( int i1 = 0; i1 < nP; i1++ )
359 : : std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
360 : : }
361 : : #endif
362 : : // first, find out if we have nodes collapsed; shrink them
363 : : // we may have to reduce nP
364 : : // it is possible that some nodes are collapsed after intersection only
365 : : // nodes will always be in order (convex intersection)
366 [ + - ]: 200 : correct_polygon( foundIds, nP );
367 : : // now we can build the triangles, from P array, with foundIds
368 : : // we will put them in the out set
369 [ + + ]: 200 : if( nP >= 3 )
370 : : {
371 : : EntityHandle polyNew;
372 [ + - ]: 198 : mb->create_element( MBPOLYGON, foundIds, nP, polyNew );
373 [ + - ]: 198 : mb->add_entities( outSet, &polyNew, 1 );
374 : :
375 : : // tag it with the index ids from tgt and src sets
376 [ + - ]: 198 : int id = rs1.index( src ); // index starts from 0
377 [ + - ]: 198 : mb->tag_set_data( srcParentTag, &polyNew, 1, &id );
378 [ + - ]: 198 : id = rs2.index( tgt );
379 [ + - ]: 198 : mb->tag_set_data( tgtParentTag, &polyNew, 1, &id );
380 : :
381 : 198 : counting++;
382 [ + - ]: 198 : mb->tag_set_data( countTag, &polyNew, 1, &counting );
383 : :
384 : : #ifdef ENABLE_DEBUG
385 : : if( dbg_1 )
386 : : {
387 : :
388 : : std::cout << "Count: " << counting + 1 << "\n";
389 : : std::cout << " polygon " << mb->id_from_handle( polyNew ) << " nodes: " << nP << " :";
390 : : for( int i1 = 0; i1 < nP; i1++ )
391 : : std::cout << " " << mb->id_from_handle( foundIds[i1] );
392 : : std::cout << "\n";
393 : : std::vector< CartVect > posi( nP );
394 : : mb->get_coords( foundIds, nP, &( posi[0][0] ) );
395 : : for( int i1 = 0; i1 < nP; i1++ )
396 : : std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << posi[i1] << "\n";
397 : :
398 : : std::stringstream fff;
399 : : fff << "file0" << counting << ".vtk";
400 : : mb->write_mesh( fff.str().c_str(), &outSet, 1 );
401 : : }
402 : : #endif
403 : : }
404 [ + - ]: 200 : delete[] foundIds;
405 : 200 : foundIds = NULL;
406 : 200 : return MB_SUCCESS;
407 : : // end copy
408 : : }
409 : :
410 [ + - ][ + - ]: 4 : } // end namespace moab
|
