1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367 | //
// Usage:
// tools/mbslavepart -d 2 -m mpas/x1.2562.grid.h5m -s mpas/x1.10242.grid.h5m -o mpas_slave.h5m -e
// 1e-8 -b 1e-6 -O
//
#include <iostream>
#include <exception>
#include <cmath>
#include <vector>
#include <string>
#include "moab/ProgOptions.hpp"
#include "moab/Core.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/BVHTree.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab_mpi.h"
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#endif
using namespace moab;
// A function to get the non-default value from a std::map
template < typename K, typename V >
static V get_map_value( const std::map< K, V >& m, const K& key, const V& defval )
{
typename std::map< K, V >::const_iterator it = m.find( key );
if( it == m.end() )
{
return defval;
}
else
{
return it->second;
}
}
int main( int argc, char* argv[] )
{
int nprocs = 1, dimension = 3;<--- 'nprocs' is assigned value '1' here.<--- 'nprocs' is assigned value '1' here.
#ifdef MOAB_HAVE_MPI
int proc_id = 0;
MPI_Init( &argc, &argv );
MPI_Comm_rank( MPI_COMM_WORLD, &proc_id );
MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
#endif
int defaultpart = 0;
double tolerance = 1e-2, treetolerance = 1e-13, btolerance = 1e-4;<--- The scope of the variable 'treetolerance' can be reduced. [+]The scope of the variable 'treetolerance' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
int i = 0;
if (x) {
// it's safe to move 'int i = 0;' here
for (int n = 0; n < 10; ++n) {
// it is possible but not safe to move 'int i = 0;' here
do_something(&i);
}
}
}
When you see this message it is always safe to reduce the variable scope 1 level.
std::string masterfile, slavefile, outfile( "slavemesh.h5m" );
bool keepsparts = false;
bool use_spherical = false;
ProgOptions opts;
opts.addOpt< std::string >( "master,m", "Master mesh filename", &masterfile );
opts.addOpt< std::string >( "slave,s", "Slave mesh filename", &slavefile );
opts.addOpt< std::string >( "output,o", "Output partitioned mesh filename", &outfile );
opts.addOpt< int >( "dim,d", "Dimension of entities to use for partitioning", &dimension );
opts.addOpt< int >( "defaultpart,p", "Default partition number if target element is not found on source grid",
&defaultpart );
opts.addOpt< double >( "eps,e", "Tolerance for the point search", &tolerance );
opts.addOpt< double >( "beps,b", "Tolerance for the bounding box search", &btolerance );
opts.addOpt< void >( "keep,K",
"Keep the existing partitions in the slave mesh (use PARALLEL_PARTITION_SLAVE instead)",
&keepsparts );
opts.addOpt< void >( "spherical", "Hint that the meshes are defined on a spherical surface (Climate problems)",
&use_spherical );
opts.parseCommandLine( argc, argv );
if( masterfile.empty() || slavefile.empty() )
{
opts.printHelp();
#ifdef MOAB_HAVE_MPI
MPI_Finalize();
#endif
exit( 1 );
}
ErrorCode error;
Core* mbCore = new Core();
// Set the read options for parallel file loading
const std::string partition_set_name = "PARALLEL_PARTITION";
const std::string global_id_name = "GLOBAL_ID";
std::vector< std::string > read_opts, write_opts;<--- Unused variable: read_opts<--- Unused variable: write_opts
std::string read_options( "" ), write_options( "" );
if( nprocs > 1 )<--- The comparison 'nprocs > 1' is always false. [+]Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
{
read_options = "PARALLEL=READ_PART;PARTITION=" + partition_set_name + ";PARALLEL_RESOLVE_SHARED_ENTS";
write_options = "PARALLEL=WRITE_PART";
}
EntityHandle masterfileset, slavefileset;
error = mbCore->create_meshset( moab::MESHSET_TRACK_OWNER | moab::MESHSET_SET, masterfileset );MB_CHK_ERR( error );
error = mbCore->create_meshset( moab::MESHSET_TRACK_OWNER | moab::MESHSET_SET, slavefileset );MB_CHK_ERR( error );
// Load file
error = mbCore->load_file( masterfile.c_str(), &masterfileset, read_options.c_str() );MB_CHK_ERR( error );
error = mbCore->load_file( slavefile.c_str(), &slavefileset, read_options.c_str() );MB_CHK_ERR( error );
// if (error != MB_SUCCESS && size > 1)
// {
// std::string newread_options = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS";
// error = mbCore->load_file(slavefile.c_str(), &slavefileset, newread_options.c_str());
// }
// else MB_CHK_ERR(error);
Tag gidtag = 0, parttag = 0, sparttag = 0;
int dum_id = -1;
error = mbCore->tag_get_handle( partition_set_name.c_str(), 1, MB_TYPE_INTEGER, parttag,
MB_TAG_SPARSE | MB_TAG_CREAT, &dum_id );MB_CHK_ERR( error );
gidtag = mbCore->globalId_tag();
if( keepsparts )
{
error = mbCore->tag_get_handle( std::string( partition_set_name + "_SLAVE" ).c_str(), 1, MB_TYPE_INTEGER,
sparttag, MB_TAG_CREAT | MB_TAG_SPARSE, &dum_id );MB_CHK_ERR( error );
}
Range melems, msets, selems, ssets;
// Get the partition sets on the master mesh
std::map< int, int > mpartvals;
error = mbCore->get_entities_by_type_and_tag( masterfileset, MBENTITYSET, &parttag, NULL, 1, msets,
moab::Interface::UNION, true );MB_CHK_ERR( error );
if( msets.size() == 0 )
{
std::cout << "No partition sets found in the master mesh. Quitting..." << std::endl;
exit( 1 );
}
for( unsigned i = 0; i < msets.size(); ++i )
{
EntityHandle mset = msets[i];
moab::Range msetelems;
error = mbCore->get_entities_by_dimension( mset, dimension, msetelems );MB_CHK_ERR( error );
melems.merge( msetelems );
int partID;
error = mbCore->tag_get_data( parttag, &mset, 1, &partID );MB_CHK_ERR( error );
// Get the global ID and use that as the indicator
std::vector< int > gidMelems( msetelems.size() );
error = mbCore->tag_get_data( gidtag, msetelems, gidMelems.data() );MB_CHK_ERR( error );
for( unsigned j = 0; j < msetelems.size(); ++j )
mpartvals[gidMelems[j]] = partID;
// mpartvals[msetelems[j]]=partID;
#ifdef VERBOSE
std::cout << "Part " << partID << " has " << msetelems.size() << " elements." << std::endl;
#endif
}
// Get information about the slave file set
error = mbCore->get_entities_by_type_and_tag( slavefileset, MBENTITYSET, &parttag, NULL, 1, ssets,
moab::Interface::UNION );MB_CHK_ERR( error );
// TODO: expand and add other dimensional elements
error = mbCore->get_entities_by_dimension( slavefileset, dimension, selems );MB_CHK_ERR( error );
std::cout << "Master (elements, parts) : (" << melems.size() << ", " << msets.size()
<< "), Slave (elements, parts) : (" << selems.size() << ", " << ssets.size() << ")" << std::endl;
double master_radius = 1.0, slave_radius = 1.0;<--- The scope of the variable 'master_radius' can be reduced. [+]The scope of the variable 'master_radius' can be reduced. Warning: Be careful when fixing this message, especially when there are inner loops. Here is an example where cppcheck will write that the scope for 'i' can be reduced:
void f(int x)
{
int i = 0;
if (x) {
// it's safe to move 'int i = 0;' here
for (int n = 0; n < 10; ++n) {
// it is possible but not safe to move 'int i = 0;' here
do_something(&i);
}
}
}
When you see this message it is always safe to reduce the variable scope 1 level.
std::vector< double > mastercoords;<--- Unused variable: mastercoords
Range masterverts, slaveverts;
{
error = mbCore->get_entities_by_dimension( masterfileset, 0, masterverts );MB_CHK_ERR( error );
error = mbCore->get_entities_by_dimension( slavefileset, 0, slaveverts );MB_CHK_ERR( error );
}
if( use_spherical )
{
double points[6];
EntityHandle mfrontback[2] = { masterverts[0], masterverts[masterverts.size() - 1] };
error = mbCore->get_coords( &mfrontback[0], 2, points );MB_CHK_ERR( error );
master_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
EntityHandle sfrontback[2] = { slaveverts[0], slaveverts[slaveverts.size() - 1] };
error = mbCore->get_coords( &sfrontback[0], 2, points );MB_CHK_ERR( error );
slave_radius = 0.5 * ( std::sqrt( points[0] * points[0] + points[1] * points[1] + points[2] * points[2] ) +
std::sqrt( points[3] * points[3] + points[4] * points[4] + points[5] * points[5] ) );
// Let us rescale both master and slave meshes to a unit sphere
error = moab::IntxUtils::ScaleToRadius( mbCore, masterfileset, 1.0 );MB_CHK_ERR( error );
error = moab::IntxUtils::ScaleToRadius( mbCore, slavefileset, 1.0 );MB_CHK_ERR( error );
}
try
{
std::map< int, moab::Range > spartvals;
int npoints_notfound = 0;
{
FileOptions fopts( ( use_spherical ? "SPHERICAL" : "" ) );
moab::EntityHandle tree_root = masterfileset;
moab::AdaptiveKDTree tree( mbCore );
error = tree.build_tree( melems, &tree_root, &fopts );
// moab::BVHTree tree(mbCore);
// error = tree.build_tree(melems, &tree_root);MB_CHK_ERR(error);
for( size_t ie = 0; ie < selems.size(); ie++ )
{
moab::EntityHandle selem, leaf;
double point[3];
selem = selems[ie];
// Get the element centroid to be queried
error = mbCore->get_coords( &selem, 1, point );MB_CHK_ERR( error );
std::vector< moab::EntityHandle > leaf_elems;
// Search for the closest source element in the master mesh corresponding
// to the target element centroid in the slave mesh
error = tree.point_search( point, leaf, treetolerance, btolerance );MB_CHK_ERR( error );
// We only care about the dimension that the user specified.
// MOAB partitions are ordered by elements anyway.
error = mbCore->get_entities_by_dimension( leaf, dimension, leaf_elems, true );MB_CHK_ERR( error );
if( leaf != 0 && leaf_elems.size() )
{
// Now get the master element centroids so that we can compute
// the minimum distance to the target point
std::vector< double > centroids( leaf_elems.size() * 3 );
error = mbCore->get_coords( &leaf_elems[0], leaf_elems.size(), ¢roids[0] );MB_CHK_ERR( error );
if( !leaf_elems.size() )
std::cout << ie << ": "
<< " No leaf elements found." << std::endl;
// else std::cout << ie << " found " << leaf_elems.size() << " leaves for
// current point " << std::endl;
double dist = 1e10;
int pinelem = -1;
for( size_t il = 0; il < leaf_elems.size(); ++il )
{
const double* centroid = ¢roids[il * 3];
const double locdist = std::pow( point[0] - centroid[0], 2 ) +
std::pow( point[1] - centroid[1], 2 ) +
std::pow( point[2] - centroid[2], 2 );
if( locdist < dist && locdist < 1.0E-2 )
{
dist = locdist;
pinelem = il;
#ifdef VERBOSE
int gidMelem;
error = mbCore->tag_get_data( gidtag, &leaf_elems[il], 1, &gidMelem );MB_CHK_ERR( error );
std::cout << "\t Trial leaf " << il << " set " << gidMelem
<< " and part = " << get_map_value( mpartvals, gidMelem, -1 )
<< " with distance = " << locdist << std::endl;
#endif
}
}
if( pinelem < 0 )
{
#ifdef VERBOSE
std::cout << ie
<< ": [Error] - Could not find a minimum distance within the "
"leaf nodes."
<< std::endl;
#endif
npoints_notfound++;
spartvals[defaultpart].insert( selems[ie] );
}
else
{
int gidMelem;
error = mbCore->tag_get_data( gidtag, &leaf_elems[pinelem], 1, &gidMelem );MB_CHK_ERR( error );
int mpartid = get_map_value( mpartvals, gidMelem, -1 );
if( mpartid < 0 )
std::cout << "[WARNING]: Part number for element " << leaf_elems[pinelem]
<< " with global ID = " << gidMelem << " not found.\n";
#ifdef VERBOSE
std::cout << "Adding element " << ie << " set " << mpartid << " with distance = " << dist
<< std::endl;
#endif
spartvals[mpartid].insert( selems[ie] );
}
}
else
{
#ifdef VERBOSE
std::cout << "[WARNING]: Adding element " << ie << " to default (part) set " << defaultpart
<< std::endl;
#endif
spartvals[defaultpart].insert( selems[ie] );
}
}
error = tree.reset_tree();MB_CHK_ERR( error );
}
if( npoints_notfound )
std::cout << "Could not find " << npoints_notfound
<< " points in the master mesh. Added to defaultpart set = " << defaultpart << std::endl;
if( use_spherical )
{
error = moab::IntxUtils::ScaleToRadius( mbCore, slavefileset, slave_radius );MB_CHK_ERR( error );
}
error = mbCore->delete_entities( &masterfileset, 1 );MB_CHK_ERR( error );
// Find parallel partition sets in the slave mesh - and delete it since we are going to
// overwrite the sets
if( !keepsparts )
{
std::cout << "Deleting " << ssets.size() << " sets in the slave mesh" << std::endl;
error = mbCore->remove_entities( slavefileset, ssets );MB_CHK_ERR( error );
ssets.clear();
}
size_t ntotslave_elems = 0, ntotslave_parts = 0;
for( std::map< int, moab::Range >::iterator it = spartvals.begin(); it != spartvals.end(); ++it )
{
int partID = it->first;
moab::EntityHandle pset;
error = mbCore->create_meshset( moab::MESHSET_SET, pset );MB_CHK_ERR( error );
error = mbCore->add_entities( pset, it->second );MB_CHK_ERR( error );
error = mbCore->add_parent_child( slavefileset, pset );MB_CHK_ERR( error );
#ifdef VERBOSE
std::cout << "Slave Part " << partID << " has " << it->second.size() << " elements." << std::endl;
#endif
ntotslave_elems += it->second.size();
ntotslave_parts++;
if( keepsparts )
{
error = mbCore->tag_set_data( sparttag, &pset, 1, &partID );MB_CHK_ERR( error );
}
else
{
error = mbCore->tag_set_data( parttag, &pset, 1, &partID );MB_CHK_ERR( error );
}
}
std::cout << "Slave mesh: given " << selems.size() << " elements, and assigned " << ntotslave_elems
<< " elements to " << ntotslave_parts << " parts.\n";
assert( ntotslave_elems == selems.size() );
// mbCore->print_database();
// Write the re-partitioned slave mesh to disk
if( nprocs == 1 )<--- The comparison 'nprocs == 1' is always true. [+]Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
{
error = mbCore->write_file( "slavemesh.vtk", "VTK", NULL, &slavefileset, 1 );MB_CHK_ERR( error );
}
error = mbCore->write_file( outfile.c_str(), NULL, write_options.c_str(), &slavefileset, 1 );MB_CHK_ERR( error );
// error = mbCore->write_file(outfile.c_str(), NULL,
// write_options.c_str());MB_CHK_ERR(error);
}
catch( std::exception& e )
{
std::cout << " exception caught during tree initialization " << e.what() << std::endl;
}
delete mbCore;
#ifdef MOAB_HAVE_MPI
MPI_Finalize();
#endif
exit( 0 );
}
|