 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 }
1.9.1.