Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
mbIntxCheck.cpp File Reference
#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[])
 

Function Documentation

◆ main()

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().