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