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