1 /*
2 * addncdata.cpp
3 * this tool will take an existing h5m file and add data from some chunk description files
4 * generated from e3sm runs;
5 * will support mainly showing the chunks in ViSit
6 *
7 * example of usage:
8 * ./mbaddchunk -i penta3d.h5m -n chunks_on_proc.txt -o penta3d_ch.h5m
9 *
10 *
11 * Basically, will output a new h5m file (penta3d_ch.h5m), which has 2 extra tags, corresponding to the chunks number
12 * and processors it sits on
13 *
14 *
15 * file penta3d.h5m is obtained from the pentagons file, with a python script
16 *
17 */
18
19
20 #include "moab/ProgOptions.hpp"
21 #include "moab/Core.hpp"
22
23 #include <math.h>
24 #include <sstream>
25 #include <iostream>
26 #include <iomanip>
27 #include <fstream>
28
29 using namespace moab;
30 using namespace std;
31
main(int argc,char * argv[])32 int main(int argc, char* argv[])
33 {
34
35 ProgOptions opts;
36
37 std::string inputfile("penta3d.h5m"), outfile("penta3d_ch.h5m"), chunkfile_name, gsmapfile;
38
39 opts.addOpt<std::string>("input,i", "input mesh filename", &inputfile);
40 opts.addOpt<std::string>("chunkFile,n", "chunk file from cam run", &chunkfile_name);
41 opts.addOpt<std::string>("gsMAPfile,g", "gsmap file", &gsmapfile);
42
43 opts.addOpt<std::string>("output,o", "output mesh filename", &outfile);
44
45
46 opts.parseCommandLine(argc, argv);
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 std::cout << " opened " << inputfile << " with initial h5m data.\n";
54 // open the netcdf file, and see if it has that variable we are looking for
55
56 Range nodes;
57 rval = mb->get_entities_by_dimension(0, 0, nodes);MB_CHK_SET_ERR(rval, "can't get nodes");
58
59 Range edges;
60 rval = mb->get_entities_by_dimension(0, 1, edges);MB_CHK_SET_ERR(rval, "can't get edges");
61
62 Range cells;
63 rval = mb->get_entities_by_dimension(0, 2, cells);MB_CHK_SET_ERR(rval, "can't get cells");
64
65 std::cout << " it has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size() << " cells\n";
66
67 // construct maps between global id and handles
68 std::map<int, EntityHandle> vGidHandle;
69 std::map<int, EntityHandle> eGidHandle;
70 std::map<int, EntityHandle> cGidHandle;
71 std::vector<int> gids;
72 Tag gid;
73 rval = mb->tag_get_handle("GLOBAL_ID", gid);MB_CHK_SET_ERR(rval, "can't get global id tag");
74 gids.resize(nodes.size());
75 rval = mb->tag_get_data(gid, nodes, &gids[0]);MB_CHK_SET_ERR(rval, "can't get global id on vertices");
76 int i=0;
77 for (Range::iterator vit=nodes.begin(); vit!=nodes.end(); vit++)
78 {
79 vGidHandle[gids[i++]] = *vit;
80 }
81
82 gids.resize(edges.size());
83 rval = mb->tag_get_data(gid, edges, &gids[0]);MB_CHK_SET_ERR(rval, "can't get global id on edges");
84 i=0;
85 for (Range::iterator vit=edges.begin(); vit!=edges.end(); vit++)
86 {
87 eGidHandle[gids[i++]] = *vit;
88 }
89
90 gids.resize(cells.size());
91 rval = mb->tag_get_data(gid, cells, &gids[0]);MB_CHK_SET_ERR(rval, "can't get global id on cells");
92 i=0;
93 for (Range::iterator vit=cells.begin(); vit!=cells.end(); vit++)
94 {
95 cGidHandle[gids[i++]] = *vit;
96 }
97
98 if (chunkfile_name.length()>0)
99 {
100
101 // Open chunk file
102 ifstream inFile;
103
104 inFile.open(chunkfile_name.c_str() );
105 if (!inFile) {
106 cout << "Unable to open chunk file";
107 exit(1); // terminate with error
108 }
109 Tag pTag, cTag;
110 int def_val=-1;
111 rval = mb->tag_get_handle("ProcID", 1, MB_TYPE_INTEGER, pTag,
112 MB_TAG_CREAT | MB_TAG_DENSE, &def_val); MB_CHK_SET_ERR(rval, "can't define processor tag");
113 rval = mb->tag_get_handle("ChunkID", 1, MB_TYPE_INTEGER, cTag,
114 MB_TAG_CREAT | MB_TAG_DENSE, &def_val); MB_CHK_SET_ERR(rval, "can't define chunk tag");
115
116 int proc, lcid, ncols ;
117 while (inFile >> proc) {
118 inFile >> lcid >> ncols;
119 int Gid;
120 for (i=0; i<ncols; i++)
121 {
122 inFile >> Gid;
123 EntityHandle cell = cGidHandle[Gid];
124 rval = mb->tag_set_data(pTag, &cell, 1, &proc);MB_CHK_SET_ERR(rval, "can't set proc tag");
125 rval = mb->tag_set_data(cTag, &cell, 1, &lcid);MB_CHK_SET_ERR(rval, "can't set chunk tag");
126
127 }
128 }
129
130 inFile.close();
131 }
132
133 if (gsmapfile.length()>0)
134 {
135
136 // Open chunk file
137 ifstream inFile;
138
139 inFile.open(gsmapfile.c_str() );
140 if (!inFile) {
141 cout << "Unable to open gsmap file";
142 exit(1); // terminate with error
143 }
144 Tag pTag, cTag;
145 int def_val=-1;
146 std::string procTagName = gsmapfile+"_proc";
147 rval = mb->tag_get_handle(procTagName.c_str(), 1, MB_TYPE_INTEGER, pTag,
148 MB_TAG_CREAT | MB_TAG_DENSE, &def_val); MB_CHK_SET_ERR(rval, "can't define processor tag");
149 std::string segTagName = gsmapfile+"_seg";
150 rval = mb->tag_get_handle(segTagName.c_str(), 1, MB_TYPE_INTEGER, cTag,
151 MB_TAG_CREAT | MB_TAG_DENSE, &def_val); MB_CHK_SET_ERR(rval, "can't define segment tag");
152
153 int compid, ngseg, gsize;
154 inFile >> compid >> ngseg >> gsize;
155 for ( i=1; i <= ngseg; i++)
156 {
157 int start, len, pe;
158 inFile >> start >> len >> pe;
159 int Gid;
160 for (int j=0; j<len; j++)
161 {
162 Gid = start + j;
163 EntityHandle cell = cGidHandle[Gid];
164 rval = mb->tag_set_data(pTag, &cell, 1, &pe);MB_CHK_SET_ERR(rval, "can't set proc tag");
165 rval = mb->tag_set_data(cTag, &cell, 1, &i);MB_CHK_SET_ERR(rval, "can't set segment tag");
166 }
167 }
168
169 inFile.close();
170 }
171
172 rval = mb->write_file(outfile.c_str()); MB_CHK_SET_ERR(rval, "can't write file");
173 std::cout << " wrote file " << outfile << "\n";
174 return 0;
175 }
176