Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
visuMapVtk.cpp
Go to the documentation of this file.
1 /*
2  * visuMapVtk.cpp
3  * this tool will take a source file, target file (h5m) and a map file in nc format,
4  *
5  * example of usage:
6  * ./mbvisumap -s source.h5m -t target.h5m -m map.nc -b startSourceID -e endSourceID -c startTargetID -f endTargetID
7  * will associate row i in map with a partial mesh
8  * will associate row j in map with a partial mesh
9  *
10  * can be built only if netcdf and hdf5 and eigen3 are available
11  *
12  *
13  */
14 #include "moab/MOABConfig.h"
15 
16 #ifndef MOAB_HAVE_EIGEN3
17 #error compareMaps tool requires eigen3 configuration
18 #endif
19 
20 #include "moab/ProgOptions.hpp"
21 #include "moab/Core.hpp"
22 #include "moab/Range.hpp"
23 
24 #include "netcdf.h"
25 #include <cmath>
26 #include <sstream>
27 #include <map>
28 #include <Eigen/Sparse>
29 
30 #define ERR_NC( e ) \
31  { \
32  printf( "Error: %s\n", nc_strerror( e ) ); \
33  exit( 2 ); \
34  }
35 
36 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
37 // ncFile is an integer initialized when opening the nc file in read mode
38 
39 int ncFile1;
40 
41 #define GET_DIM1( ncdim, name, val ) \
42  { \
43  int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
44  if( NC_NOERR == gdfail ) \
45  { \
46  size_t tmp_val; \
47  gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
48  if( NC_NOERR != gdfail ) \
49  { \
50  ERR_NC( gdfail ) \
51  } \
52  else \
53  ( val ) = tmp_val; \
54  } \
55  else \
56  ( val ) = 0; \
57  }
58 
59 #define GET_VAR1( name, id, dims ) \
60  { \
61  ( id ) = -1; \
62  int gvfail = nc_inq_varid( ncFile1, name, &( id ) ); \
63  if( NC_NOERR == gvfail ) \
64  { \
65  int ndims; \
66  gvfail = nc_inq_varndims( ncFile1, id, &ndims ); \
67  if( NC_NOERR == gvfail ) \
68  { \
69  ( dims ).resize( ndims ); \
70  gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
71  if( NC_NOERR != gvfail ) \
72  { \
73  ERR_NC( gvfail ) \
74  } \
75  } \
76  } \
77  }
78 
79 #define GET_1D_INT_VAR1( name, id, vals ) \
80  { \
81  GET_VAR1( name, id, vals ); \
82  if( -1 != ( id ) ) \
83  { \
84  size_t ntmp; \
85  int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp ); \
86  if( NC_NOERR != ivfail ) \
87  { \
88  ERR_NC( ivfail ) \
89  } \
90  ( vals ).resize( ntmp ); \
91  size_t ntmp1 = 0; \
92  ivfail = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
93  if( NC_NOERR != ivfail ) \
94  { \
95  ERR_NC( ivfail ) \
96  } \
97  } \
98  }
99 
100 #define GET_1D_DBL_VAR1( name, id, vals ) \
101  { \
102  std::vector< int > dum_dims; \
103  GET_VAR1( name, id, dum_dims ); \
104  if( -1 != ( id ) ) \
105  { \
106  size_t ntmp; \
107  int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp ); \
108  if( NC_NOERR != dvfail ) \
109  { \
110  ERR_NC( dvfail ) \
111  } \
112  ( vals ).resize( ntmp ); \
113  size_t ntmp1 = 0; \
114  dvfail = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
115  if( NC_NOERR != dvfail ) \
116  { \
117  ERR_NC( dvfail ) \
118  } \
119  } \
120  }
121 
122 using namespace moab;
123 using namespace std;
124 
125 int main( int argc, char* argv[] )
126 {
127 
128  ProgOptions opts;
129  int dimSource = 2; // for FV meshes is 2; for SE meshes, use fine mesh, dim will be 0
130  int dimTarget = 2; //
131  int otype = 0;
132 
133  std::string inputfile1, inputSource, inputTarget;
134  opts.addOpt< std::string >( "map,m", "input map ", &inputfile1 );
135  opts.addOpt< std::string >( "source,s", "source mesh", &inputSource );
136  opts.addOpt< std::string >( "target,t", "target mesh", &inputTarget );
137  opts.addOpt< int >( "dimSource,d", "dimension of source ", &dimSource );
138  opts.addOpt< int >( "dimTarget,g", "dimension of target ", &dimTarget );
139  opts.addOpt< int >( "typeOutput,o", " output type vtk(0), h5m(1)", &otype );
140 
141  int startSourceID = -1, endSourceID = -1, startTargetID = -1, endTargetID = -1;
142  opts.addOpt< int >( "startSourceID,b", "start source id ", &startSourceID );
143  opts.addOpt< int >( "endSourceID,e", "end source id ", &endSourceID );
144  opts.addOpt< int >( "startTargetID,c", "start target id ", &startTargetID );
145  opts.addOpt< int >( "endTargetID,f", "end target id ", &endTargetID );
146  // -b startSourceID -e endSourceID -c startTargetID -f endTargetID
147 
148  opts.parseCommandLine( argc, argv );
149 
150  std::string extension = ".vtk";
151  if( 1 == otype ) extension = ".h5m";
152  // Open netcdf/exodus file
153  int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
154  if( NC_NOWRITE != fail )
155  {
156  ERR_NC( fail )
157  }
158 
159  std::cout << " opened " << inputfile1 << " for map 1 \n";
160 
161  int temp_dim;
162  int na1, nb1, ns1;
163  GET_DIM1( temp_dim, "n_a", na1 );
164  GET_DIM1( temp_dim, "n_b", nb1 );
165  GET_DIM1( temp_dim, "n_s", ns1 );
166  std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
167  std::vector< int > col1( ns1 ), row1( ns1 );
168  std::vector< double > val1( ns1 );
169  int idrow1, idcol1, ids1;
170  GET_1D_INT_VAR1( "row", idrow1, row1 );
171  GET_1D_INT_VAR1( "col", idcol1, col1 );
172  GET_1D_DBL_VAR1( "S", ids1, val1 );
173  // first matrix
174  typedef Eigen::Triplet< double > Triplet;
175  std::vector< Triplet > tripletList;
176  tripletList.reserve( ns1 );
177  for( int iv = 0; iv < ns1; iv++ )
178  {
179  // all row and column indices are 1-based in the map file.
180  // inside Eigen3, we will use 0-based; then we will have to add back 1 when
181  // we dump out the files
182  tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
183  }
184  Eigen::SparseMatrix< double > weight1( nb1, na1 );
185 
186  weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
187  weight1.makeCompressed();
188  // we read the matrix; now read moab source and target
189  Core core;
190  Interface* mb = &core;
191  ErrorCode rval;
192  Tag gtag = mb->globalId_tag();
193  EntityHandle sourceSet, targetSet;
194  // those are the maps from global ids to the moab entity handles corresponding to those global ids
195  // which are corresponding to the global DOFs
196  map< int, EntityHandle > sourceHandles;
197  map< int, EntityHandle > targetHandles;
198  rval = mb->create_meshset( MESHSET_SET, sourceSet );MB_CHK_SET_ERR( rval, "can't create source mesh set" );
199  rval = mb->create_meshset( MESHSET_SET, targetSet );MB_CHK_SET_ERR( rval, "can't create target mesh set" );
200  const char* readopts = "";
201  rval = mb->load_file( inputSource.c_str(), &sourceSet, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
202  rval = mb->load_file( inputTarget.c_str(), &targetSet, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
203  Range sRange;
204  rval = mb->get_entities_by_dimension( sourceSet, dimSource, sRange );MB_CHK_SET_ERR( rval, "Failed to get sRange" );
205  vector< int > sids;
206  sids.resize( sRange.size() );
207  rval = mb->tag_get_data( gtag, sRange, &sids[0] );MB_CHK_SET_ERR( rval, "Failed to get ids for srange" );
208  // all global ids are 1 based in the source file, and they correspond to the dofs in the map file
209  for( size_t i = 0; i < sids.size(); i++ )
210  {
211  EntityHandle eh = sRange[i];
212  int gid = sids[i];
213  sourceHandles[gid] = eh;
214  }
215  Range tRange;
216  rval = mb->get_entities_by_dimension( targetSet, dimTarget, tRange );MB_CHK_SET_ERR( rval, "Failed to get tRange" );
217  vector< int > tids;
218  tids.resize( tRange.size() );
219  rval = mb->tag_get_data( gtag, tRange, &tids[0] );MB_CHK_SET_ERR( rval, "Failed to get ids for trange" );
220  // all global ids are 1 based in the target file, and they correspond to the dofs in the map file
221  for( size_t i = 0; i < tids.size(); i++ )
222  {
223  EntityHandle eh = tRange[i];
224  int gid = tids[i];
225  targetHandles[gid] = eh;
226  }
227  // a dense tag for weights
228  Tag wtag;
229  double defVal = 0;
230 
231  std::string name_map = inputfile1;
232  // strip last 3 chars (.nc extension)
233  name_map.erase( name_map.begin() + name_map.length() - 3, name_map.end() );
234  // if path , remove from name
235  size_t pos = name_map.rfind( '/', name_map.length() );
236  if( pos != std::string::npos ) name_map = name_map.erase( 0, pos + 1 );
237 
238  rval = mb->tag_get_handle( "weight", 1, MB_TYPE_DOUBLE, wtag, MB_TAG_CREAT | MB_TAG_DENSE, &defVal );MB_CHK_SET_ERR( rval, "Failed to create weight" );
239  EntityHandle partialSet;
240  rval = mb->create_meshset( MESHSET_SET, partialSet );MB_CHK_SET_ERR( rval, "can't create partial set" );
241  // how to get a complete row in sparse matrix? Or a complete column ?
242  for( int col = startSourceID - 1; col <= endSourceID - 1; col++ )
243  {
244  Range targetEnts; // will find entries for column col-1 in sparse matrix weight1, and its entries
245  // will assign a dense tag with values, and write out the file
246  if( col < 0 ) continue;
247  Eigen::SparseVector< double > colVect = weight1.col( col );
248  // the row indices correspond to target cells
249  for( Eigen::SparseVector< double >::InnerIterator it( colVect ); it; ++it )
250  {
251  double weight = it.value(); // == vec[ it.index() ]
252  // we add back the 1 that we subtract
253  int globalIdRow = it.index() + 1;
254  EntityHandle th = targetHandles[globalIdRow];
255  targetEnts.insert( th );
256  rval = mb->tag_set_data( wtag, &th, 1, &weight );MB_CHK_SET_ERR( rval, "Failed to set weight tag on target" );
257  }
258 
259  if( dimTarget == 0 )
260  {
261  Range adjCells;
262  rval = mb->get_adjacencies( targetEnts, 2, false, adjCells, Interface::UNION );MB_CHK_SET_ERR( rval, " can't get adj cells " );
263  targetEnts.merge( adjCells );
264  }
265 
266  rval = mb->add_entities( partialSet, targetEnts );MB_CHK_SET_ERR( rval, "Failed to add target entities to partial set" );
267  // write now the set in a numbered file
268  std::stringstream fff;
269  fff << name_map << "_column" << col + 1 << extension;
270  rval = mb->write_mesh( fff.str().c_str(), &partialSet, 1 );MB_CHK_ERR( rval );
271  // remove from partial set the entities it has
272  rval = mb->clear_meshset( &partialSet, 1 );MB_CHK_SET_ERR( rval, "Failed to empty partial set" );
273  }
274 
275  // how to get a complete row in sparse matrix?
276  for( int row = startTargetID - 1; row <= endTargetID - 1; row++ )
277  {
278  Range sourceEnts; // will find entries for row in sparse matrix weight1, and its entries
279  // will assign a dense tag with values, and write out the file
280  if( row < 0 ) continue;
281  Eigen::SparseVector< double > rowVect = weight1.row( row );
282  // the row indices correspond to target cells
283  for( Eigen::SparseVector< double >::InnerIterator it( rowVect ); it; ++it )
284  {
285  double weight = it.value(); // == vec[ it.index() ]
286  int globalIdCol = it.index() + 1;
287  EntityHandle sh = sourceHandles[globalIdCol];
288  sourceEnts.insert( sh );
289  rval = mb->tag_set_data( wtag, &sh, 1, &weight );MB_CHK_SET_ERR( rval, "Failed to set weight tag on source" );
290  }
291  if( dimSource == 0 )
292  {
293  Range adjCells;
294  rval = mb->get_adjacencies( sourceEnts, 2, false, adjCells, Interface::UNION );MB_CHK_SET_ERR( rval, " can't get adj cells " );
295  sourceEnts.merge( adjCells );
296  }
297  rval = mb->add_entities( partialSet, sourceEnts );MB_CHK_SET_ERR( rval, "Failed to add source entities" );
298  // write now the set in a numbered file
299  std::stringstream fff;
300  fff << name_map << "_row" << row + 1 << extension;
301  rval = mb->write_mesh( fff.str().c_str(), &partialSet, 1 );MB_CHK_ERR( rval );
302  rval = mb->clear_meshset( &partialSet, 1 );MB_CHK_SET_ERR( rval, "Failed to empty partial set" );
303  }
304  return 0;
305 }