 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#include <iostream>#include <cstdlib>#include <vector>#include <string>#include <sstream>#include <cassert>#include "moab/Core.hpp"#include "moab/IntxMesh/IntxUtils.hpp"#include "moab/ProgOptions.hpp"#include "moab/CpuTimer.hpp"#include "DebugOutput.hpp"
Include dependency graph for mbIntxCheck.cpp:Go to the source code of this file.
Functions | |
| int | main (int argc, char *argv[]) |
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 34 of file mbIntxCheck.cpp.
35 {
36 std::stringstream sstr;
37 // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
38 const IntxAreaUtils::AreaMethod areaMethod = IntxAreaUtils::GaussQuadrature;
39
40 int rank = 0, size = 1;
41 #ifdef MOAB_HAVE_MPI
42 MPI_Init( &argc, &argv );
43 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
44 MPI_Comm_size( MPI_COMM_WORLD, &size );
45 #endif
46
47 std::string sourceFile, targetFile, intxFile;
48 std::string source_verif( "outS.h5m" ), target_verif( "outt.h5m" );
49 int sphere = 1;
50 int oldNamesParents = 0;
51 double areaErrSource = -1;
52 double areaErrTarget = -1;
53 ProgOptions opts;
54
55 opts.addOpt< std::string >( "source,s", "source file ", &sourceFile );
56 opts.addOpt< std::string >( "target,t", "target file ", &targetFile );
57 opts.addOpt< std::string >( "intersection,i", "intersection file ", &intxFile );
58 opts.addOpt< std::string >( "verif_source,v", "output source verification ", &source_verif );
59 opts.addOpt< std::string >( "verif_target,w", "output target verification ", &target_verif );
60 opts.addOpt< double >( "threshold_source,m", "error source threshold ", &areaErrSource );
61 opts.addOpt< double >( "threshold_target,q", "error target threshold ", &areaErrTarget );
62
63 opts.addOpt< int >( "sphere,p", "mesh on a sphere", &sphere );
64 opts.addOpt< int >( "old_convention,n", "old names for parent tags", &oldNamesParents );
65
66 opts.parseCommandLine( argc, argv );
67 // load meshes in parallel if needed
68 std::string opts_read = ( size == 1 ? ""
69 : std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
70 std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) );
71
72 // read meshes in 3 file sets
73 ErrorCode rval;
74 Core moab;
75 Interface* mb = &moab; // global
76 EntityHandle sset, tset, ixset;
77
78 // create meshsets and load files
79 rval = mb->create_meshset( MESHSET_SET, sset );MB_CHK_ERR( rval );
80 rval = mb->create_meshset( MESHSET_SET, tset );MB_CHK_ERR( rval );
81 rval = mb->create_meshset( MESHSET_SET, ixset );MB_CHK_ERR( rval );
82 if( 0 == rank ) std::cout << "Loading source file " << sourceFile << "\n";
83 rval = mb->load_file( sourceFile.c_str(), &sset, opts_read.c_str() );MB_CHK_SET_ERR( rval, "failed reading source file" );
84 if( 0 == rank ) std::cout << "Loading target file " << targetFile << "\n";
85 rval = mb->load_file( targetFile.c_str(), &tset, opts_read.c_str() );MB_CHK_SET_ERR( rval, "failed reading target file" );
86
87 if( 0 == rank ) std::cout << "Loading intersection file " << intxFile << "\n";
88 rval = mb->load_file( intxFile.c_str(), &ixset, opts_read.c_str() );MB_CHK_SET_ERR( rval, "failed reading intersection file" );
89 double R = 1.;
90 if( sphere )
91 {
92 // fix radius of both meshes, to be consistent with radius 1
93 rval = moab::IntxUtils::ScaleToRadius( mb, sset, R );MB_CHK_ERR( rval );
94 rval = moab::IntxUtils::ScaleToRadius( mb, tset, R );MB_CHK_ERR( rval );
95 }
96 Range intxCells;
97 rval = mb->get_entities_by_dimension( ixset, 2, intxCells );MB_CHK_ERR( rval );
98
99 Range sourceCells;
100 rval = mb->get_entities_by_dimension( sset, 2, sourceCells );MB_CHK_ERR( rval );
101
102 Range targetCells;
103 rval = mb->get_entities_by_dimension( tset, 2, targetCells );MB_CHK_ERR( rval );
104
105 Tag sourceParentTag;
106 Tag targetParentTag;
107 if( oldNamesParents )
108 {
109 rval = mb->tag_get_handle( "RedParent", targetParentTag );MB_CHK_SET_ERR( rval, "can't find target parent tag" );
110 rval = mb->tag_get_handle( "BlueParent", sourceParentTag );MB_CHK_SET_ERR( rval, "can't find source parent tag" );
111 }
112 else
113 {
114 rval = mb->tag_get_handle( "TargetParent", targetParentTag );MB_CHK_SET_ERR( rval, "can't find target parent tag" );
115 rval = mb->tag_get_handle( "SourceParent", sourceParentTag );MB_CHK_SET_ERR( rval, "can't find source parent tag" );
116 }
117
118 // error sets, for better visualization
119 EntityHandle errorSourceSet;
120 rval = mb->create_meshset( MESHSET_SET, errorSourceSet );MB_CHK_ERR( rval );
121 EntityHandle errorTargetSet;
122 rval = mb->create_meshset( MESHSET_SET, errorTargetSet );MB_CHK_ERR( rval );
123
124 std::map< int, double > sourceAreas;
125 std::map< int, double > targetAreas;
126
127 std::map< int, double > sourceAreasIntx;
128 std::map< int, double > targetAreasIntx;
129
130 std::map< int, int > sourceNbIntx;
131 std::map< int, int > targetNbIntx;
132
133 std::map< int, EntityHandle > sourceMap;
134 std::map< int, EntityHandle > targetMap;
135
136 Tag gidTag = mb->globalId_tag();
137
138 Tag areaTag;
139 rval = mb->tag_get_handle( "OrigArea", 1, MB_TYPE_DOUBLE, areaTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
140
141 moab::IntxAreaUtils areaAdaptor( areaMethod );
142 Range non_convex_intx_cells;
143 for( Range::iterator eit = sourceCells.begin(); eit != sourceCells.end(); ++eit )
144 {
145 EntityHandle cell = *eit;
146 const EntityHandle* verts;
147 int num_nodes;
148 rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_ERR( rval );
149 if( MB_SUCCESS != rval ) return -1;
150 std::vector< double > coords( 3 * num_nodes );
151 // get coordinates
152 rval = mb->get_coords( verts, num_nodes, &coords[0] );
153 if( MB_SUCCESS != rval ) return -1;
154 double area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, R );
155 int sourceID;
156 rval = mb->tag_get_data( gidTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
157 sourceAreas[sourceID] = area;
158 rval = mb->tag_set_data( areaTag, &cell, 1, &area );MB_CHK_ERR( rval );
159 sourceMap[sourceID] = cell;
160 }
161 for( Range::iterator eit = targetCells.begin(); eit != targetCells.end(); ++eit )
162 {
163 EntityHandle cell = *eit;
164 const EntityHandle* verts;
165 int num_nodes;
166 rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_ERR( rval );
167 if( MB_SUCCESS != rval ) return -1;
168 std::vector< double > coords( 3 * num_nodes );
169 // get coordinates
170 rval = mb->get_coords( verts, num_nodes, &coords[0] );
171 if( MB_SUCCESS != rval ) return -1;
172 double area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, R );
173 int targetID;
174 rval = mb->tag_get_data( gidTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
175 targetAreas[targetID] = area;
176 rval = mb->tag_set_data( areaTag, &cell, 1, &area );MB_CHK_ERR( rval );
177 targetMap[targetID] = cell;
178 }
179
180 for( Range::iterator eit = intxCells.begin(); eit != intxCells.end(); ++eit )
181 {
182 EntityHandle cell = *eit;
183 const EntityHandle* verts;
184 int num_nodes;
185 rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_ERR( rval );
186 if( MB_SUCCESS != rval ) return -1;
187 std::vector< double > coords( 3 * num_nodes );
188 // get coordinates
189 rval = mb->get_coords( verts, num_nodes, &coords[0] );
190 if( MB_SUCCESS != rval ) return -1;
191 int check_sign = 1;
192 double intx_area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, R, &check_sign );
193
194 rval = mb->tag_set_data( areaTag, &cell, 1, &intx_area );
195 ;MB_CHK_ERR( rval );
196 int sourceID, targetID;
197 rval = mb->tag_get_data( sourceParentTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
198 rval = mb->tag_get_data( targetParentTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
199
200 std::map< int, double >::iterator sit = sourceAreasIntx.find( sourceID );
201 if( sit == sourceAreasIntx.end() )
202 {
203 sourceAreasIntx[sourceID] = intx_area;
204 sourceNbIntx[sourceID] = 1;
205 }
206 else
207 {
208 sourceAreasIntx[sourceID] += intx_area;
209 sourceNbIntx[sourceID]++;
210 }
211
212 std::map< int, double >::iterator tit = targetAreasIntx.find( targetID );
213 if( tit == targetAreasIntx.end() )
214 {
215 targetAreasIntx[targetID] = intx_area;
216 targetNbIntx[targetID] = 1;
217 }
218 else
219 {
220 targetAreasIntx[targetID] += intx_area;
221 targetNbIntx[targetID]++;
222 }
223 if( intx_area < 0 )
224 {
225 std::cout << "negative intx area: " << intx_area << " n:" << num_nodes << " parents: " << sourceID << " "
226 << targetID << "\n";
227 }
228 if( check_sign < 0 )
229 {
230 non_convex_intx_cells.insert( cell );
231 std::cout << " non-convex polygon: " << intx_area << " n:" << num_nodes << " parents: " << sourceID << " "
232 << targetID << "\n";
233 }
234 }
235 Tag diffTag;
236 rval = mb->tag_get_handle( "AreaDiff", 1, MB_TYPE_DOUBLE, diffTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
237
238 Tag countIntxCellsTag;
239 rval = mb->tag_get_handle( "CountIntx", 1, MB_TYPE_INTEGER, countIntxCellsTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
240
241 for( Range::iterator eit = sourceCells.begin(); eit != sourceCells.end(); ++eit )
242 {
243 EntityHandle cell = *eit;
244
245 int sourceID;
246 rval = mb->tag_get_data( gidTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
247 double areaDiff = sourceAreas[sourceID];
248 std::map< int, double >::iterator sit = sourceAreasIntx.find( sourceID );
249 int countIntxCells = 0;
250 if( sit != sourceAreasIntx.end() )
251 {
252 areaDiff -= sourceAreasIntx[sourceID];
253 countIntxCells = sourceNbIntx[sourceID];
254 }
255 rval = mb->tag_set_data( diffTag, &cell, 1, &areaDiff );
256 rval = mb->tag_set_data( countIntxCellsTag, &cell, 1, &countIntxCells );
257 // add to errorSourceSet set if needed
258 if( ( areaErrSource > 0 ) && ( fabs( areaDiff ) > areaErrSource ) )
259 {
260 rval = mb->add_entities( errorSourceSet, &cell, 1 );MB_CHK_ERR( rval );
261 }
262 }
263 if( 0 == rank ) std::cout << "write source verification file " << source_verif << "\n";
264 rval = mb->write_file( source_verif.c_str(), 0, 0, &sset, 1 );MB_CHK_ERR( rval );
265 if( areaErrSource > 0 )
266 {
267 Range sourceErrorCells;
268 rval = mb->get_entities_by_handle( errorSourceSet, sourceErrorCells );MB_CHK_ERR( rval );
269 EntityHandle errorSourceIntxSet;
270 rval = mb->create_meshset( MESHSET_SET, errorSourceIntxSet );MB_CHK_ERR( rval );
271 if( !sourceErrorCells.empty() )
272 {
273 // add the intx cells that have these as source parent
274 std::vector< int > sourceIDs;
275 sourceIDs.resize( sourceErrorCells.size() );
276 rval = mb->tag_get_data( gidTag, sourceErrorCells, &sourceIDs[0] );MB_CHK_SET_ERR( rval, "can't get source IDs" );
277 std::sort( sourceIDs.begin(), sourceIDs.end() );
278 for( Range::iterator eit = intxCells.begin(); eit != intxCells.end(); ++eit )
279 {
280 EntityHandle cell = *eit;
281 int sourceID;
282 rval = mb->tag_get_data( sourceParentTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
283 std::vector< int >::iterator j = std::lower_bound( sourceIDs.begin(), sourceIDs.end(), sourceID );
284 if( ( j != sourceIDs.end() ) && ( *j == sourceID ) )
285 {
286 rval = mb->add_entities( errorSourceIntxSet, &cell, 1 );
287 ;MB_CHK_ERR( rval );
288 }
289 }
290 std::string filtersource = std::string( "filt_" ) + source_verif;
291 rval = mb->write_file( filtersource.c_str(), 0, 0, &errorSourceSet, 1 );
292 std::string filterIntxsource = std::string( "filtIntx_" ) + source_verif;
293 rval = mb->write_file( filterIntxsource.c_str(), 0, 0, &errorSourceIntxSet, 1 );
294 }
295 }
296
297 for( Range::iterator eit = targetCells.begin(); eit != targetCells.end(); ++eit )
298 {
299 EntityHandle cell = *eit;
300
301 int targetID;
302 rval = mb->tag_get_data( gidTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
303 double areaDiff = targetAreas[targetID];
304 int countIntxCells = 0;
305 std::map< int, double >::iterator sit = targetAreasIntx.find( targetID );
306 if( sit != targetAreasIntx.end() )
307 {
308 areaDiff -= targetAreasIntx[targetID];
309 countIntxCells = targetNbIntx[targetID];
310 }
311
312 rval = mb->tag_set_data( diffTag, &cell, 1, &areaDiff );
313 rval = mb->tag_set_data( countIntxCellsTag, &cell, 1, &countIntxCells );
314 // add to errorTargetSet set if needed
315 if( ( areaErrTarget > 0 ) && ( fabs( areaDiff ) > areaErrTarget ) )
316 {
317 rval = mb->add_entities( errorTargetSet, &cell, 1 );MB_CHK_ERR( rval );
318 }
319 }
320 if( 0 == rank ) std::cout << "write target verification file " << target_verif << "\n";
321 rval = mb->write_file( target_verif.c_str(), 0, 0, &tset, 1 );MB_CHK_ERR( rval );
322 if( areaErrTarget > 0 )
323 {
324 Range targetErrorCells;
325 rval = mb->get_entities_by_handle( errorTargetSet, targetErrorCells );MB_CHK_ERR( rval );
326 if( !targetErrorCells.empty() )
327 {
328 EntityHandle errorTargetIntxSet;
329 rval = mb->create_meshset( MESHSET_SET, errorTargetIntxSet );MB_CHK_ERR( rval );
330 // add the intx cells that have these as target parent
331 std::vector< int > targetIDs;
332 targetIDs.resize( targetErrorCells.size() );
333 rval = mb->tag_get_data( gidTag, targetErrorCells, &targetIDs[0] );MB_CHK_SET_ERR( rval, "can't get target IDs" );
334 std::sort( targetIDs.begin(), targetIDs.end() );
335 for( Range::iterator eit = intxCells.begin(); eit != intxCells.end(); ++eit )
336 {
337 EntityHandle cell = *eit;
338 int targetID;
339 rval = mb->tag_get_data( targetParentTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
340 std::vector< int >::iterator j = std::lower_bound( targetIDs.begin(), targetIDs.end(), targetID );
341 if( ( j != targetIDs.end() ) && ( *j == targetID ) )
342 {
343 rval = mb->add_entities( errorTargetIntxSet, &cell, 1 );
344 ;MB_CHK_ERR( rval );
345 }
346 }
347 std::string filterTarget = std::string( "filt_" ) + target_verif;
348 rval = mb->write_file( filterTarget.c_str(), 0, 0, &errorTargetSet, 1 );
349 std::string filterIntxtarget = std::string( "filtIntx_" ) + target_verif;
350 rval = mb->write_file( filterIntxtarget.c_str(), 0, 0, &errorTargetIntxSet, 1 );
351 }
352 }
353
354 if( !non_convex_intx_cells.empty() )
355 {
356
357 sourceCells.clear();
358 targetCells.clear();
359 for( Range::iterator it = non_convex_intx_cells.begin(); it != non_convex_intx_cells.end(); it++ )
360 {
361 EntityHandle cellIntx = *it;
362 int sourceID, targetID;
363 rval = mb->tag_get_data( sourceParentTag, &cellIntx, 1, &sourceID );MB_CHK_ERR( rval );
364 rval = mb->tag_get_data( targetParentTag, &cellIntx, 1, &targetID );MB_CHK_ERR( rval );
365 sourceCells.insert( sourceMap[sourceID] );
366 targetCells.insert( targetMap[targetID] );
367 }
368 EntityHandle nonConvexSet;
369 rval = mb->create_meshset( MESHSET_SET, nonConvexSet );MB_CHK_ERR( rval );
370 rval = mb->add_entities( nonConvexSet, non_convex_intx_cells );
371 rval = mb->write_file( "nonConvex.h5m", 0, 0, &nonConvexSet, 1 );MB_CHK_ERR( rval );
372
373 EntityHandle sSet;
374 rval = mb->create_meshset( MESHSET_SET, sSet );MB_CHK_ERR( rval );
375 rval = mb->add_entities( sSet, sourceCells );
376 rval = mb->write_file( "nonConvexSource.h5m", 0, 0, &sSet, 1 );MB_CHK_ERR( rval );
377 EntityHandle tSet;
378 rval = mb->create_meshset( MESHSET_SET, tSet );MB_CHK_ERR( rval );
379 rval = mb->add_entities( tSet, targetCells );
380 rval = mb->write_file( "nonConvexTarget.h5m", 0, 0, &tSet, 1 );MB_CHK_ERR( rval );
381 rval = mb->add_entities( nonConvexSet, sourceCells );
382 rval = mb->add_entities( nonConvexSet, targetCells );
383 rval = mb->write_file( "nonConvexAll.h5m", 0, 0, &nonConvexSet, 1 );MB_CHK_ERR( rval );
384 }
385 return 0;
386 }
References moab::Core::add_entities(), ProgOptions::addOpt(), moab::IntxAreaUtils::area_spherical_polygon(), moab::Range::begin(), moab::Range::clear(), moab::Core::create_meshset(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::IntxAreaUtils::GaussQuadrature, moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::Core::get_entities_by_handle(), moab::Core::globalId_tag(), moab::Range::insert(), moab::Core::load_file(), mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MESHSET_SET, ProgOptions::parseCommandLine(), moab::R, moab::IntxUtils::ScaleToRadius(), size, moab::Range::size(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), and moab::Core::write_file().
1.9.1.