Branch data Line data Source code
1 : : /*
2 : : * Intx2Mesh.cpp
3 : : *
4 : : * Created on: Oct 2, 2012
5 : : */
6 : :
7 : : #include "moab/IntxMesh/Intx2Mesh.hpp"
8 : : #ifdef MOAB_HAVE_MPI
9 : : #include "moab/ParallelComm.hpp"
10 : : #include "MBParallelConventions.h"
11 : : #include "moab/ParallelMergeMesh.hpp"
12 : : #endif /* MOAB_HAVE_MPI */
13 : : #include "MBTagConventions.hpp"
14 : : // this is for DBL_MAX
15 : : #include <float.h>
16 : : #include <queue>
17 : : #include <sstream>
18 : : #include "moab/GeomUtil.hpp"
19 : : #include "moab/AdaptiveKDTree.hpp"
20 : :
21 : : namespace moab
22 : : {
23 : :
24 : : #ifdef ENABLE_DEBUG
25 : : int Intx2Mesh::dbg_1 = 0;
26 : : #endif
27 : :
28 : 3 : Intx2Mesh::Intx2Mesh( Interface* mbimpl )
29 : : : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
30 : : countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), tgtConn( NULL ),
31 : : srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ), my_rank( 0 )
32 : : #ifdef MOAB_HAVE_MPI
33 : : ,
34 : : parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
35 : : #endif
36 : : ,
37 [ + - ][ + - ]: 63 : max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
[ + - ][ + + ]
[ + - ][ + + ]
[ + - ][ + - ]
[ + - ][ + - ]
38 : : {
39 [ + - ]: 3 : gid = mbimpl->globalId_tag();
40 : 3 : }
41 : :
42 : 6 : Intx2Mesh::~Intx2Mesh()
43 : : {
44 : : // TODO Auto-generated destructor stub
45 : : #ifdef MOAB_HAVE_MPI
46 [ - + ]: 3 : if( remote_cells )
47 : : {
48 [ # # ]: 0 : delete remote_cells;
49 : 0 : remote_cells = NULL;
50 : : }
51 : : #endif
52 [ - + ]: 3 : }
53 : 6 : ErrorCode Intx2Mesh::FindMaxEdgesInSet( EntityHandle eset, int& max_edges )
54 : : {
55 [ + - ]: 6 : Range cells;
56 [ + - ][ - + ]: 6 : ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval );
[ # # ][ # # ]
57 : :
58 : 6 : max_edges = 0; // can be 0 for point clouds
59 [ + - ][ + - ]: 2906 : for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
[ + - ][ + - ]
[ + + ]
60 : : {
61 [ + - ]: 2900 : EntityHandle cell = *cit;
62 : : const EntityHandle* conn4;
63 : 2900 : int nnodes = 3;
64 [ + - ][ - + ]: 2900 : rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
65 [ + + ]: 2900 : if( nnodes > max_edges ) max_edges = nnodes;
66 : : }
67 : : // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
68 : : // global) but it is better to be consistent
69 : : #ifdef MOAB_HAVE_MPI
70 [ - + ]: 6 : if( parcomm )
71 : : {
72 : 0 : int local_max_edges = max_edges;
73 : : // now reduce max_edges over all processors
74 : : int mpi_err =
75 [ # # ][ # # ]: 0 : MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
[ # # ]
76 [ # # ]: 0 : if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
77 : : }
78 : : #endif
79 : :
80 : 6 : return MB_SUCCESS;
81 : : }
82 : 3 : ErrorCode Intx2Mesh::FindMaxEdges( EntityHandle set1, EntityHandle set2 )
83 : : {
84 [ - + ][ # # ]: 3 : ErrorCode rval = FindMaxEdgesInSet( set1, max_edges_1 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 1" );
[ # # ][ # # ]
[ # # ][ # # ]
85 [ - + ][ # # ]: 3 : rval = FindMaxEdgesInSet( set2, max_edges_2 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 2" );
[ # # ][ # # ]
[ # # ][ # # ]
86 : :
87 : 3 : return MB_SUCCESS;
88 : : }
89 : :
90 : 3 : ErrorCode Intx2Mesh::createTags()
91 : : {
92 [ - + ][ # # ]: 3 : if( tgtParentTag ) mb->tag_delete( tgtParentTag );
93 [ - + ][ # # ]: 3 : if( srcParentTag ) mb->tag_delete( srcParentTag );
94 [ - + ][ # # ]: 3 : if( countTag ) mb->tag_delete( countTag );
95 : :
96 : 3 : unsigned char def_data_bit = 0; // unused by default
97 : : // maybe the tgt tag is better to be deleted every time, and recreated;
98 : : // or is it easy to set all values to something again? like 0?
99 [ + - ][ - + ]: 3 : ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "can't get tgt flag tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
100 : :
101 : : // create tgt edges if they do not exist yet; so when they are looked upon, they are found
102 : : // this is the only call that is potentially NlogN, in the whole method
103 [ + - ][ - + ]: 3 : rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
104 : :
105 : : // now, create a map from each edge to a list of potential new nodes on a tgt edge
106 : : // this memory has to be cleaned up
107 : : // change it to a vector, and use the index in range of tgt edges
108 : 3 : int indx = 0;
109 [ + - ][ + - ]: 3 : extraNodesVec.resize( TgtEdges.size() );
110 [ + - ][ + - ]: 3217 : for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
[ + - ][ + - ]
[ + + ]
111 : : {
112 [ + - ][ + - ]: 3214 : std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
113 [ + - ]: 3214 : extraNodesVec[indx] = nv;
114 : : }
115 : :
116 : 3 : int defaultInt = -1;
117 : :
118 : : rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
119 [ + - ][ - + ]: 3 : &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
120 : :
121 : : rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
122 [ + - ][ - + ]: 3 : &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
123 : :
124 [ + - ][ - + ]: 3 : rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
125 : :
126 : : // for each cell in set 1, determine its neigh in set 1 (could be null too)
127 : : // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
128 [ + - ][ - + ]: 3 : rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
129 [ + - ][ - + ]: 3 : rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 2" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
130 : :
131 : : // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
132 : : // them each time edges were for sure created before (tgtEdges)
133 [ + - ]: 3 : std::vector< EntityHandle > zeroh( max_edges_2, 0 );
134 : : // if we have a tag with this name, it could be of a different size, so delete it
135 [ + - ]: 3 : rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
136 [ - + ][ # # ]: 3 : if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
[ # # ]
137 : : rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
138 [ + - ][ + - ]: 3 : MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
139 [ + - ][ + - ]: 1617 : for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
[ + - ][ + - ]
[ + + ]
140 : : {
141 [ + - ]: 1614 : EntityHandle tgtCell = *rit;
142 : 1614 : int num_nodes = 0;
143 [ + - ][ - + ]: 1614 : rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
144 : : // account for padded polygons
145 [ - + ][ # # ]: 1614 : while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
146 : 0 : num_nodes--;
147 : :
148 : 1614 : int i = 0;
149 [ + + ]: 8022 : for( i = 0; i < num_nodes; i++ )
150 : : {
151 : 6408 : EntityHandle v[2] = { tgtConn[i],
152 : 6408 : tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
153 [ + - ]: 6408 : std::vector< EntityHandle > adj_entities;
154 [ + - ]: 6408 : rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
155 [ + - ][ - + ]: 6408 : if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
[ - + ]
156 [ + - ][ + - ]: 6408 : zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
[ + - ]
157 : : // also, even if number of edges is less than max_edges_2, they will be ignored, even if
158 : : // the tag is dense
159 : 6408 : }
160 : : // now set the value of the tag
161 [ + - ][ + - ]: 1614 : rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
162 : : }
163 : 3 : return MB_SUCCESS;
164 : : }
165 : :
166 : 6 : ErrorCode Intx2Mesh::DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag )
167 : : {
168 [ + - ]: 6 : Range cells;
169 [ + - ][ - + ]: 6 : ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells in set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
170 : :
171 [ + - ]: 12 : std::vector< EntityHandle > neighbors( max_edges );
172 [ + - ]: 12 : std::vector< EntityHandle > zeroh( max_edges, 0 );
173 : : // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
174 : : // one on src mesh
175 [ + - ][ + - ]: 6 : rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create neighbors tag" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
176 : :
177 [ + - ][ + - ]: 2906 : for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
[ + - ][ + - ]
[ + + ]
178 : : {
179 [ + - ]: 2900 : EntityHandle cell = *cit;
180 : 2900 : int nnodes = 3;
181 : : // will get the nnodes ordered neighbors;
182 : : // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
183 : : const EntityHandle* conn4;
184 [ + - ][ - + ]: 2900 : rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
185 : 2900 : int nsides = nnodes;
186 : : // account for possible padded polygons
187 [ - + ][ # # ]: 2900 : while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
188 : 0 : nsides--;
189 : :
190 [ + + ]: 14340 : for( int i = 0; i < nsides; i++ )
191 : : {
192 : : EntityHandle v[2];
193 : 11440 : v[0] = conn4[i];
194 : 11440 : v[1] = conn4[( i + 1 ) % nsides];
195 : : // get all cells adjacent to these 2 vertices on the edge
196 [ + - ]: 11440 : std::vector< EntityHandle > adjcells;
197 [ + - ]: 22880 : std::vector< EntityHandle > cellsInSet;
[ + - + ]
198 [ + - ][ - + ]: 11440 : rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't adjacency to 2 verts" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
199 : : // now look for the cells contained in the input set;
200 : : // the input set should be a correct mesh, not overlapping cells, and manifold
201 : 11440 : size_t siz = adjcells.size();
202 [ + + ]: 34284 : for( size_t j = 0; j < siz; j++ )
203 [ + - ][ + - ]: 22844 : if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
[ + - ][ + - ]
[ + - ]
204 : 11440 : siz = cellsInSet.size();
205 : :
206 [ - + ]: 11440 : if( siz > 2 )
207 : : {
208 [ # # ][ # # ]: 0 : std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
[ # # ][ # # ]
209 [ # # ][ # # ]: 0 : << "\n";MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" ); // non-manifold
[ # # ][ # # ]
[ # # ][ # # ]
210 : : }
211 [ + + ]: 11440 : if( siz == 1 )
212 : : {
213 : : // it must be the border of the input mesh;
214 [ + - ]: 36 : neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border
215 : : // borders do not appear for a sphere in serial, but they do appear for
216 : : // parallel processing anyway
217 : 36 : continue;
218 : : }
219 : : // here siz ==2, it is either the first or second
220 [ + - ][ + + ]: 11404 : if( cell == cellsInSet[0] )
221 [ + - ][ + - ]: 5702 : neighbors[i] = cellsInSet[1];
222 : : else
223 [ + - ][ + - ]: 11440 : neighbors[i] = cellsInSet[0];
[ + - + ]
224 : 11440 : }
225 : : // fill the rest with 0
226 [ + + ]: 2980 : for( int i = nsides; i < max_edges; i++ )
227 [ + - ]: 80 : neighbors[i] = 0;
228 : : // now simply set the neighbors tag; the last few positions will not be used, but for
229 : : // simplicity will keep them all (MAXEDGES)
230 [ + - ][ + - ]: 2900 : rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" );
[ - + ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
231 : : }
232 : 12 : return MB_SUCCESS;
233 : : }
234 : :
235 : : // slow interface; this will not do the advancing front trick
236 : : // some are triangles, some are quads, some are polygons ...
237 : 0 : ErrorCode Intx2Mesh::intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
238 : : {
239 : : ErrorCode rval;
240 : 0 : mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
241 : 0 : mbs2 = mbset2;
242 : 0 : outSet = outputSet;
243 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
244 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
[ # # ][ # # ]
245 : : // from create tags, copy relevant ones
246 [ # # ][ # # ]: 0 : if( tgtParentTag ) mb->tag_delete( tgtParentTag );
247 [ # # ][ # # ]: 0 : if( srcParentTag ) mb->tag_delete( srcParentTag );
248 [ # # ][ # # ]: 0 : if( countTag ) mb->tag_delete( countTag );
249 : :
250 : : // create tgt edges if they do not exist yet; so when they are looked upon, they are found
251 : : // this is the only call that is potentially NlogN, in the whole method
252 [ # # ][ # # ]: 0 : rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
253 : :
254 : 0 : int indx = 0;
255 [ # # ][ # # ]: 0 : extraNodesVec.resize( TgtEdges.size() );
256 [ # # ][ # # ]: 0 : for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
[ # # ][ # # ]
[ # # ]
257 : : {
258 [ # # ][ # # ]: 0 : std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
259 [ # # ]: 0 : extraNodesVec[indx] = nv;
260 : : }
261 : :
262 : 0 : int defaultInt = -1;
263 : :
264 : : rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
265 [ # # ][ # # ]: 0 : &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
266 : :
267 : : rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
268 [ # # ][ # # ]: 0 : &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
269 : :
270 [ # # ][ # # ]: 0 : rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
271 : :
272 : : // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
273 : : // them each time edges were for sure created before (tgtEdges)
274 : : // if we have a tag with this name, it could be of a different size, so delete it
275 [ # # ]: 0 : rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
276 [ # # ][ # # ]: 0 : if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
[ # # ]
277 [ # # ]: 0 : std::vector< EntityHandle > zeroh( max_edges_2, 0 );
278 : : rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
279 [ # # ][ # # ]: 0 : MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
280 [ # # ][ # # ]: 0 : for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
[ # # ][ # # ]
[ # # ]
281 : : {
282 [ # # ]: 0 : EntityHandle tgtCell = *rit;
283 : 0 : int num_nodes = 0;
284 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
285 : : // account for padded polygons
286 [ # # ][ # # ]: 0 : while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
287 : 0 : num_nodes--;
288 : :
289 : 0 : int i = 0;
290 [ # # ]: 0 : for( i = 0; i < num_nodes; i++ )
291 : : {
292 : 0 : EntityHandle v[2] = { tgtConn[i],
293 : 0 : tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons
294 [ # # ]: 0 : std::vector< EntityHandle > adj_entities;
295 [ # # ]: 0 : rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
296 [ # # ][ # # ]: 0 : if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error
[ # # ]
297 [ # # ][ # # ]: 0 : zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes
[ # # ]
298 : : // also, even if number of edges is less than max_edges_2, they will be ignored, even if
299 : : // the tag is dense
300 : 0 : }
301 : : // now set the value of the tag
302 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
303 : : }
304 : :
305 : : // create the kd tree on source cells, and intersect all targets in an expensive loop
306 : : // build a kd tree with the rs1 (source) cells
307 [ # # ]: 0 : AdaptiveKDTree kd( mb );
308 : 0 : EntityHandle tree_root = 0;
309 [ # # ][ # # ]: 0 : rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
[ # # ][ # # ]
310 : :
311 : : // find out max edge on source mesh;
312 : 0 : double max_length = 0;
313 : : {
314 [ # # ]: 0 : std::vector< double > coords;
315 [ # # ]: 0 : coords.resize( 3 * max_edges_1 );
316 [ # # ][ # # ]: 0 : for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ )
[ # # ][ # # ]
[ # # ][ # # ]
317 : : {
318 : 0 : const EntityHandle* conn = NULL;
319 : : int nnodes;
320 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
321 [ # # ][ # # ]: 0 : while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
322 : 0 : nnodes--;
323 [ # # ][ # # ]: 0 : rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coordinates" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
324 [ # # ]: 0 : for( int j = 0; j < nnodes; j++ )
325 : : {
326 : 0 : int next = ( j + 1 ) % nnodes;
327 : : double leng;
328 [ # # ][ # # ]: 0 : leng = ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
[ # # ][ # # ]
329 [ # # ][ # # ]: 0 : ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
[ # # ][ # # ]
330 [ # # ][ # # ]: 0 : ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
[ # # ][ # # ]
331 : 0 : leng = sqrt( leng );
332 [ # # ]: 0 : if( leng > max_length ) max_length = leng;
333 : : }
334 : 0 : }
335 : : }
336 : : // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
337 : 0 : double tolerance = 1.e-15;
338 [ # # ]: 0 : if( max_length < 1. )
339 : : {
340 : : // basically, the sag for an arc of length max_length on a circle of radius 1
341 : 0 : tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
342 : 0 : tolerance = 3 * tolerance; // why ?
343 [ # # ]: 0 : if( !my_rank )
344 [ # # ][ # # ]: 0 : std::cout << " max edge length: " << max_length << " tolerance for kd tree: " << tolerance << "\n";
[ # # ][ # # ]
[ # # ]
345 : : }
346 [ # # ][ # # ]: 0 : for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
[ # # ][ # # ]
[ # # ]
347 : : {
348 [ # # ]: 0 : EntityHandle tcell = *it;
349 : : // find vertex positions
350 : 0 : const EntityHandle* conn = NULL;
351 : 0 : int nnodes = 0;
352 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
353 : : // find leaves close to those positions
354 [ # # ]: 0 : double areaTgtCell = setup_tgt_cell( tcell, nnodes ); // this is the area in the gnomonic plane
355 : 0 : double recoveredArea = 0;
356 [ # # ]: 0 : std::vector< double > positions;
357 [ # # ]: 0 : positions.resize( nnodes * 3 );
358 [ # # ][ # # ]: 0 : rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
359 : :
360 : : // distance to search will be based on average edge length
361 : 0 : double av_len = 0;
362 [ # # ]: 0 : for( int k = 0; k < nnodes; k++ )
363 : : {
364 : 0 : int ik = ( k + 1 ) % nnodes;
365 : 0 : double len1 = 0;
366 [ # # ]: 0 : for( int j = 0; j < 3; j++ )
367 : : {
368 [ # # ][ # # ]: 0 : double len2 = positions[3 * k + j] - positions[3 * ik + j];
369 : 0 : len1 += len2 * len2;
370 : : }
371 : 0 : av_len += sqrt( len1 );
372 : : }
373 [ # # ]: 0 : if( nnodes > 0 ) av_len /= nnodes;
374 : : // find leaves within a distance from each vertex of target
375 : : // in those leaves, collect all cells; we will try for an intx in there
376 [ # # ][ # # ]: 0 : Range close_source_cells;
377 [ # # ][ # # ]: 0 : std::vector< EntityHandle > leaves;
378 [ # # ]: 0 : for( int i = 0; i < nnodes; i++ )
379 : : {
380 : :
381 : 0 : leaves.clear();
382 [ # # ][ # # ]: 0 : rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
383 : :
384 [ # # ][ # # ]: 0 : for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
[ # # ]
385 : : {
386 [ # # ]: 0 : Range tmp;
387 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
388 : :
389 [ # # ][ # # ]: 0 : close_source_cells.merge( tmp.begin(), tmp.end() );
[ # # ][ # # ]
390 : 0 : }
391 : : }
392 : : #ifdef VERBOSE
393 : : if( close_source_cells.empty() )
394 : : {
395 : : std::cout << " there are no close source cells to target cell " << tcell << " id from handle "
396 : : << mb->id_from_handle( tcell ) << "\n";
397 : : }
398 : : #endif
399 [ # # ][ # # ]: 0 : for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
[ # # ][ # # ]
[ # # ]
400 : : {
401 [ # # ]: 0 : EntityHandle startSrc = *it2;
402 : 0 : double area = 0;
403 : : // if area is > 0 , we have intersections
404 : : double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
405 : : //
406 : 0 : int nP = 0;
407 : : int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
408 : : int nsTgt, nsSrc;
409 [ # # ][ # # ]: 0 : rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
[ # # ][ # # ]
410 [ # # ]: 0 : if( area > 0 )
411 : : {
412 [ # # ]: 0 : if( nP > 1 )
413 : : { // this will also construct triangles/polygons in the new mesh, if needed
414 [ # # ][ # # ]: 0 : rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
[ # # ][ # # ]
415 : : }
416 : 0 : recoveredArea += area;
417 : : }
418 : : }
419 [ # # ]: 0 : recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fract
420 : 0 : }
421 : : // before cleaning up , we need to settle the position of the intersection points
422 : : // on the boundary edges
423 : : // this needs to be collective, so we should maybe wait something
424 : : #ifdef MOAB_HAVE_MPI
425 [ # # ][ # # ]: 0 : rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
426 : : #endif
427 : :
428 [ # # ]: 0 : this->clean();
429 : 0 : return MB_SUCCESS;
430 : : }
431 : : // main interface; this will do the advancing front trick
432 : : // some are triangles, some are quads, some are polygons ...
433 : 3 : ErrorCode Intx2Mesh::intersect_meshes( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
434 : : {
435 : : ErrorCode rval;
436 : 3 : mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
437 : 3 : mbs2 = mbset2;
438 : 3 : outSet = outputSet;
439 : : #ifdef VERBOSE
440 : : std::stringstream ffs, fft;
441 : : ffs << "source_rank0" << my_rank << ".vtk";
442 : : rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
443 : : fft << "target_rank0" << my_rank << ".vtk";
444 : : rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );
445 : :
446 : : #endif
447 : : // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
448 : :
449 : 3 : EntityHandle startSrc = 0, startTgt = 0;
450 : :
451 [ + - ][ - + ]: 3 : rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
452 [ + - ][ - + ]: 3 : rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
[ # # ][ # # ]
453 : : // std::cout << "rs1.size() = " << rs1.size() << " and rs2.size() = " << rs2.size() << "\n";
454 : : // std::cout.flush();
455 : :
456 [ + - ]: 3 : createTags(); // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
457 : :
458 [ + - ]: 3 : Range rs22 = rs2; // a copy of the initial range; we will remove from it elements as we
459 : : // advance ; rs2 is needed for marking the polygon to the tgt parent
460 : :
461 : : // create the local kdd tree with source elements; will use it to search
462 : : // more efficiently for the seeds in advancing front;
463 : : // some of the target cells will not be covered by source cells, and they need to be eliminated
464 : : // early from contention
465 : :
466 : : // build a kd tree with the rs1 (source) cells
467 [ + - ]: 6 : AdaptiveKDTree kd( mb );
468 : 3 : EntityHandle tree_root = 0;
469 [ + - ][ - + ]: 3 : rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
[ # # ][ # # ]
470 : :
471 [ + - ][ + + ]: 6 : while( !rs22.empty() )
472 : : {
473 : : #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
474 : : if( rs22.size() < rs2.size() )
475 : : {
476 : : std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
477 : : << "\n";
478 : : std::stringstream ffo;
479 : : ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
480 : : rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
481 : : }
482 : : #endif
483 : 3 : bool seedFound = false;
484 [ + - ][ # # ]: 3 : for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
[ + - ][ + - ]
[ + - ]
485 : : {
486 [ + - ]: 3 : startTgt = *it;
487 : 3 : int found = 0;
488 : : // find vertex positions
489 : 3 : const EntityHandle* conn = NULL;
490 : 3 : int nnodes = 0;
491 [ + - ][ - + ]: 6 : rval = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
492 : : // find leaves close to those positions
493 [ + - ]: 3 : std::vector< double > positions;
494 [ + - ]: 3 : positions.resize( nnodes * 3 );
495 [ + - ][ + - ]: 3 : rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
496 : : // find leaves within a distance from each vertex of target
497 : : // in those leaves, collect all cells; we will try for an intx in there, instead of
498 : : // looping over all rs1 cells, as before
499 [ + - ]: 6 : Range close_source_cells;
[ - - + ]
500 [ + - ]: 6 : std::vector< EntityHandle > leaves;
[ - - + ]
501 [ + + ]: 14 : for( int i = 0; i < nnodes; i++ )
502 : : {
503 : :
504 : 11 : leaves.clear();
505 [ + - ][ + - ]: 11 : rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
506 : :
507 [ + - ][ + - ]: 22 : for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
[ + + ]
508 : : {
509 [ + - ]: 11 : Range tmp;
510 [ + - ][ + - ]: 11 : rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
511 : :
512 [ + - ][ + - ]: 11 : close_source_cells.merge( tmp.begin(), tmp.end() );
[ + - ][ + - ]
513 : 11 : }
514 : : }
515 : :
516 [ + - ][ + - ]: 7 : for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
[ + - # # ]
517 : : {
518 [ + - ]: 7 : startSrc = *it2;
519 : 7 : double area = 0;
520 : : // if area is > 0 , we have intersections
521 : : double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon
522 : : //
523 : 7 : int nP = 0;
524 : : int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
525 : : int nsTgt, nsSrc;
526 : : rval =
527 [ + - ][ - + ]: 7 : computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
[ # # ][ # # ]
528 [ + + ]: 7 : if( area > 0 )
529 : : {
530 : 3 : found = 1;
531 : 3 : seedFound = true;
532 : 3 : break; // found 2 elements that intersect; these will be the seeds
533 : : }
534 : : }
535 [ + - ]: 3 : if( found )
536 : 3 : break;
537 : : else
538 : : {
539 : : #if defined( VERBOSE )
540 : : std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
541 : : << " not intx with any source\n";
542 : : #endif
543 [ # # ]: 3 : rs22.erase( startTgt );
[ - - + ]
544 : : }
545 : 0 : }
546 [ - + ]: 3 : if( !seedFound ) continue; // continue while(!rs22.empty())
547 : :
548 [ + - ][ + - ]: 3 : std::queue< EntityHandle > srcQueue; // these are corresponding to Ta,
549 [ + - ]: 3 : srcQueue.push( startSrc );
550 [ + - ][ + - ]: 6 : std::queue< EntityHandle > tgtQueue;
[ + - # # ]
551 [ + - ]: 3 : tgtQueue.push( startTgt );
552 : :
553 [ + - ][ + - ]: 6 : Range toResetSrcs; // will be used to reset src flags for every tgt quad
554 : : // processed
555 : :
556 : : /*if (my_rank==0)
557 : : dbg_1 = 1;*/
558 : 3 : unsigned char used = 1;
559 : : // mark the start tgt quad as used, so it will not come back again
560 [ + - ][ - + ]: 3 : rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
[ # # ][ # # ]
561 [ + - ][ + + ]: 1617 : while( !tgtQueue.empty() )
[ + - ]
562 : : {
563 : : // flags for the side : 0 means a src cell not found on side
564 : : // a paired src not found yet for the neighbors of tgt
565 [ + - ][ + + ]: 35508 : Range nextSrc[MAXEDGES]; // there are new ranges of possible next src cells for
[ + - # # ]
566 : : // seeding the side j of tgt cell
567 : :
568 [ + - ]: 1614 : EntityHandle currentTgt = tgtQueue.front();
569 [ + - ]: 1614 : tgtQueue.pop();
570 : : int nsidesTgt; // will be initialized now
571 [ + - ]: 1614 : double areaTgtCell = setup_tgt_cell( currentTgt, nsidesTgt ); // this is the area in the gnomonic plane
572 : 1614 : double recoveredArea = 0;
573 : : // get the neighbors of tgt, and if they are solved already, do not bother with that
574 : : // side of tgt
575 : 1614 : EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
576 [ + - ][ - + ]: 1614 : rval = mb->tag_get_data( tgtNeighTag, ¤tTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
577 : : #ifdef ENABLE_DEBUG
578 : : if( dbg_1 )
579 : : {
580 : : std::cout << "Next: neighbors for current tgt ";
581 : : for( int kk = 0; kk < nsidesTgt; kk++ )
582 : : {
583 : : if( tgtNeighbors[kk] > 0 )
584 : : std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
585 : : else
586 : : std::cout << 0 << " ";
587 : : }
588 : : std::cout << std::endl;
589 : : }
590 : : #endif
591 : : // now get the status of neighbors; if already solved, make them 0, so not to bother
592 : : // anymore on that side of tgt
593 [ + + ]: 8022 : for( int j = 0; j < nsidesTgt; j++ )
594 : : {
595 : 6408 : EntityHandle tgtNeigh = tgtNeighbors[j];
596 : 6408 : unsigned char status = 1;
597 [ + + ]: 6408 : if( tgtNeigh == 0 ) continue;
598 [ + - ][ - + ]: 6388 : rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval ); // status 0 is unused
[ # # ][ # # ]
599 [ + + ]: 6388 : if( 1 == status ) tgtNeighbors[j] = 0; // so will not look anymore on this side of tgt
600 : : }
601 : :
602 : : #ifdef ENABLE_DEBUG
603 : : if( dbg_1 )
604 : : {
605 : : std::cout << "reset sources: ";
606 : : for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
607 : : std::cout << mb->id_from_handle( *itr ) << " ";
608 : : std::cout << std::endl;
609 : : }
610 : : #endif
611 [ + - ]: 1614 : EntityHandle currentSrc = srcQueue.front();
612 : : // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
613 : : // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
614 : : // pairs;
615 : : // at every moment, the queue contains pairs of cells that intersect, and they form the
616 : : // "advancing front"
617 [ + - ]: 1614 : srcQueue.pop();
618 [ + - ]: 1614 : toResetSrcs.clear(); // empty the range of used srcs, will have to be set unused again,
619 : : // at the end of tgt element processing
620 [ + - ]: 1614 : toResetSrcs.insert( currentSrc );
621 : : // mb2->set_tag_data
622 [ + - ][ + - ]: 19368 : std::queue< EntityHandle > localSrc;
[ + + ]
623 [ + - ]: 1614 : localSrc.push( currentSrc );
624 : : #ifdef VERBOSE
625 : : int countingStart = counting;
626 : : #endif
627 : : // will advance-front search in the neighborhood of tgt cell, until we finish processing
628 : : // all
629 : : // possible src cells; localSrc queue will contain all possible src cells that cover
630 : : // the current tgt cell
631 [ + - ][ + + ]: 8365 : while( !localSrc.empty() )
632 : : {
633 : : //
634 [ + - ]: 6751 : EntityHandle srcT = localSrc.front();
635 [ + - ]: 6751 : localSrc.pop();
636 : : double P[10 * MAXEDGES], area; //
637 : 6751 : int nP = 0;
638 : 6751 : int nb[MAXEDGES] = { 0 };
639 : 6751 : int nr[MAXEDGES] = { 0 };
640 : :
641 : : int nsidesSrc; ///
642 : : // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
643 : : // intersection points could include the vertices of initial elements
644 : : // nb [j] = 0 means no intersection on the side j for element src (markers)
645 : : // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
646 : : // tgt poly. A potential next poly in the tgt queue is the tgt poly that is
647 : : // adjacent to this side
648 : : rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
649 [ + - ][ - + ]: 6751 : nsidesTgt );MB_CHK_ERR( rval );
[ # # ][ # # ]
650 [ + - ]: 6751 : if( nP > 0 )
651 : : {
652 : : #ifdef ENABLE_DEBUG
653 : : if( dbg_1 )
654 : : {
655 : : for( int k = 0; k < 3; k++ )
656 : : {
657 : : std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
658 : : }
659 : : }
660 : : #endif
661 : :
662 : : // intersection found: output P and original triangles if nP > 2
663 : 6751 : EntityHandle neighbors[MAXEDGES] = { 0 };
664 [ + - ]: 6751 : rval = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
665 [ - + ]: 6751 : if( rval != MB_SUCCESS )
666 : : {
667 [ # # ][ # # ]: 0 : std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
[ # # ]
668 : 6751 : return MB_FAILURE;
669 : : }
670 : :
671 : : // add neighbors to the localSrc queue, if they are not marked
672 [ + + ]: 33352 : for( int nn = 0; nn < nsidesSrc; nn++ )
673 : : {
674 : 26601 : EntityHandle neighbor = neighbors[nn];
675 [ + + ][ + + ]: 26601 : if( neighbor > 0 && nb[nn] > 0 ) // advance across src boundary nn
676 : : {
677 [ + - ][ + - ]: 13255 : if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
[ + - ][ + + ]
678 : : {
679 [ + - ]: 5137 : localSrc.push( neighbor );
680 : : #ifdef ENABLE_DEBUG
681 : : if( dbg_1 )
682 : : {
683 : : std::cout << " local src elem " << mb->id_from_handle( neighbor )
684 : : << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
685 : : mb->list_entities( &neighbor, 1 );
686 : : }
687 : : #endif
688 [ + - ]: 5137 : toResetSrcs.insert( neighbor );
689 : : }
690 : : }
691 : : }
692 : : // n(find(nc>0))=ac; % ac is starting candidate for neighbor
693 [ + + ]: 33554 : for( int nn = 0; nn < nsidesTgt; nn++ )
694 : : {
695 [ + + ][ + + ]: 26803 : if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
696 [ + - ]: 3978 : nextSrc[nn].insert( srcT ); // potential src cell that can intersect
697 : : // the tgt neighbor nn
698 : : }
699 [ + + ]: 6751 : if( nP > 1 )
700 : : { // this will also construct triangles/polygons in the new mesh, if needed
701 [ + - ][ - + ]: 6671 : rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
[ # # ][ # # ]
702 : : }
703 : :
704 : 6751 : recoveredArea += area;
705 : : }
706 : : #ifdef ENABLE_DEBUG
707 : : else if( dbg_1 )
708 : : {
709 : : std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
710 : : << mb->id_from_handle( srcT ) << "\n";
711 : : }
712 : : #endif
713 : : } // end while (!localSrc.empty())
714 : 1614 : recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fraction
715 : : #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
716 : : if( fabs( recoveredArea ) > epsilon_1 )
717 : : {
718 : : #ifdef VERBOSE
719 : : std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
720 : : << " fraction error recovery:" << recoveredArea
721 : : << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart
722 : : << "\n";
723 : : #endif
724 : : }
725 : : #endif
726 : : // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
727 : : // mesh)
728 [ + - ]: 1614 : rs22.erase( currentTgt );
729 : : // also, look at its neighbors, and add to the seeds a next one
730 : :
731 [ + + ][ + - ]: 8022 : for( int j = 0; j < nsidesTgt; j++ )
732 : : {
733 : 6408 : EntityHandle tgtNeigh = tgtNeighbors[j];
734 [ + + ][ + - ]: 6408 : if( tgtNeigh == 0 || nextSrc[j].size() == 0 ) // if tgt is bigger than src, there could be no src
[ - + ][ + + ]
735 : : // to advance on that side
736 : 4778 : continue;
737 : 1630 : int nsidesTgt2 = 0;
738 : : setup_tgt_cell( tgtNeigh,
739 [ + - ]: 1630 : nsidesTgt2 ); // find possible intersection with src cell from nextSrc
740 [ + - ][ + - ]: 3437 : for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
[ + - ][ + - ]
[ + + ]
741 : : {
742 [ + - ]: 1807 : EntityHandle nextB = *nit;
743 : : // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
744 : : // src quad
745 : : double P[10 * MAXEDGES], area; //
746 : 1807 : int nP = 0;
747 : 1807 : int nb[MAXEDGES] = { 0 };
748 : 1807 : int nr[MAXEDGES] = { 0 };
749 : :
750 : : int nsidesSrc; ///
751 : : rval = computeIntersectionBetweenTgtAndSrc(
752 [ + - ][ - + ]: 1807 : /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
[ # # ][ # # ]
753 [ + + ]: 1807 : if( area > 0 )
754 : : {
755 [ + - ]: 1611 : tgtQueue.push( tgtNeigh );
756 [ + - ]: 1611 : srcQueue.push( nextB );
757 : : #ifdef ENABLE_DEBUG
758 : : if( dbg_1 )
759 : : std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
760 : : << mb->id_from_handle( nextB ) << std::endl;
761 : : #endif
762 [ + - ][ - + ]: 1611 : rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
[ # # ][ # # ]
763 : 1611 : break; // so we are done with this side of tgt, we have found a proper next
764 : : // seed
765 : : }
766 : : }
767 : : }
768 : :
769 [ # # ]: 1614 : } // end while (!tgtQueue.empty())
770 : 3 : }
771 : : #ifdef ENABLE_DEBUG
772 : : if( dbg_1 )
773 : : {
774 : : for( int k = 0; k < 6; k++ )
775 : : mout_1[k].close();
776 : : }
777 : : #endif
778 : : // before cleaning up , we need to settle the position of the intersection points
779 : : // on the boundary edges
780 : : // this needs to be collective, so we should maybe wait something
781 : : #ifdef MOAB_HAVE_MPI
782 [ + - ][ - + ]: 3 : rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
783 : : #endif
784 : :
785 [ + - ]: 3 : this->clean();
786 : 6 : return MB_SUCCESS;
787 : : }
788 : :
789 : : // clean some memory allocated
790 : 3 : void Intx2Mesh::clean()
791 : : {
792 : : //
793 : 3 : int indx = 0;
794 [ + - ][ + - ]: 3217 : for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
[ + - ][ + - ]
[ + + ]
795 : : {
796 [ + - ][ + - ]: 3214 : delete extraNodesVec[indx];
797 : : }
798 : : // extraNodesMap.clear();
799 : 3 : extraNodesVec.clear();
800 : : // also, delete some bit tags, used to mark processed tgts and srcs
801 : 3 : mb->tag_delete( TgtFlagTag );
802 : 3 : counting = 0; // reset counting to original value
803 : 3 : }
804 : : // this method will reduce number of nodes, collapse edges that are of length 0
805 : : // so a polygon like 428 431 431 will become a line 428 431
806 : : // or something like 428 431 431 531 -> 428 431 531
807 : 6671 : void Intx2Mesh::correct_polygon( EntityHandle* nodes, int& nP )
808 : : {
809 : 6671 : int i = 0;
810 [ + + ]: 30877 : while( i < nP )
811 : : {
812 : 24206 : int nextIndex = ( i + 1 ) % nP;
813 [ + + ]: 24206 : if( nodes[i] == nodes[nextIndex] )
814 : : {
815 : : #ifdef ENABLE_DEBUG
816 : : // we need to reduce nP, and collapse nodes
817 : : if( dbg_1 )
818 : : {
819 : : std::cout << " nodes duplicated in list: ";
820 : : for( int j = 0; j < nP; j++ )
821 : : std::cout << nodes[j] << " ";
822 : : std::cout << "\n";
823 : : std::cout << " node " << nodes[i] << " at index " << i << " is duplicated"
824 : : << "\n";
825 : : }
826 : : #endif
827 : : // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
828 : : // then next thing does nothing
829 : : // (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
830 [ + + ]: 12 : for( int k = i; k < nP - 1; k++ )
831 : 8 : nodes[k] = nodes[k + 1];
832 : 4 : nP--; // decrease the number of nodes; also, decrease i, just if we may need to check
833 : : // again
834 : 4 : i--;
835 : : }
836 : 24206 : i++;
837 : : }
838 : 6671 : return;
839 : : }
840 : : #ifdef MOAB_HAVE_MPI
841 : 0 : ErrorCode Intx2Mesh::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts )
842 : : {
843 [ # # ]: 0 : localEnts.clear();
844 [ # # ][ # # ]: 0 : ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
[ # # ][ # # ]
845 : :
846 [ # # ]: 0 : rval = mb->get_connectivity( localEnts, local_verts );
847 [ # # ][ # # ]: 0 : int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
[ # # ][ # # ]
848 : :
849 [ # # ]: 0 : assert( parcomm != NULL );
850 : :
851 : : // get the position of local vertices, and decide local boxes (allBoxes...)
852 : 0 : double bmin[3] = { DBL_MAX, DBL_MAX, DBL_MAX };
853 : 0 : double bmax[3] = { -DBL_MAX, -DBL_MAX, -DBL_MAX };
854 : :
855 [ # # ]: 0 : std::vector< double > coords( 3 * num_local_verts );
856 [ # # ][ # # ]: 0 : rval = mb->get_coords( local_verts, &coords[0] );ERRORR( rval, "can't get coords of vertices " );
[ # # ][ # # ]
[ # # ]
857 : :
858 [ # # ]: 0 : for( int i = 0; i < num_local_verts; i++ )
859 : : {
860 [ # # ]: 0 : for( int k = 0; k < 3; k++ )
861 : : {
862 [ # # ]: 0 : double val = coords[3 * i + k];
863 [ # # ]: 0 : if( val < bmin[k] ) bmin[k] = val;
864 [ # # ]: 0 : if( val > bmax[k] ) bmax[k] = val;
865 : : }
866 : : }
867 [ # # ][ # # ]: 0 : int numprocs = parcomm->proc_config().proc_size();
868 [ # # ]: 0 : allBoxes.resize( 6 * numprocs );
869 : :
870 [ # # ][ # # ]: 0 : my_rank = parcomm->proc_config().proc_rank();
871 [ # # ]: 0 : for( int k = 0; k < 3; k++ )
872 : : {
873 [ # # ]: 0 : allBoxes[6 * my_rank + k] = bmin[k];
874 [ # # ]: 0 : allBoxes[6 * my_rank + 3 + k] = bmax[k];
875 : : }
876 : :
877 : : // now communicate to get all boxes
878 : : int mpi_err;
879 : : #if( MPI_VERSION >= 2 )
880 : : // use "in place" option
881 [ # # ]: 0 : mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
882 [ # # ][ # # ]: 0 : parcomm->proc_config().proc_comm() );
[ # # ]
883 : : #else
884 : : {
885 : : std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() );
886 : : mpi_err = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
887 : : parcomm->proc_config().proc_comm() );
888 : : allBoxes = allBoxes_tmp;
889 : : }
890 : : #endif
891 [ # # ]: 0 : if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
892 : :
893 [ # # ]: 0 : if( my_rank == 0 )
894 : : {
895 [ # # ][ # # ]: 0 : std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
[ # # ][ # # ]
896 [ # # ]: 0 : << " on second mesh \n";
897 [ # # ]: 0 : for( int i = 0; i < numprocs; i++ )
898 : : {
899 [ # # ][ # # ]: 0 : std::cout << "proc: " << i << " box min: " << allBoxes[6 * i] << " " << allBoxes[6 * i + 1] << " "
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
900 [ # # ][ # # ]: 0 : << allBoxes[6 * i + 2] << " \n";
[ # # ]
901 [ # # ][ # # ]: 0 : std::cout << " box max: " << allBoxes[6 * i + 3] << " " << allBoxes[6 * i + 4] << " "
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
902 [ # # ][ # # ]: 0 : << allBoxes[6 * i + 5] << " \n";
[ # # ]
903 : : }
904 : : }
905 : :
906 : 0 : return MB_SUCCESS;
907 : : }
908 : 0 : ErrorCode Intx2Mesh::create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set )
909 : : {
910 : : // compute the bounding box on each proc
911 [ # # ]: 0 : assert( parcomm != NULL );
912 : :
913 [ # # ]: 0 : localEnts.clear();
914 [ # # ][ # # ]: 0 : ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
[ # # ][ # # ]
915 : :
916 : 0 : Tag dpTag = 0;
917 [ # # ]: 0 : std::string tag_name( "DP" );
918 [ # # ][ # # ]: 0 : rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );ERRORR( rval, "can't get DP tag" );
[ # # ][ # # ]
919 : :
920 : 0 : EntityHandle dum = 0;
921 : : Tag corrTag;
922 [ # # ][ # # ]: 0 : rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );ERRORR( rval, "can't get CORR tag" );
[ # # ][ # # ]
923 : : // get all local verts
924 [ # # ]: 0 : Range local_verts;
925 [ # # ]: 0 : rval = mb->get_connectivity( localEnts, local_verts );
926 [ # # ][ # # ]: 0 : int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
[ # # ][ # # ]
927 : :
928 [ # # ][ # # ]: 0 : rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );ERRORR( rval, "can't build processor boxes" );
[ # # ][ # # ]
929 : :
930 [ # # ]: 0 : std::vector< int > gids( num_local_verts );
931 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
[ # # ][ # # ]
[ # # ]
932 : :
933 : : // now see the departure points; to what boxes should we send them?
934 [ # # ]: 0 : std::vector< double > dep_points( 3 * num_local_verts );
935 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( dpTag, local_verts, (void*)&dep_points[0] );ERRORR( rval, "can't get DP tag values" );
[ # # ][ # # ]
[ # # ]
936 : : // ranges to send to each processor; will hold vertices and elements (quads?)
937 : : // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
938 [ # # ]: 0 : std::map< int, Range > Rto;
939 [ # # ][ # # ]: 0 : int numprocs = parcomm->proc_config().proc_size();
940 : :
941 [ # # ][ # # ]: 0 : for( Range::iterator eit = localEnts.begin(); eit != localEnts.end(); ++eit )
[ # # ][ # # ]
[ # # ]
942 : : {
943 [ # # ]: 0 : EntityHandle q = *eit;
944 : : const EntityHandle* conn4;
945 : : int num_nodes;
946 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
[ # # ][ # # ]
947 [ # # ]: 0 : CartVect qbmin( DBL_MAX );
948 [ # # ]: 0 : CartVect qbmax( -DBL_MAX );
949 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
950 : : {
951 : 0 : EntityHandle v = conn4[i];
952 [ # # ][ # # ]: 0 : size_t index = local_verts.find( v ) - local_verts.begin();
[ # # ]
953 [ # # ][ # # ]: 0 : CartVect dp( &dep_points[3 * index] ); // will use constructor
954 [ # # ]: 0 : for( int j = 0; j < 3; j++ )
955 : : {
956 [ # # ][ # # ]: 0 : if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
[ # # ][ # # ]
[ # # ]
957 [ # # ][ # # ]: 0 : if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
[ # # ][ # # ]
[ # # ]
958 : : }
959 : : }
960 [ # # ]: 0 : for( int p = 0; p < numprocs; p++ )
961 : : {
962 [ # # ][ # # ]: 0 : CartVect bbmin( &allBoxes[6 * p] );
963 [ # # ][ # # ]: 0 : CartVect bbmax( &allBoxes[6 * p + 3] );
964 [ # # ][ # # ]: 0 : if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) ) { Rto[p].insert( q ); }
[ # # ][ # # ]
965 : : }
966 : : }
967 : :
968 : : // now, build TLv and TLq, for each p
969 : 0 : size_t numq = 0;
970 : 0 : size_t numv = 0;
971 [ # # ]: 0 : for( int p = 0; p < numprocs; p++ )
972 : : {
973 [ # # ]: 0 : if( p == (int)my_rank ) continue; // do not "send" it, because it is already here
974 [ # # ]: 0 : Range& range_to_P = Rto[p];
975 : : // add the vertices to it
976 [ # # ][ # # ]: 0 : if( range_to_P.empty() ) continue; // nothing to send to proc p
977 [ # # ]: 0 : Range vertsToP;
978 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
[ # # ][ # # ]
979 [ # # ]: 0 : numq = numq + range_to_P.size();
980 [ # # ]: 0 : numv = numv + vertsToP.size();
981 [ # # ][ # # ]: 0 : range_to_P.merge( vertsToP );
982 : 0 : }
983 [ # # ]: 0 : TupleList TLv;
984 [ # # ]: 0 : TupleList TLq;
985 [ # # ]: 0 : TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, DP points
986 [ # # ]: 0 : TLv.enableWriteAccess();
987 : :
988 : 0 : int sizeTuple = 2 + max_edges_1; // determined earlier, for src, first mesh
989 : 0 : TLq.initialize( 2 + max_edges_1, 0, 1, 0,
990 [ # # ]: 0 : numq ); // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh
991 [ # # ]: 0 : TLq.enableWriteAccess();
992 : : #ifdef VERBOSE
993 : : std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
994 : : #endif
995 [ # # ]: 0 : for( int to_proc = 0; to_proc < numprocs; to_proc++ )
996 : : {
997 [ # # ]: 0 : if( to_proc == (int)my_rank ) continue;
998 [ # # ]: 0 : Range& range_to_P = Rto[to_proc];
999 [ # # ]: 0 : Range V = range_to_P.subset_by_type( MBVERTEX );
1000 : :
1001 [ # # ][ # # ]: 0 : for( Range::iterator it = V.begin(); it != V.end(); ++it )
[ # # ][ # # ]
[ # # ]
1002 : : {
1003 [ # # ]: 0 : EntityHandle v = *it;
1004 [ # # ][ # # ]: 0 : unsigned int index = local_verts.find( v ) - local_verts.begin();
[ # # ]
1005 [ # # ]: 0 : int n = TLv.get_n();
1006 : 0 : TLv.vi_wr[2 * n] = to_proc; // send to processor
1007 [ # # ]: 0 : TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range
1008 [ # # ]: 0 : TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i]
1009 [ # # ]: 0 : TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
1010 [ # # ]: 0 : TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
1011 [ # # ]: 0 : TLv.inc_n();
1012 : : }
1013 : : // also, prep the quad for sending ...
1014 [ # # ][ # # ]: 0 : Range Q = range_to_P.subset_by_dimension( 2 );
1015 [ # # ][ # # ]: 0 : for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
1016 : : {
1017 [ # # ]: 0 : EntityHandle q = *it;
1018 : : int global_id;
1019 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
[ # # ][ # # ]
1020 [ # # ]: 0 : int n = TLq.get_n();
1021 : 0 : TLq.vi_wr[sizeTuple * n] = to_proc; //
1022 : 0 : TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ...
1023 : : const EntityHandle* conn4;
1024 : : int num_nodes;
1025 : : rval = mb->get_connectivity( q, conn4,
1026 [ # # ]: 0 : num_nodes ); // could be up to MAXEDGES, but it is limited by max_edges_1
1027 [ # # ][ # # ]: 0 : ERRORR( rval, "can't get connectivity for cell" );
[ # # ]
1028 [ # # ][ # # ]: 0 : if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
[ # # ]
1029 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
1030 : : {
1031 : 0 : EntityHandle v = conn4[i];
1032 [ # # ][ # # ]: 0 : unsigned int index = local_verts.find( v ) - local_verts.begin();
[ # # ]
1033 [ # # ]: 0 : TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1034 : : }
1035 [ # # ]: 0 : for( int k = num_nodes; k < max_edges_1; k++ )
1036 : : {
1037 : 0 : TLq.vi_wr[sizeTuple * n + 2 + k] =
1038 : 0 : 0; // fill the rest of node ids with 0; we know that the node ids start from 1!
1039 : : }
1040 : 0 : TLq.vul_wr[n] = q; // save here the entity handle, it will be communicated back
1041 : : // maybe we should forget about global ID
1042 [ # # ]: 0 : TLq.inc_n();
1043 : : }
1044 : 0 : }
1045 : :
1046 : : // now we are done populating the tuples; route them to the appropriate processors
1047 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
[ # # ]
1048 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
[ # # ]
1049 : : // the elements are already in localEnts;
1050 : :
1051 : : // maps from global ids to new vertex and quad handles, that are added
1052 [ # # ]: 0 : std::map< int, EntityHandle > globalID_to_handle;
1053 : : /*std::map<int, EntityHandle> globalID_to_eh;*/
1054 : 0 : globalID_to_eh.clear(); // need for next iteration
1055 : : // now, look at every TLv, and see if we have to create a vertex there or not
1056 [ # # ]: 0 : int n = TLv.get_n(); // the size of the points received
1057 [ # # ]: 0 : for( int i = 0; i < n; i++ )
1058 : : {
1059 : 0 : int globalId = TLv.vi_rd[2 * i + 1];
1060 [ # # ][ # # ]: 0 : if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
[ # # ]
1061 : : {
1062 : : EntityHandle new_vert;
1063 : 0 : double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1064 [ # # ][ # # ]: 0 : rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
[ # # ][ # # ]
1065 [ # # ]: 0 : globalID_to_handle[globalId] = new_vert;
1066 : : }
1067 : : }
1068 : :
1069 : : // now, all dep points should be at their place
1070 : : // look in the local list of q for this proc, and create all those quads and vertices if needed
1071 : : // it may be an overkill, but because it does not involve communication, we do it anyway
1072 [ # # ]: 0 : Range& local = Rto[my_rank];
1073 [ # # ]: 0 : Range local_q = local.subset_by_dimension( 2 );
1074 : : // the local should have all the vertices in local_verts
1075 [ # # ][ # # ]: 0 : for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
[ # # ][ # # ]
[ # # ]
1076 : : {
1077 [ # # ]: 0 : EntityHandle q = *it;
1078 : : int nnodes;
1079 : : const EntityHandle* conn4;
1080 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( q, conn4, nnodes );ERRORR( rval, "can't get connectivity of local q " );
[ # # ][ # # ]
1081 : : EntityHandle new_conn[MAXEDGES];
1082 [ # # ]: 0 : for( int i = 0; i < nnodes; i++ )
1083 : : {
1084 : 0 : EntityHandle v1 = conn4[i];
1085 [ # # ][ # # ]: 0 : unsigned int index = local_verts.find( v1 ) - local_verts.begin();
[ # # ]
1086 [ # # ]: 0 : int globalId = gids[index];
1087 [ # # ][ # # ]: 0 : if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
[ # # ]
1088 : : {
1089 : : // we need to create that vertex, at this position dep_points
1090 [ # # ][ # # ]: 0 : double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] };
[ # # ]
1091 : : EntityHandle new_vert;
1092 [ # # ][ # # ]: 0 : rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
[ # # ][ # # ]
1093 [ # # ]: 0 : globalID_to_handle[globalId] = new_vert;
1094 : : }
1095 [ # # ][ # # ]: 0 : new_conn[i] = globalID_to_handle[gids[index]];
1096 : : }
1097 : : EntityHandle new_element;
1098 : : //
1099 : 0 : EntityType entType = MBQUAD;
1100 [ # # ]: 0 : if( nnodes > 4 ) entType = MBPOLYGON;
1101 [ # # ]: 0 : if( nnodes < 4 ) entType = MBTRI;
1102 : :
1103 [ # # ][ # # ]: 0 : rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new quad " );
[ # # ][ # # ]
1104 [ # # ][ # # ]: 0 : rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
[ # # ][ # # ]
1105 : : int gid_el;
1106 : : // get the global ID of the initial quad
1107 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
[ # # ][ # # ]
1108 [ # # ]: 0 : globalID_to_eh[gid_el] = new_element;
1109 : : // is this redundant or not?
1110 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( corrTag, &new_element, 1, &q );ERRORR( rval, "can't set corr tag on new el" );
[ # # ][ # # ]
1111 : : // set the global id on new elem
1112 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &new_element, 1, &gid_el );ERRORR( rval, "can't set global id tag on new el" );
[ # # ][ # # ]
1113 : : }
1114 : : // now look at all elements received through; we do not want to duplicate them
1115 [ # # ]: 0 : n = TLq.get_n(); // number of elements received by this processor
1116 : : // form the remote cells, that will be used to send the tracer info back to the originating proc
1117 [ # # ][ # # ]: 0 : remote_cells = new TupleList();
1118 [ # # ]: 0 : remote_cells->initialize( 2, 0, 1, 0, n ); // will not have tracer data anymore
1119 [ # # ]: 0 : remote_cells->enableWriteAccess();
1120 [ # # ]: 0 : for( int i = 0; i < n; i++ )
1121 : : {
1122 : 0 : int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1123 : 0 : int from_proc = TLq.vi_wr[sizeTuple * i];
1124 : : // do we already have a quad with this global ID, represented?
1125 [ # # ][ # # ]: 0 : if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
[ # # ]
1126 : : {
1127 : : // construct the conn quad
1128 : : EntityHandle new_conn[MAXEDGES];
1129 : 0 : int nnodes = -1;
1130 [ # # ]: 0 : for( int j = 0; j < max_edges_1; j++ )
1131 : : {
1132 : 0 : int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID
1133 [ # # ]: 0 : if( vgid == 0 )
1134 : 0 : new_conn[j] = 0;
1135 : : else
1136 : : {
1137 [ # # ][ # # ]: 0 : assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
[ # # ]
1138 [ # # ]: 0 : new_conn[j] = globalID_to_handle[vgid];
1139 : 0 : nnodes = j + 1; // nodes are at the beginning, and are variable number
1140 : : }
1141 : : }
1142 : : EntityHandle new_element;
1143 : : //
1144 : 0 : EntityType entType = MBQUAD;
1145 [ # # ]: 0 : if( nnodes > 4 ) entType = MBPOLYGON;
1146 [ # # ]: 0 : if( nnodes < 4 ) entType = MBTRI;
1147 [ # # ][ # # ]: 0 : rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
[ # # ][ # # ]
1148 [ # # ]: 0 : globalID_to_eh[globalIdEl] = new_element;
1149 [ # # ][ # # ]: 0 : rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
[ # # ][ # # ]
1150 : : /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q);ERRORR(rval, "can't set corr tag on new el");*/
1151 : 0 : remote_cells->vi_wr[2 * i] = from_proc;
1152 : 0 : remote_cells->vi_wr[2 * i + 1] = globalIdEl;
1153 : : // remote_cells->vr_wr[i] = 0.; // no contribution yet sent back
1154 : 0 : remote_cells->vul_wr[i] = TLq.vul_rd[i]; // this is the corresponding tgt cell (arrival)
1155 [ # # ]: 0 : remote_cells->inc_n();
1156 : : // set the global id on new elem
1157 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set global id tag on new el" );
[ # # ][ # # ]
1158 : : }
1159 : : }
1160 : : // order the remote cells tuple list, with the global id, because we will search in it
1161 : : // remote_cells->print("remote_cells before sorting");
1162 [ # # ]: 0 : moab::TupleList::buffer sort_buffer;
1163 [ # # ]: 0 : sort_buffer.buffer_init( n );
1164 [ # # ]: 0 : remote_cells->sort( 1, &sort_buffer );
1165 [ # # ]: 0 : sort_buffer.reset();
1166 : 0 : return MB_SUCCESS;
1167 : : }
1168 : :
1169 : : // this algorithm assumes lagr set is already created, and some elements will be coming from
1170 : : // other procs, and populate the covering_set
1171 : : // we need to keep in a tuple list the remote cells from other procs, because we need to send back
1172 : : // the intersection info (like area of the intx polygon, and the current concentration) maybe total
1173 : : // mass in that intx
1174 : 0 : ErrorCode Intx2Mesh::create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set )
1175 : : {
1176 : 0 : EntityHandle dum = 0;
1177 : :
1178 : : Tag corrTag;
1179 [ # # ]: 0 : ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
1180 : : // start copy from 2nd alg
1181 : : // compute the bounding box on each proc
1182 [ # # ]: 0 : assert( parcomm != NULL );
1183 [ # # ][ # # ]: 0 : if( 1 == parcomm->proc_config().proc_size() )
[ # # ]
1184 : : {
1185 : 0 : covering_set = lagr_set; // nothing to communicate, it must be serial
1186 : 0 : return MB_SUCCESS;
1187 : : }
1188 : :
1189 : : // get all local verts
1190 [ # # ]: 0 : Range local_verts;
1191 [ # # ]: 0 : rval = mb->get_connectivity( localEnts, local_verts );
1192 [ # # ][ # # ]: 0 : int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
[ # # ][ # # ]
1193 : :
1194 [ # # ]: 0 : std::vector< int > gids( num_local_verts );
1195 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
[ # # ][ # # ]
[ # # ]
1196 : :
1197 [ # # ]: 0 : Range localDepCells;
1198 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( lagr_set, 2, localDepCells );ERRORR( rval, "can't get ents by dimension from lagr set" );
[ # # ][ # # ]
1199 : :
1200 : : // get all lagr verts (departure vertices)
1201 [ # # ]: 0 : Range lagr_verts;
1202 [ # # ]: 0 : rval = mb->get_connectivity( localDepCells, lagr_verts ); // they should be created in
1203 : : // the same order as the euler vertices
1204 [ # # ][ # # ]: 0 : int num_lagr_verts = (int)lagr_verts.size();ERRORR( rval, "can't get local lagr vertices" );
[ # # ][ # # ]
1205 : :
1206 : : // now see the departure points position; to what boxes should we send them?
1207 [ # # ]: 0 : std::vector< double > dep_points( 3 * num_lagr_verts );
1208 [ # # ][ # # ]: 0 : rval = mb->get_coords( lagr_verts, &dep_points[0] );ERRORR( rval, "can't get departure points position" );
[ # # ][ # # ]
[ # # ]
1209 : : // ranges to send to each processor; will hold vertices and elements (quads?)
1210 : : // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
1211 [ # # ]: 0 : std::map< int, Range > Rto;
1212 [ # # ][ # # ]: 0 : int numprocs = parcomm->proc_config().proc_size();
1213 : :
1214 [ # # ][ # # ]: 0 : for( Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit )
[ # # ][ # # ]
[ # # ]
1215 : : {
1216 [ # # ]: 0 : EntityHandle q = *eit;
1217 : : const EntityHandle* conn4;
1218 : : int num_nodes;
1219 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
[ # # ][ # # ]
1220 [ # # ]: 0 : CartVect qbmin( DBL_MAX );
1221 [ # # ]: 0 : CartVect qbmax( -DBL_MAX );
1222 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
1223 : : {
1224 : 0 : EntityHandle v = conn4[i];
1225 [ # # ]: 0 : int index = lagr_verts.index( v );
1226 [ # # ]: 0 : assert( -1 != index );
1227 [ # # ][ # # ]: 0 : CartVect dp( &dep_points[3 * index] ); // will use constructor
1228 [ # # ]: 0 : for( int j = 0; j < 3; j++ )
1229 : : {
1230 [ # # ][ # # ]: 0 : if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
[ # # ][ # # ]
[ # # ]
1231 [ # # ][ # # ]: 0 : if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
[ # # ][ # # ]
[ # # ]
1232 : : }
1233 : : }
1234 [ # # ]: 0 : for( int p = 0; p < numprocs; p++ )
1235 : : {
1236 [ # # ][ # # ]: 0 : CartVect bbmin( &allBoxes[6 * p] );
1237 [ # # ][ # # ]: 0 : CartVect bbmax( &allBoxes[6 * p + 3] );
1238 [ # # ][ # # ]: 0 : if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) ) { Rto[p].insert( q ); }
[ # # ][ # # ]
1239 : : }
1240 : : }
1241 : :
1242 : : // now, build TLv and TLq, for each p
1243 : 0 : size_t numq = 0;
1244 : 0 : size_t numv = 0;
1245 [ # # ]: 0 : for( int p = 0; p < numprocs; p++ )
1246 : : {
1247 [ # # ]: 0 : if( p == (int)my_rank ) continue; // do not "send" it, because it is already here
1248 [ # # ]: 0 : Range& range_to_P = Rto[p];
1249 : : // add the vertices to it
1250 [ # # ][ # # ]: 0 : if( range_to_P.empty() ) continue; // nothing to send to proc p
1251 [ # # ]: 0 : Range vertsToP;
1252 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
[ # # ][ # # ]
1253 [ # # ]: 0 : numq = numq + range_to_P.size();
1254 [ # # ]: 0 : numv = numv + vertsToP.size();
1255 [ # # ][ # # ]: 0 : range_to_P.merge( vertsToP );
1256 : 0 : }
1257 [ # # ]: 0 : TupleList TLv;
1258 [ # # ]: 0 : TupleList TLq;
1259 [ # # ]: 0 : TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, DP points
1260 [ # # ]: 0 : TLv.enableWriteAccess();
1261 : :
1262 : 0 : int sizeTuple = 2 + max_edges_1; // max edges could be up to MAXEDGES :) for polygons
1263 : 0 : TLq.initialize( 2 + max_edges_1, 0, 1, 0,
1264 [ # # ]: 0 : numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v)
1265 : : // send also the corresponding tgt cell it will come to
1266 [ # # ]: 0 : TLq.enableWriteAccess();
1267 : : #ifdef VERBOSE
1268 : : std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
1269 : : #endif
1270 : :
1271 [ # # ]: 0 : for( int to_proc = 0; to_proc < numprocs; to_proc++ )
1272 : : {
1273 [ # # ]: 0 : if( to_proc == (int)my_rank ) continue;
1274 [ # # ]: 0 : Range& range_to_P = Rto[to_proc];
1275 [ # # ]: 0 : Range V = range_to_P.subset_by_type( MBVERTEX );
1276 : :
1277 [ # # ][ # # ]: 0 : for( Range::iterator it = V.begin(); it != V.end(); ++it )
[ # # ][ # # ]
[ # # ]
1278 : : {
1279 [ # # ]: 0 : EntityHandle v = *it;
1280 [ # # ]: 0 : int index = lagr_verts.index( v ); // will be the same index as the corresponding vertex in euler verts
1281 [ # # ]: 0 : assert( -1 != index );
1282 [ # # ]: 0 : int n = TLv.get_n();
1283 : 0 : TLv.vi_wr[2 * n] = to_proc; // send to processor
1284 [ # # ]: 0 : TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range
1285 [ # # ]: 0 : TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i]
1286 [ # # ]: 0 : TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
1287 [ # # ]: 0 : TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
1288 [ # # ]: 0 : TLv.inc_n();
1289 : : }
1290 : : // also, prep the 2d cells for sending ...
1291 [ # # ][ # # ]: 0 : Range Q = range_to_P.subset_by_dimension( 2 );
1292 [ # # ][ # # ]: 0 : for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
1293 : : {
1294 [ # # ]: 0 : EntityHandle q = *it; // this is a src cell
1295 : : int global_id;
1296 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
[ # # ][ # # ]
1297 [ # # ]: 0 : int n = TLq.get_n();
1298 : 0 : TLq.vi_wr[sizeTuple * n] = to_proc; //
1299 : 0 : TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ...
1300 : : const EntityHandle* conn4;
1301 : : int num_nodes;
1302 : : rval = mb->get_connectivity(
1303 [ # # ]: 0 : q, conn4, num_nodes ); // could be up to 10;ERRORR( rval, "can't get connectivity for quad" );
1304 [ # # ][ # # ]: 0 : if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
[ # # ]
1305 [ # # ]: 0 : for( int i = 0; i < num_nodes; i++ )
1306 : : {
1307 : 0 : EntityHandle v = conn4[i];
1308 [ # # ]: 0 : int index = lagr_verts.index( v );
1309 [ # # ]: 0 : assert( -1 != index );
1310 [ # # ]: 0 : TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1311 : : }
1312 [ # # ]: 0 : for( int k = num_nodes; k < max_edges_1; k++ )
1313 : : {
1314 : 0 : TLq.vi_wr[sizeTuple * n + 2 + k] =
1315 : 0 : 0; // fill the rest of node ids with 0; we know that the node ids start from 1!
1316 : : }
1317 : : EntityHandle tgtCell;
1318 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( corrTag, &q, 1, &tgtCell );ERRORR( rval, "can't get corresponding tgt cell for dep cell" );
[ # # ][ # # ]
1319 : 0 : TLq.vul_wr[n] = tgtCell; // this will be sent to remote_cells, to be able to come back
1320 [ # # ]: 0 : TLq.inc_n();
1321 : : }
1322 : 0 : }
1323 : : // now we can route them to each processor
1324 : : // now we are done populating the tuples; route them to the appropriate processors
1325 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
[ # # ]
1326 [ # # ][ # # ]: 0 : ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
[ # # ]
1327 : : // the elements are already in localEnts;
1328 : :
1329 : : // maps from global ids to new vertex and quad handles, that are added
1330 [ # # ]: 0 : std::map< int, EntityHandle > globalID_to_handle;
1331 : : // we already have vertices from lagr set; they are already in the processor, even before
1332 : : // receiving other verts from neighbors
1333 : 0 : int k = 0;
1334 [ # # ][ # # ]: 0 : for( Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ )
[ # # ][ # # ]
[ # # ]
1335 : : {
1336 [ # # ][ # # ]: 0 : globalID_to_handle[gids[k]] = *vit; // a little bit of overkill
[ # # ]
1337 : : // we do know that the global ids between euler and lagr verts are parallel
1338 : : }
1339 : : /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one?
1340 : 0 : globalID_to_eh.clear();
1341 : : // now, look at every TLv, and see if we have to create a vertex there or not
1342 [ # # ]: 0 : int n = TLv.get_n(); // the size of the points received
1343 [ # # ]: 0 : for( int i = 0; i < n; i++ )
1344 : : {
1345 : 0 : int globalId = TLv.vi_rd[2 * i + 1];
1346 [ # # ][ # # ]: 0 : if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
[ # # ]
1347 : : {
1348 : : EntityHandle new_vert;
1349 : 0 : double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1350 [ # # ][ # # ]: 0 : rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
[ # # ][ # # ]
1351 [ # # ]: 0 : globalID_to_handle[globalId] = new_vert;
1352 : : }
1353 : : }
1354 : :
1355 : : // now, all dep points should be at their place
1356 : : // look in the local list of 2d cells for this proc, and create all those cells if needed
1357 : : // it may be an overkill, but because it does not involve communication, we do it anyway
1358 [ # # ]: 0 : Range& local = Rto[my_rank];
1359 [ # # ]: 0 : Range local_q = local.subset_by_dimension( 2 );
1360 : : // the local should have all the vertices in lagr_verts
1361 [ # # ][ # # ]: 0 : for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
[ # # ][ # # ]
[ # # ]
1362 : : {
1363 [ # # ]: 0 : EntityHandle q = *it; // these are from lagr cells, local
1364 : : int gid_el;
1365 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
[ # # ][ # # ]
1366 [ # # ]: 0 : globalID_to_eh[gid_el] = q; // do we need this? maybe to just mark the ones on this processor
1367 : : // maybe a range of global cell ids is fine?
1368 : : }
1369 : : // now look at all elements received through; we do not want to duplicate them
1370 [ # # ]: 0 : n = TLq.get_n(); // number of elements received by this processor
1371 : : // a cell should be received from one proc only; so why are we so worried about duplicated
1372 : : // elements? a vertex can be received from multiple sources, that is fine
1373 : :
1374 [ # # ][ # # ]: 0 : remote_cells = new TupleList();
1375 [ # # ]: 0 : remote_cells->initialize( 2, 0, 1, 0, n ); // no tracers anymore in these tuples
1376 [ # # ]: 0 : remote_cells->enableWriteAccess();
1377 [ # # ]: 0 : for( int i = 0; i < n; i++ )
1378 : : {
1379 : 0 : int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1380 : 0 : int from_proc = TLq.vi_rd[sizeTuple * i];
1381 : : // do we already have a quad with this global ID, represented?
1382 [ # # ][ # # ]: 0 : if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
[ # # ]
1383 : : {
1384 : : // construct the conn quad
1385 : : EntityHandle new_conn[MAXEDGES];
1386 : 0 : int nnodes = -1;
1387 [ # # ]: 0 : for( int j = 0; j < max_edges_1; j++ )
1388 : : {
1389 : 0 : int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID
1390 [ # # ]: 0 : if( vgid == 0 )
1391 : 0 : new_conn[j] = 0;
1392 : : else
1393 : : {
1394 [ # # ][ # # ]: 0 : assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
[ # # ]
1395 [ # # ]: 0 : new_conn[j] = globalID_to_handle[vgid];
1396 : 0 : nnodes = j + 1; // nodes are at the beginning, and are variable number
1397 : : }
1398 : : }
1399 : : EntityHandle new_element;
1400 : : //
1401 : 0 : EntityType entType = MBQUAD;
1402 [ # # ]: 0 : if( nnodes > 4 ) entType = MBPOLYGON;
1403 [ # # ]: 0 : if( nnodes < 4 ) entType = MBTRI;
1404 [ # # ][ # # ]: 0 : rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
[ # # ][ # # ]
1405 [ # # ]: 0 : globalID_to_eh[globalIdEl] = new_element;
1406 [ # # ]: 0 : local_q.insert( new_element );
1407 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set gid on new element " );
[ # # ][ # # ]
1408 : : }
1409 : 0 : remote_cells->vi_wr[2 * i] = from_proc;
1410 : 0 : remote_cells->vi_wr[2 * i + 1] = globalIdEl;
1411 : : // remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
1412 : 0 : remote_cells->vul_wr[i] = TLq.vul_rd[i]; // this is the corresponding tgt cell (arrival)
1413 [ # # ]: 0 : remote_cells->inc_n();
1414 : : }
1415 : : // now, create a new set, covering_set
1416 [ # # ][ # # ]: 0 : rval = mb->create_meshset( MESHSET_SET, covering_set );ERRORR( rval, "can't create new mesh set " );
[ # # ][ # # ]
1417 [ # # ][ # # ]: 0 : rval = mb->add_entities( covering_set, local_q );ERRORR( rval, "can't add entities to new mesh set " );
[ # # ][ # # ]
1418 : : // order the remote cells tuple list, with the global id, because we will search in it
1419 : : // remote_cells->print("remote_cells before sorting");
1420 [ # # ]: 0 : moab::TupleList::buffer sort_buffer;
1421 [ # # ]: 0 : sort_buffer.buffer_init( n );
1422 [ # # ]: 0 : remote_cells->sort( 1, &sort_buffer );
1423 [ # # ]: 0 : sort_buffer.reset();
1424 : 0 : return MB_SUCCESS;
1425 : : // end copy
1426 : : }
1427 : :
1428 : 3 : ErrorCode Intx2Mesh::resolve_intersection_sharing()
1429 : : {
1430 [ - + ][ # # ]: 3 : if( parcomm && parcomm->size() > 1 )
[ - + ]
1431 : : {
1432 : : /*
1433 : : moab::ParallelMergeMesh pm(parcomm, epsilon_1);
1434 : : ErrorCode rval = pm.merge(outSet, false, 2); // resolve only the output set, do not skip
1435 : : local merge, use dim 2 ERRORR(rval, "can't merge intersection ");
1436 : : */
1437 : : // look at non-owned shared vertices, that could be part of original source set
1438 : : // they should be removed from intx set reference, because they might not have a
1439 : : // correspondent on the other task
1440 [ # # ]: 0 : Range nonOwnedVerts;
1441 [ # # ][ # # ]: 0 : Range vertsInIntx;
1442 [ # # ][ # # ]: 0 : Range intxCells;
1443 [ # # ][ # # ]: 0 : ErrorCode rval = mb->get_entities_by_dimension( outSet, 2, intxCells );MB_CHK_ERR( rval );
[ # # ][ # # ]
1444 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( intxCells, vertsInIntx );MB_CHK_ERR( rval );
[ # # ][ # # ]
1445 : :
1446 [ # # ][ # # ]: 0 : rval = parcomm->filter_pstatus( vertsInIntx, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonOwnedVerts );MB_CHK_ERR( rval );
[ # # ][ # # ]
1447 : :
1448 : : // some of these vertices can be in original set 1, which was covered, transported;
1449 : : // but they should not be "shared" from the intx point of view, because they are not shared
1450 : : // with another task they might have come from coverage as a plain vertex, so losing the
1451 : : // sharing property ?
1452 : :
1453 [ # # ][ # # ]: 0 : Range coverVerts;
1454 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( rs1, coverVerts );MB_CHK_ERR( rval );
[ # # ][ # # ]
1455 : : // find out those that are on the interface
1456 [ # # ][ # # ]: 0 : Range vertsCovInterface;
1457 [ # # ][ # # ]: 0 : rval = parcomm->filter_pstatus( coverVerts, PSTATUS_INTERFACE, PSTATUS_AND, -1, &vertsCovInterface );MB_CHK_ERR( rval );
[ # # ][ # # ]
1458 : : // how many of these are in
1459 [ # # ][ # # ]: 0 : Range nodesToDuplicate = intersect( vertsCovInterface, nonOwnedVerts );
1460 : : // first, get all cells connected to these vertices, from intxCells
1461 : :
1462 [ # # ][ # # ]: 0 : Range connectedCells;
1463 [ # # ][ # # ]: 0 : rval = mb->get_adjacencies( nodesToDuplicate, 2, false, connectedCells, Interface::UNION );MB_CHK_ERR( rval );
[ # # ][ # # ]
1464 : : // only those in intx set:
1465 [ # # ][ # # ]: 0 : connectedCells = intersect( connectedCells, intxCells );
1466 : : // first duplicate vertices in question:
1467 [ # # ][ # # ]: 0 : std::map< EntityHandle, EntityHandle > duplicatedVerticesMap;
1468 [ # # ][ # # ]: 0 : for( Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ )
[ # # ][ # # ]
[ # # ]
1469 : : {
1470 [ # # ]: 0 : EntityHandle vertex = *vit;
1471 : : double coords[3];
1472 [ # # ][ # # ]: 0 : rval = mb->get_coords( &vertex, 1, coords );MB_CHK_ERR( rval );
[ # # ][ # # ]
1473 : : EntityHandle newVertex;
1474 [ # # ][ # # ]: 0 : rval = mb->create_vertex( coords, newVertex );MB_CHK_ERR( rval );
[ # # ][ # # ]
1475 [ # # ]: 0 : duplicatedVerticesMap[vertex] = newVertex;
1476 : : }
1477 : :
1478 : : // look now at connectedCells, and change their connectivities:
1479 [ # # ][ # # ]: 0 : for( Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ )
[ # # ][ # # ]
[ # # ][ # # ]
1480 : : {
1481 [ # # ]: 0 : EntityHandle intxCell = *eit;
1482 : : // replace connectivity
1483 [ # # ]: 0 : std::vector< EntityHandle > connectivity;
1484 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( &intxCell, 1, connectivity );MB_CHK_ERR( rval );
[ # # ][ # # ]
1485 [ # # ]: 0 : for( size_t i = 0; i < connectivity.size(); i++ )
1486 : : {
1487 [ # # ]: 0 : EntityHandle currentVertex = connectivity[i];
1488 [ # # ]: 0 : std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex );
1489 [ # # ][ # # ]: 0 : if( mit != duplicatedVerticesMap.end() )
1490 : : {
1491 [ # # ][ # # ]: 0 : connectivity[i] = mit->second; // replace connectivity directly
1492 : : }
1493 : : }
1494 : 0 : int nnodes = (int)connectivity.size();
1495 [ # # ][ # # ]: 0 : rval = mb->set_connectivity( intxCell, &connectivity[0], nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ][ # # ]
1496 : 0 : }
1497 : : }
1498 : 3 : return MB_SUCCESS;
1499 : : }
1500 : : #endif /* MOAB_HAVE_MPI */
1501 [ + - ][ + - ]: 12 : } /* namespace moab */
|
