Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
addConnec.cpp
Go to the documentation of this file.
1 /*
2  *
3  *
4 /** @example addConnec.cpp Add connectivity to point cloud data, for better view in visit
5  * this tool will take an existing h5m fine point cloud data (phys grid or land pc)
6  * and add 2d cells from a fine atm mesh
7  *
8  * example of usage:
9  * ./addConnec -i wholeFineATM.h5m -s wholeLND_proj01.h5m -o wholeLndM.h5m
10  *
11  * Basically, will output a new h5m file (wholeLndM.h5m), which has also cell connectivity copied from fine atm
12  * matching is based on the global ids
13  * between what we think is the order on the original file (wholeFineATM.h5m)
14  *
15  * file wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4,
16  *
17  */
18 #include "moab/ProgOptions.hpp"
19 #include "moab/Core.hpp"
20 
21 #include <cmath>
22 #include <sstream>
23 
24 using namespace moab;
25 
26 int main( int argc, char* argv[] )
27 {
28 
29  ProgOptions opts;
30 
31  std::string inputfile, outfile( "out.h5m" ), sourcefile;
32 
33  opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
34  opts.addOpt< std::string >( "source,s", "h5m file aligned with the mesh input file", &sourcefile );
35  opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
36 
37  opts.parseCommandLine( argc, argv );
38 
39  std::cout << "input mesh file: " << inputfile << "\n";
40  std::cout << "source file: " << sourcefile << "\n";
41  std::cout << "output file: " << outfile << "\n";
42 
43  if( inputfile.empty() )
44  {
45  opts.printHelp();
46  return 0;
47  }
48  ErrorCode rval;
49  Core* mb = new Core();
50 
51  rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );
52 
53  Core* mb2 = new Core();
54  rval = mb2->load_file( sourcefile.c_str() );MB_CHK_SET_ERR( rval, "can't load source file" );
55 
56  // get vertices on ini mesh; get global id on ini mesh
57  // get global id on source mesh
58  Tag gid;
59  Tag gid2;
60  rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on ini mesh " );
61  rval = mb2->tag_get_handle( "GLOBAL_ID", gid2 );MB_CHK_SET_ERR( rval, "can't get GLOBAL_ID tag on source mesh " );
62 
63  // get vertices on ini mesh; build
64  Range iniVerts;
65  rval = mb->get_entities_by_dimension( 0, 0, iniVerts );MB_CHK_SET_ERR( rval, "can't get verts on initial mesh " );
66 
67  std::vector< int > gids;
68  gids.resize( iniVerts.size() );
69  rval = mb->tag_get_data( gid, iniVerts, &( gids[0] ) );MB_CHK_SET_ERR( rval, "can't get gid on initial verts " );
70  // build now the map
71  std::map< int, EntityHandle > fromGidToEh;
72  int i = 0;
73  for( Range::iterator vit = iniVerts.begin(); vit != iniVerts.end(); ++vit, i++ )
74  {
75  fromGidToEh[gids[i]] = *vit;
76  }
77  // now get the source verts, and tags, and set it on new mesh
78 
79  Range sourceVerts;
80  rval = mb2->get_entities_by_dimension( 0, 0, sourceVerts );MB_CHK_SET_ERR( rval, "can't get verts on source mesh " );
81  std::vector< int > gids2;
82  gids2.resize( sourceVerts.size() );
83  rval = mb2->tag_get_data( gid2, sourceVerts, &( gids2[0] ) );MB_CHK_SET_ERR( rval, "can't get gid2 on cloud mesh " );
84  std::map< int, EntityHandle > fromGid2ToEh;
85  i = 0;
86  for( Range::iterator vit = sourceVerts.begin(); vit != sourceVerts.end(); ++vit, i++ )
87  {
88  fromGid2ToEh[gids2[i]] = *vit;
89  }
90 
91  Range usedVerts;
92  for( size_t i = 0; i < gids2.size(); i++ )
93  {
94  int globalId = gids2[i];
95  EntityHandle usedVert = fromGidToEh[globalId];
96  usedVerts.insert( usedVert );
97  }
98  Range cells;
99  rval = mb->get_adjacencies( usedVerts, 2, false, cells, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adj cells " );
100 
101  // now create a new cell in mb2 for each one in mb
102  for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
103  {
104  EntityHandle cell = *cit;
105  int nnodes = 0;
106  const EntityHandle* conn = NULL;
107  rval = mb->get_connectivity( cell, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity " );
108 
109  std::vector< EntityHandle > nconn( nnodes );
110  bool goodCell = true;
111  for( int i = 0; i < nnodes; i++ )
112  {
113  EntityHandle v = conn[i];
114  int index = iniVerts.index( v );
115  if( index < 0 ) goodCell = false;
116  int id = gids[index];
117  if( fromGid2ToEh.find( id ) == fromGid2ToEh.end() )
118  goodCell = false;
119  else
120  nconn[i] = fromGid2ToEh[id];
121  }
122  if( !goodCell ) continue;
123  EntityType type = MBTRI;
124  if( nnodes == 3 )
125  type = MBTRI;
126  else if( nnodes == 4 )
127  type = MBQUAD;
128  else if( nnodes > 4 )
129  type = MBPOLYGON;
130  EntityHandle nel;
131  rval = mb2->create_element( type, &nconn[0], nnodes, nel );MB_CHK_SET_ERR( rval, "can't create new cell " );
132  }
133  // save file
134  rval = mb2->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
135 
136  delete mb;
137  delete mb2;
138 
139  return 0;
140 }