Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
VisTags.cpp
Go to the documentation of this file.
1 /** @example VisTags.cpp \n
2  * \brief tool for visualizing multi level tags \n
3  * <b>To run</b>: VisTags <inp_file> <outfile> -O <read_opts> -t <tags> -l <levels> -d <dim> \n
4  *
5  * In this example, it is shown how to create some simple tags for those tags that come from
6  * climate data, multiple levels.
7  * you can read directly nc data, or *.h5m file that will have the tag with multi levels
8  * output will be a vtk file with dense tags of form tag_name_<level>
9  * the tag name might contain a time index too, like T0 or U0
10  * <tag> is a list of tags, separated by commas, no spaces
11  * <levels> is a list of levels, separated by commas, no spaces
12  * dimension of entities with the tags will be specified with -d (default 2)
13  *
14  * an example of use
15  *
16  * VisTags gcrm_r3.nc out.vtk -O VARIABLE=u -t u0,u1 -l 0,1,2 -d 2
17  * (we knew that it had variable u in the file, that it had 256 levels, that there are 2 time
18  * steps, etc)
19  *
20  * or
21  * VisTags gcrm_r3.nc out.vtk -t u0 -l 0,1,2 -d 2
22  * (it will read all variables, but we need to know that u0 will be created as a tag)
23  *
24  * the out.vtk file will contain u0_0, u0_1, as simple dense double tags
25  */
26 
27 #include <iostream>
28 #include <vector>
29 #include <sstream>
30 #include <string>
31 
32 // Include header for MOAB instance and tag conventions
33 #include "moab/Core.hpp"
34 #include "MBTagConventions.hpp"
35 #include "moab/FileOptions.hpp"
36 #ifdef MOAB_HAVE_MPI
37 #include "moab/ParallelComm.hpp"
38 #endif
39 using namespace moab;
40 using namespace std;
41 
42 int main( int argc, char** argv )
43 {
44 #ifdef MOAB_HAVE_NETCDF
45 
46 #ifdef MOAB_HAVE_MPI
47  MPI_Init( &argc, &argv );
48 #endif
49 
50  ErrorCode rval;
51  string file_input, file_output;
52  string read_opts, tags; // Tags to write, separated by commas; it is the name of the tag
53  if( argc < 2 )
54  {
55  file_input = string( MESH_DIR ) + string( "/io/gcrm_r3.nc" );
56  file_output = "VisTagsOut.vtk";
57  }
58  else
59  {
60  file_input = argv[1];
61  file_output = argv[2];
62  }
63  read_opts = "";
64  tags = "";
65 
66  // Instantiate
67  Interface* mb = new( std::nothrow ) Core;
68  if( NULL == mb ) return 1;
69 
70  int dimension = 2;
71  // In MOAB, it may have index after reading (T0, T1, etc)
72  char* levels = NULL; // Levels, separated by commas, no spaces (like 0, 1, 19)
73  if( argc > 3 )
74  {
75  int index = 3;
76  while( index < argc )
77  {
78  if( !strcmp( argv[index], "-O" ) ) // This is for reading options, optional
79  read_opts = argv[++index];
80  if( !strcmp( argv[index], "-t" ) ) tags = argv[++index];
81  if( !strcmp( argv[index], "-l" ) ) levels = argv[++index];
82  if( !strcmp( argv[index], "-d" ) ) dimension = atoi( argv[++index] );
83  index++;
84  }
85  }
86 
87  ostringstream opts;
88  opts << ";;TAGS=" << tags << ";LEVELS=" << levels << "\0";
89  FileOptions fo( opts.str().c_str() );
90 
91  vector< string > tagsNames;
92  vector< int > levelsArray;
93  fo.get_strs_option( "TAGS", tagsNames );
94  fo.get_ints_option( "LEVELS", levelsArray );
95 
96  // Load the input file with the specified options
97  rval = mb->load_file( file_input.c_str(), 0, read_opts.c_str() );MB_CHK_SET_ERR( rval, "not loading file" );
98 
99  Range ents;
100  rval = mb->get_entities_by_dimension( 0, dimension, ents );MB_CHK_SET_ERR( rval, "not getting ents" );
101 
102  // Now create double tags for entities of dimension
103  for( size_t i = 0; i < tagsNames.size(); i++ )
104  {
105  string tagName = tagsNames[i];
106  Tag tagh;
107  rval = mb->tag_get_handle( tagName.c_str(), tagh );
108  if( MB_SUCCESS != rval )
109  {
110  cout << "not getting tag " << tagName.c_str() << "\n";
111  continue;
112  }
113 
114  int len = 0;
115  rval = mb->tag_get_length( tagh, len );
116  if( MB_SUCCESS != rval )
117  {
118  cout << "not getting tag len " << tagName.c_str() << "\n";
119  continue;
120  }
121 
122  DataType type;
123  rval = mb->tag_get_data_type( tagh, type );
124  if( MB_SUCCESS != rval )
125  {
126  cout << "not getting tag type " << tagName.c_str() << "\n";
127  continue;
128  }
129 
130  int count;
131  void* dataptr; // Assume double tags, for simplicity
132  rval = mb->tag_iterate( tagh, ents.begin(), ents.end(), count, dataptr );
133  if( MB_SUCCESS != rval || count != (int)ents.size() )
134  {
135  cout << "not getting tag iterate right " << tagName.c_str() << "\n";
136  continue;
137  }
138 
139  // Now create a new tag, with a new name, concatenated, and copy data there , for each level
140  for( size_t j = 0; j < levelsArray.size(); j++ )
141  {
142  int level = levelsArray[j];
143  if( level >= len )
144  {
145  cout << "level too big at " << level << "\n";
146  continue;
147  }
148 
149  ostringstream newTagName;
150  newTagName << tagName << "_" << level;
151  Tag newTagh;
152  rval = mb->tag_get_handle( newTagName.str().c_str(), 1, type, newTagh, MB_TAG_DENSE | MB_TAG_CREAT );
153  if( MB_SUCCESS != rval )
154  {
155  cout << "not getting new tag " << newTagName.str() << "\n";
156  continue;
157  }
158 
159  void* newDataPtr;
160  rval = mb->tag_iterate( newTagh, ents.begin(), ents.end(), count, newDataPtr );
161  if( MB_SUCCESS != rval || count != (int)ents.size() )
162  {
163  cout << "not getting new tag iterate " << newTagName.str() << "\n";
164  continue;
165  }
166 
167  if( MB_TYPE_DOUBLE == type )
168  {
169  double* ptrD = (double*)newDataPtr;
170  double* oldData = (double*)dataptr;
171  for( int k = 0; k < count; k++, ptrD++ )
172  *ptrD = oldData[level + count * k];
173  }
174  } // for (size_t j = 0; j < levelsArray.size(); j++)
175 
176  mb->tag_delete( tagh ); // No need for the tag anymore, write it to the new file
177  } // for (size_t i = 0; i < tagsNames.size(); i++)
178 
179  rval = mb->write_file( file_output.c_str() );MB_CHK_SET_ERR( rval, "Can't write file " << file_output );
180  cout << "Successfully wrote file " << file_output << "\n";
181 
182  delete mb;
183 
184 #ifdef MOAB_HAVE_MPI
185  MPI_Finalize();
186 #endif
187 #else
188  std::cout << " configure with netcdf for this example to work\n";
189 #endif
190  return 0;
191 }