Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
visuMap.cpp
Go to the documentation of this file.
1 /*
2  * visuMap.cpp
3  * this tool will take a source file, target file (h5m) and a map file in nc format, and will visualize weights
4  *
5  * example of usage:
6  * ./mbvisumap -s source.h5m -t target.h5m -m map.nc -b startSourceID \
7  * -e endSourceID -c startTargetID -f endTargetID -o 1
8  * will associate row i, corresponding to target DOF i, in the map, with a partial mesh with entities from source mesh that
9  * target the DOF i; i.e. the weights w(i,j)!=0 , j=1,n_b, will be displayed on source cells with global DOF j
10  * will associate column j corresponding to source DOF j, in the map, with a partial mesh with entities from the target mesh
11  * that are affected by the source DOF j; i.e., the weights w(i,j)!=0, i=1,n_a, will be displayed on target cells with
12  * global DOF i
13  *
14  * The option -o controls if the row and columns files are output in vtk or in h5m format
15  *
16  * can be built only if netcdf and hdf5 and eigen3 are available
17  *
18  * default option is now -o 2, which will create an h5m edge mesh file, with the each edge corresponding to w(i,j)!=0 in the
19  * map file, connecting source center i with target center j. The map will be displayed on a sphere of radius 1,
20  * with the source centers highly elevated from the surfaces, to differentiate them from the target vertices, which
21  * stay on the sphere of radius 1; the elevation is controlled by a new option, -r, with a default value of .05
22  * which means that the source vertices will be put on a sphere of radius 1.05, creating an umbrella for each source center
23  *
24  * only the map file is needed, positions for source and target centers are taken from the map file itself
25  * example of usage:
26  * ./mbvisumap -m map.nc -r 0.01
27  *
28  */
29 #include "moab/MOABConfig.h"
30 
31 #ifndef MOAB_HAVE_EIGEN3
32 #error mbvisumap tool requires eigen3 configuration
33 #endif
34 
35 #ifndef MOAB_HAVE_HDF5
36 #error mbvisumap tool requires hdf5 configuration
37 #endif
38 
39 #include "moab/ProgOptions.hpp"
40 #include "moab/Core.hpp"
41 #include "moab/Range.hpp"
43 #include "moab/ReadUtilIface.hpp"
44 
45 #include "netcdf.h"
46 #include <cmath>
47 #include <sstream>
48 #include <map>
49 #include <Eigen/Sparse>
50 
51 #define ERR_NC( e ) \
52  { \
53  printf( "Error: %s\n", nc_strerror( e ) ); \
54  exit( 2 ); \
55  }
56 
57 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file
58 // ncFile1 is an integer initialized when opening the nc file in read mode
59 
60 int ncFile1;
61 
62 #define GET_DIM1( ncdim, name, val ) \
63  { \
64  int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
65  if( NC_NOERR == gdfail ) \
66  { \
67  size_t tmp_val; \
68  gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
69  if( NC_NOERR != gdfail ) \
70  { \
71  ERR_NC( gdfail ) \
72  } \
73  else \
74  ( val ) = tmp_val; \
75  } \
76  else \
77  ( val ) = 0; \
78  }
79 
80 #define GET_VAR1( name, id, dims ) \
81  { \
82  ( id ) = -1; \
83  int gvfail = nc_inq_varid( ncFile1, name, &( id ) ); \
84  if( NC_NOERR == gvfail ) \
85  { \
86  int ndims; \
87  gvfail = nc_inq_varndims( ncFile1, id, &ndims ); \
88  if( NC_NOERR == gvfail ) \
89  { \
90  ( dims ).resize( ndims ); \
91  gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
92  if( NC_NOERR != gvfail ) \
93  { \
94  ERR_NC( gvfail ) \
95  } \
96  } \
97  } \
98  }
99 
100 #define GET_1D_INT_VAR1( name, id, vals ) \
101  { \
102  GET_VAR1( name, id, vals ); \
103  if( -1 != ( id ) ) \
104  { \
105  size_t ntmp; \
106  int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp ); \
107  if( NC_NOERR != ivfail ) \
108  { \
109  ERR_NC( ivfail ) \
110  } \
111  ( vals ).resize( ntmp ); \
112  size_t ntmp1 = 0; \
113  ivfail = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
114  if( NC_NOERR != ivfail ) \
115  { \
116  ERR_NC( ivfail ) \
117  } \
118  } \
119  }
120 
121 #define GET_1D_DBL_VAR1( name, id, vals ) \
122  { \
123  std::vector< int > dum_dims; \
124  GET_VAR1( name, id, dum_dims ); \
125  if( -1 != ( id ) ) \
126  { \
127  size_t ntmp; \
128  int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp ); \
129  if( NC_NOERR != dvfail ) \
130  { \
131  ERR_NC( dvfail ) \
132  } \
133  ( vals ).resize( ntmp ); \
134  size_t ntmp1 = 0; \
135  dvfail = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
136  if( NC_NOERR != dvfail ) \
137  { \
138  ERR_NC( dvfail ) \
139  } \
140  } \
141  }
142 
143 using namespace moab;
144 using namespace std;
145 
146 int main( int argc, char* argv[] )
147 {
148 
149  ProgOptions opts;
150  int dimSource = 2; // for FV meshes is 2; for SE meshes, use fine mesh, dim will be 0
151  int dimTarget = 2; //
152  int otype = 2;
153  double fraction = 0.05;
154  std::string inputfile1, inputSource, inputTarget;
155  opts.addOpt< std::string >( "map,m", "input map ", &inputfile1 );
156  opts.addOpt< std::string >( "source,s", "source mesh", &inputSource );
157  opts.addOpt< std::string >( "target,t", "target mesh", &inputTarget );
158  opts.addOpt< int >( "dimSource,d", "dimension of source ", &dimSource );
159  opts.addOpt< int >( "dimTarget,g", "dimension of target ", &dimTarget );
160  opts.addOpt< int >( "typeOutput,o", " output type vtk(0), h5m(1), view(default = 2) ", &otype );
161 
162  int startSourceID = -1, endSourceID = -1, startTargetID = -1, endTargetID = -1;
163  opts.addOpt< int >( "startSourceID,b", "start source id ", &startSourceID );
164  opts.addOpt< int >( "endSourceID,e", "end source id ", &endSourceID );
165  opts.addOpt< int >( "startTargetID,c", "start target id ", &startTargetID );
166  opts.addOpt< int >( "endTargetID,f", "end target id ", &endTargetID );
167  opts.addOpt< double >( "raiseFraction,r", "fraction for raising source points height (default 0.05)", &fraction );
168  // -b startSourceID -e endSourceID -c startTargetID -f endTargetID
169 
170  opts.parseCommandLine( argc, argv );
171 
172  std::string extension = ".vtk";
173  if( 1 <= otype ) extension = ".h5m";
174 
175  // Open netcdf map file
176  int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
177  if( NC_NOWRITE != fail )
178  {
179  ERR_NC( fail )
180  }
181 
182  std::cout << " opened " << inputfile1 << " for map 1 \n";
183 
184  int temp_dim;
185  int na1, nb1, ns1;
186  GET_DIM1( temp_dim, "n_a", na1 );
187  GET_DIM1( temp_dim, "n_b", nb1 );
188  GET_DIM1( temp_dim, "n_s", ns1 );
189  std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
190  std::vector< int > col1( ns1 ), row1( ns1 );
191  std::vector< double > val1( ns1 );
192  int idrow1, idcol1, ids1;
193  GET_1D_INT_VAR1( "row", idrow1, row1 );
194  GET_1D_INT_VAR1( "col", idcol1, col1 );
195  GET_1D_DBL_VAR1( "S", ids1, val1 );
196 
197  // we read the matrix; now read moab source and target
198  Core core;
199  Interface* mb = &core;
200  ErrorCode rval;
201  Tag gtag = mb->globalId_tag();
202 
203  // a dense tag for weights
204  Tag wtag;
205  double defVal = 0;
206 
207  std::string name_map = inputfile1;
208  // strip last 3 chars (.nc extension)
209  name_map.erase( name_map.begin() + name_map.length() - 3, name_map.end() );
210  // if path , remove from name
211  size_t pos = name_map.rfind( '/', name_map.length() );
212  if( pos != std::string::npos ) name_map = name_map.erase( 0, pos + 1 );
213 
214  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" );
215 
216  if( 2 == otype )
217  {
218  // create a view of the full map, in which each weight is shown on an edge that starts at the source cell center
219  // and ends at the target cell center
220  // source cell centers are raised a little, let's say a fraction 0.05 * radius, which is 1
221  // each edge gets the associated weight as a tag
222  // first read the cell centers for source and target meshes, directly from the map file
223 
224  std::vector< double > xc_a( na1 ), yc_a( na1 );
225  std::vector< double > xc_b( nb1 ), yc_b( nb1 );
226  int idxc_a, idxc_b, idyc_a, idyc_b;
227  GET_1D_DBL_VAR1( "xc_a", idxc_a, xc_a );
228  GET_1D_DBL_VAR1( "xc_b", idxc_b, xc_b );
229  GET_1D_DBL_VAR1( "yc_a", idyc_a, yc_a );
230  GET_1D_DBL_VAR1( "yc_b", idyc_b, yc_b );
231  // create source vertices, and target vertices
232  std::vector< double > vertex_coords_src( 3 * na1 );
233  // xc_a and yc_a are in degrees, usually
234  for( int i = 0; i < na1; i++ )
235  {
237  sph.R = 1 + fraction; // slightly higher
238  sph.lon = xc_a[i] * M_PI / 180;
239  sph.lat = yc_a[i] * M_PI / 180;
241  vertex_coords_src[3 * i] = pos[0];
242  vertex_coords_src[3 * i + 1] = pos[1];
243  vertex_coords_src[3 * i + 2] = pos[2];
244  }
245  Range source_verts;
246  rval = mb->create_vertices( &vertex_coords_src[0], na1, source_verts );MB_CHK_SET_ERR( rval, "can't create source vertices" );
247  // create a set with source vertices
248  EntityHandle srcSet;
249  rval = mb->create_meshset( MESHSET_SET, srcSet );MB_CHK_SET_ERR( rval, "can't create source set for vertices" );
250  rval = mb->add_entities( srcSet, source_verts );MB_CHK_SET_ERR( rval, "can't add vertices" );
251  std::vector< int > vgid( na1 );
252  for( int i = 0; i < na1; i++ )
253  vgid[i] = i + 1;
254  rval = mb->tag_set_data( gtag, source_verts, &vgid[0] );MB_CHK_SET_ERR( rval, "can't set global id on source verts" );
255 
256  std::vector< double > vertex_coords_tgt( 3 * nb1 );
257  // xc_a and yc_a are in degrees, usually
258  for( int i = 0; i < nb1; i++ )
259  {
261  sph.R = 1; //
262  sph.lon = xc_b[i] * M_PI / 180;
263  sph.lat = yc_b[i] * M_PI / 180;
265  vertex_coords_tgt[3 * i] = pos[0];
266  vertex_coords_tgt[3 * i + 1] = pos[1];
267  vertex_coords_tgt[3 * i + 2] = pos[2];
268  }
269  Range target_verts;
270  rval = mb->create_vertices( &vertex_coords_tgt[0], nb1, target_verts );MB_CHK_SET_ERR( rval, "can't create target vertices" );
271  EntityHandle tgtSet;
272  rval = mb->create_meshset( MESHSET_SET, tgtSet );MB_CHK_SET_ERR( rval, "can't create target set for vertices" );
273  rval = mb->add_entities( tgtSet, target_verts );MB_CHK_SET_ERR( rval, "can't add vertices" );
274  vgid.resize( nb1 );
275  for( int i = 0; i < nb1; i++ )
276  vgid[i] = i + 1;
277  rval = mb->tag_set_data( gtag, target_verts, &vgid[0] );MB_CHK_SET_ERR( rval, "can't set global id on target verts" );
278  // create ns1 edges
279 
280  ReadUtilIface* read_iface;
281  rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
282 
283  EntityHandle actual_start_handle;
284  EntityHandle* array = nullptr;
285  rval = read_iface->get_element_connect( ns1, 2, MBEDGE, 1, actual_start_handle, array );MB_CHK_ERR( rval );
286 
287  for( int i = 0; i < ns1; i++ )
288  {
289  array[2 * i] = source_verts[col1[i] - 1]; // 1 based to 0 based index
290  array[2 * i + 1] = target_verts[row1[i] - 1]; // 1 based to 0 based index
291  }
292  Range edges( actual_start_handle, actual_start_handle + ns1 - 1 );
293 
294  rval = mb->tag_set_data( wtag, edges, &val1[0] );MB_CHK_SET_ERR( rval, "can't set tag on edges" );
295 
296  vgid.resize( ns1 );
297  for( int i = 0; i < ns1; i++ )
298  vgid[i] = i + 1;
299  rval = mb->tag_set_data( gtag, edges, &vgid[0] );MB_CHK_SET_ERR( rval, "can't set global id on edges" );
300 
301  std::string name_file = name_map + extension;
302  rval = mb->write_mesh( name_file.c_str() );MB_CHK_ERR( rval );
303  std::cout << " wrote view map file " << name_file << " with source fraction height: " << fraction << "\n";
304 
305  return 0; // do not bother with other files created, just one file with weights on edges
306  }
307  // first matrix
308  typedef Eigen::Triplet< double > Triplet;
309  std::vector< Triplet > tripletList;
310  tripletList.reserve( ns1 );
311  for( int iv = 0; iv < ns1; iv++ )
312  {
313  // all row and column indices are 1-based in the map file.
314  // inside Eigen3, we will use 0-based; then we will have to add back 1 when
315  // we dump out the files
316  tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
317  }
318  Eigen::SparseMatrix< double > weight1( nb1, na1 );
319 
320  weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
321  weight1.makeCompressed();
322  EntityHandle sourceSet, targetSet;
323  // those are the maps from global ids to the moab entity handles corresponding to those global ids
324  // which are corresponding to the global DOFs
325  map< int, EntityHandle > sourceHandles;
326  map< int, EntityHandle > targetHandles;
327  rval = mb->create_meshset( MESHSET_SET, sourceSet );MB_CHK_SET_ERR( rval, "can't create source mesh set" );
328  rval = mb->create_meshset( MESHSET_SET, targetSet );MB_CHK_SET_ERR( rval, "can't create target mesh set" );
329  const char* readopts = "";
330  rval = mb->load_file( inputSource.c_str(), &sourceSet, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
331  rval = mb->load_file( inputTarget.c_str(), &targetSet, readopts );MB_CHK_SET_ERR( rval, "Failed to read" );
332  Range sRange;
333  rval = mb->get_entities_by_dimension( sourceSet, dimSource, sRange );MB_CHK_SET_ERR( rval, "Failed to get sRange" );
334  vector< int > sids;
335  sids.resize( sRange.size() );
336  rval = mb->tag_get_data( gtag, sRange, &sids[0] );MB_CHK_SET_ERR( rval, "Failed to get ids for srange" );
337  // all global ids are 1 based in the source file, and they correspond to the dofs in the map file
338  for( size_t i = 0; i < sids.size(); i++ )
339  {
340  EntityHandle eh = sRange[i];
341  int gid = sids[i];
342  sourceHandles[gid] = eh;
343  }
344  Range tRange;
345  rval = mb->get_entities_by_dimension( targetSet, dimTarget, tRange );MB_CHK_SET_ERR( rval, "Failed to get tRange" );
346  vector< int > tids;
347  tids.resize( tRange.size() );
348  rval = mb->tag_get_data( gtag, tRange, &tids[0] );MB_CHK_SET_ERR( rval, "Failed to get ids for trange" );
349  // all global ids are 1 based in the target file, and they correspond to the dofs in the map file
350  for( size_t i = 0; i < tids.size(); i++ )
351  {
352  EntityHandle eh = tRange[i];
353  int gid = tids[i];
354  targetHandles[gid] = eh;
355  }
356  EntityHandle partialSet;
357  rval = mb->create_meshset( MESHSET_SET, partialSet );MB_CHK_SET_ERR( rval, "can't create partial set" );
358  // how to get a complete row in sparse matrix? Or a complete column ?
359  for( int col = startSourceID - 1; col <= endSourceID - 1; col++ )
360  {
361  Range targetEnts; // will find entries for column col-1 in sparse matrix weight1, and its entries
362  // will assign a dense tag with values, and write out the file
363  if( col < 0 ) continue;
364  Eigen::SparseVector< double > colVect = weight1.col( col );
365  // the row indices correspond to target cells
366  for( Eigen::SparseVector< double >::InnerIterator it( colVect ); it; ++it )
367  {
368  double weight = it.value(); // == vec[ it.index() ]
369  // we add back the 1 that we subtract
370  int globalIdRow = it.index() + 1;
371  EntityHandle th = targetHandles[globalIdRow];
372  targetEnts.insert( th );
373  rval = mb->tag_set_data( wtag, &th, 1, &weight );MB_CHK_SET_ERR( rval, "Failed to set weight tag on target" );
374  }
375 
376  if( dimTarget == 0 )
377  {
378  Range adjCells;
379  rval = mb->get_adjacencies( targetEnts, 2, false, adjCells, Interface::UNION );MB_CHK_SET_ERR( rval, " can't get adj cells " );
380  targetEnts.merge( adjCells );
381  }
382 
383  rval = mb->add_entities( partialSet, targetEnts );MB_CHK_SET_ERR( rval, "Failed to add target entities to partial set" );
384  // write now the set in a numbered file
385  std::stringstream fff;
386  fff << name_map << "_column" << col + 1 << extension;
387  rval = mb->write_mesh( fff.str().c_str(), &partialSet, 1 );MB_CHK_ERR( rval );
388  // remove from partial set the entities it has
389  rval = mb->clear_meshset( &partialSet, 1 );MB_CHK_SET_ERR( rval, "Failed to empty partial set" );
390  }
391 
392  // how to get a complete row in sparse matrix?
393  for( int row = startTargetID - 1; row <= endTargetID - 1; row++ )
394  {
395  Range sourceEnts; // will find entries for row in sparse matrix weight1, and its entries
396  // will assign a dense tag with values, and write out the file
397  if( row < 0 ) continue;
398  Eigen::SparseVector< double > rowVect = weight1.row( row );
399  // the row indices correspond to target cells
400  for( Eigen::SparseVector< double >::InnerIterator it( rowVect ); it; ++it )
401  {
402  double weight = it.value(); // == vec[ it.index() ]
403  int globalIdCol = it.index() + 1;
404  EntityHandle sh = sourceHandles[globalIdCol];
405  sourceEnts.insert( sh );
406  rval = mb->tag_set_data( wtag, &sh, 1, &weight );MB_CHK_SET_ERR( rval, "Failed to set weight tag on source" );
407  }
408  if( dimSource == 0 )
409  {
410  Range adjCells;
411  rval = mb->get_adjacencies( sourceEnts, 2, false, adjCells, Interface::UNION );MB_CHK_SET_ERR( rval, " can't get adj cells " );
412  sourceEnts.merge( adjCells );
413  }
414  rval = mb->add_entities( partialSet, sourceEnts );MB_CHK_SET_ERR( rval, "Failed to add source entities" );
415  // write now the set in a numbered file
416  std::stringstream fff;
417  fff << name_map << "_row" << row + 1 << extension;
418  rval = mb->write_mesh( fff.str().c_str(), &partialSet, 1 );MB_CHK_ERR( rval );
419  rval = mb->clear_meshset( &partialSet, 1 );MB_CHK_SET_ERR( rval, "Failed to empty partial set" );
420  }
421  return 0;
422 }