Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
t3d2mesh.C
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4 #include <string>
5 #include <list>
6 #include <vector>
7 #include <algorithm>
8 
9 #include "Mesh.H"
10 #include "Profiler.H"
11 #include "Global.H"
12 #include "primitive_utilities.H"
13 
14 class T3DGeomEnt : public std::pair<unsigned int,unsigned int>
15 {
16 public:
17  T3DGeomEnt() : std::pair<unsigned int,unsigned int>()
18  {};
19  T3DGeomEnt(const T3DGeomEnt &ing)
20  {
21  this->first = ing.first;
22  this->second = ing.second;
23  };
24  const T3DGeomEnt &operator=(const T3DGeomEnt &ing)
25  {
26  this->first = ing.first;
27  this->second = ing.second;
28  return(*this);
29  };
30  bool operator<(const T3DGeomEnt &ge) const { return (first < ge.first ? true :
31  first > ge.first ? false : second < ge.second); };
32  bool operator==(const T3DGeomEnt &ge) const { return (first == ge.first ? second==ge.second : false); };
33 };
34 
35 
36 int
37 main(int argc,char *argv[])
38 {
39  std::string line;
40  std::istringstream Istr;
41  unsigned int mesh_type = 0;
42  unsigned int elem_order = 0;
43  unsigned int nnodes = 0;
44  unsigned int nedges = 0;
45  unsigned int ntris = 0;
46  unsigned int nquads = 0;
47  unsigned int ntets = 0;
48  unsigned int nhexs = 0;
49  unsigned int npyrs = 0;
50  unsigned int npriss = 0;
51 
52 
53  // Read Block 1
54  std::getline(std::cin,line);
55  Istr.str(line);
56  Istr >> mesh_type >> elem_order;
57  std::getline(std::cin,line);
58  Istr.str(line);
59  Istr >> nnodes >> nedges;
60  switch(mesh_type){
61  case 3:
62  Istr >> ntris >> ntets;
63  break;
64  case 4:
65  Istr >> nquads >> nhexs;
66  break;
67  case 7:
68  Istr >> ntris >> nquads >> ntets >> npyrs >> npriss >> nhexs;
69  break;
70  default:
71  std::cerr << "Fatal Error: Unsupported mesh type, "
72  << mesh_type << "." << std::endl;
73  return(1);
74  break;
75 
76  }
77  std::getline(std::cin,line);
78 
79  // Calculate some mesh parameters
80  unsigned int nelems;
81  unsigned int n_2d_elem = ntris + nquads;
82  unsigned int n_3d_elem = ntets + npyrs + npriss + nhexs;
83  if(n_2d_elem > 0 && n_3d_elem > 0){
84  std::cerr << "Fatal Error: Mixing 2d and 3d elements not supported at this time."
85  << std::endl;
86  return(1);
87  }
88  nelems = n_2d_elem > 0 ? n_2d_elem : n_3d_elem;
89 
90  std::cout << nnodes << std::endl;
91  // Read Block 2 (nodes)
92  // Open a file to store properties
93  std::ofstream PropOut;
94  std::list<T3DGeomEnt> geomlist;
95  std::list<int> nodal_properties;
96  std::list<int> elemental_properties;
97  std::vector<T3DGeomEnt> nodal_geoms(nnodes);
98  std::vector<T3DGeomEnt> elem_geoms(nelems);
99  for(unsigned int i = 0;i < nnodes;i++){
100  T3DGeomEnt ge;
101  std::getline(std::cin,line);
102  int superfluous_id;
103  Istr.clear();
104  Istr.str(line);
105  Istr >> superfluous_id;
106  double x,y,z;
107  Istr >> x >> y >> z;
108  std::cout << x << "\t" << y << "\t" << z << std::endl;
109  int property = 0;
110  Istr >> nodal_geoms[i].first >> nodal_geoms[i].second >> property;
111  // std::getline(Istr,line);
112  // std::string::size_type xpos = line.find_first_not_of(" ");
113  // line = line.substr(xpos);
114  // PropOut << line << std::endl;
115  nodal_properties.push_back(property);
116  }
117  // add a blank line between nodal and elemental properties
118  // PropOut << std::endl;
119  std::getline(std::cin,line);
120  std::cout << nelems << std::endl;
121 
122  // Skip Block 3 (edges)
123  if(nedges > 0){
124  std::cerr << "Warning: skipping edges." << std::endl;
125  for(unsigned int i = 0;i < nedges;i++)
126  std::getline(std::cin,line);
127  std::getline(std::cin,line);
128  }
129 
130  // Read Block 4 (triangle elements)
131  if(mesh_type != 4 && ntris > 0){
132  for(unsigned int i = 0;i < ntris;i++){
133  std::getline(std::cin,line);
134  std::string::size_type xpos = line.find_first_not_of(" ");
135  line = line.substr(xpos);
136  Istr.clear();
137  Istr.str(line);
138  unsigned int superfluous_id;
139  Istr >> superfluous_id;
140  unsigned int node_id;
141  for(unsigned int j = 0;j < elem_order*3;j++){
142  Istr >> node_id;
143  std::cout << node_id;
144  if(j != (elem_order*3 - 1)) std::cout << "\t";
145  }
146  std::cout << std::endl;
147  int property = 0;
148  Istr >> elem_geoms[i].first >> elem_geoms[i].second >> property;
149  elemental_properties.push_back(property);
150  // std::getline(Istr,line);
151  // std::string::size_type x = line.find_first_not_of(" ");
152  // line = line.substr(x);
153  // PropOut << line << std::endl;
154  }
155  std::getline(std::cin,line);
156  }
157  Mesh::IndexType egoffset = ntris;
158 
159  // Read Block 5 (quad elements)
160  if(mesh_type != 3 && nquads > 0){
161  for(unsigned int i = 0;i < nquads;i++){
162  std::getline(std::cin,line);
163  std::string::size_type xpos = line.find_first_not_of(" ");
164  line = line.substr(xpos);
165  Istr.clear();
166  Istr.str(line);
167  unsigned int superfluous_id;
168  Istr >> superfluous_id;
169  unsigned int node_id;
170  for(unsigned int j = 0;j < elem_order*4;j++){
171  Istr >> node_id;
172  std::cout << node_id;
173  if(j != (elem_order*4 - 1)) std::cout << "\t";
174  }
175  std::cout << std::endl;
176  int property = 0;
177  Istr >> elem_geoms[i+egoffset].first >> elem_geoms[i+egoffset].second >> property;
178  elemental_properties.push_back(property);
179  }
180  std::getline(std::cin,line);
181  }
182  egoffset += nquads;
183  // Read Block 6 (tet elements)
184  if(mesh_type != 4 && ntets > 0){
185  unsigned int elem_size = elem_order==1 ? 4 : 10;
186  for(unsigned int i = 0;i < ntets;i++){
187  std::getline(std::cin,line);
188  std::string::size_type xpos = line.find_first_not_of(" ");
189  line = line.substr(xpos);
190  Istr.clear();
191  Istr.str(line);
192  unsigned int superfluous_id;
193  Istr >> superfluous_id;
194  unsigned int node_id;
195  for(unsigned int j = 0;j < elem_size;j++){
196  Istr >> node_id;
197  std::cout << node_id;
198  if(j != (elem_size - 1)) std::cout << "\t";
199  }
200  std::cout << std::endl;
201  int property = 0;
202  Istr >> elem_geoms[i+egoffset].first >> elem_geoms[i+egoffset].second >> property;
203  elemental_properties.push_back(property);
204  }
205  std::getline(std::cin,line);
206  }
207  egoffset += ntets;
208  // Read Block 7 (pyr elements)
209  if(mesh_type == 7 && npyrs > 0){
210  unsigned int elem_size = elem_order==1 ? 5 : 8;
211  for(unsigned int i = 0;i < npyrs;i++){
212  std::getline(std::cin,line);
213  std::string::size_type xpos = line.find_first_not_of(" ");
214  line = line.substr(xpos);
215  Istr.clear();
216  Istr.str(line);
217  unsigned int superfluous_id;
218  Istr >> superfluous_id;
219  unsigned int node_id;
220  for(unsigned int j = 0;j < elem_size;j++){
221  Istr >> node_id;
222  std::cout << node_id;
223  if(j != (elem_size - 1)) std::cout << "\t";
224  }
225  std::cout << std::endl;
226  int property = 0;
227  Istr >> elem_geoms[i+egoffset].first >> elem_geoms[i+egoffset].second >> property;
228  elemental_properties.push_back(property);
229  }
230  std::getline(std::cin,line);
231  }
232 
233  egoffset += npyrs;
234 
235  // Read Block 8 (pris elements)
236  if(mesh_type == 7 && npriss > 0){
237  unsigned int elem_size = elem_order==1 ? 6 : 15;
238  for(unsigned int i = 0;i < npriss;i++){
239  std::getline(std::cin,line);
240  std::string::size_type xpos = line.find_first_not_of(" ");
241  line = line.substr(xpos);
242  Istr.clear();
243  Istr.str(line);
244  unsigned int superfluous_id;
245  Istr >> superfluous_id;
246  unsigned int node_id;
247  for(unsigned int j = 0;j < elem_size;j++){
248  Istr >> node_id;
249  std::cout << node_id;
250  if(j != (elem_size - 1)) std::cout << "\t";
251  }
252  std::cout << std::endl;
253  int property = 0;
254  Istr >> elem_geoms[i+egoffset].first >> elem_geoms[i+egoffset].second >> property;
255  elemental_properties.push_back(property);
256  }
257  std::getline(std::cin,line);
258  }
259  egoffset += npriss;
260  // Read Block 9 (hex elements)
261  if(mesh_type != 3 && nhexs > 0){
262  unsigned int elem_size = elem_order==1 ? 8 : 20;
263  for(unsigned int i = 0;i < nhexs;i++){
264  std::getline(std::cin,line);
265  std::string::size_type xpos = line.find_first_not_of(" ");
266  line = line.substr(xpos);
267  Istr.clear();
268  Istr.str(line);
269  unsigned int superfluous_id;
270  Istr >> superfluous_id;
271  unsigned int node_id;
272  for(unsigned int j = 0;j < elem_size;j++){
273  Istr >> node_id;
274  std::cout << node_id;
275  if(j != (elem_size - 1)) std::cout << "\t";
276  }
277  std::cout << std::endl;
278  int property = 0;
279  Istr >> elem_geoms[i+egoffset].first >> elem_geoms[i+egoffset].second >> property;
280  elemental_properties.push_back(property);
281  }
282  std::getline(std::cin,line);
283  }
284  // PropOut.close();
285  geomlist.resize(0);
286  // Now sort out the geometry information - real slow
287  std::vector<T3DGeomEnt>::iterator gi = nodal_geoms.begin();
288  while(gi != nodal_geoms.end()){
289  std::list<T3DGeomEnt>::iterator gli = std::find(geomlist.begin(),geomlist.end(),*gi);
290  if(gli == geomlist.end()){
291  geomlist.push_back(*gi);
292  }
293  gi++;
294  }
295  gi = elem_geoms.begin();
296  while(gi != elem_geoms.end())
297  {
298  std::list<T3DGeomEnt>::iterator gli = std::find(geomlist.begin(),geomlist.end(),*gi);
299  if(gli == geomlist.end())
300  geomlist.push_back(*gi);
301  gi++;
302  }
303  geomlist.sort();
304  std::list<T3DGeomEnt>::iterator gli = geomlist.begin();
305  std::cerr << "GEOMETRIES:" << std::endl;
306  while(gli != geomlist.end())
307  {
308  std::cerr << "Dimension: " << gli->first << ", Id = " << gli->second << std::endl;
309  gli++;
310  }
311  Mesh::IndexType ngeoms = geomlist.size();
312  Mesh::IndexType nvertex = 0;
313  Mesh::IndexType ncurve = 0;
314  Mesh::IndexType nsurf = 0;
315  Mesh::IndexType nregion = 0;
316  Mesh::IndexType npatch = 0;
317  // Mesh::IndexType nshell = 0;
318  // Mesh::IndexType ninter = 0;
319  std::vector<Mesh::GeometricEntity> geometries(ngeoms);
320  std::vector<T3DGeomEnt> t3d_geovec(ngeoms);
321  gli = geomlist.begin();
322  while(gli != geomlist.end()){
323  ubyte dim = (gli->first == 1 ? 0 :
324  gli->first == 2 ? 1 :
325  gli->first == 3 ? 2 :
326  gli->first == 4 ? 3 :
327  gli->first == 5 ? 2 :
328  gli->first == 6 ? 2 :
329  gli->first);
330  if(dim > 3){
331  std::cerr << "Error: Cannot deal with interface mesh type"
332  << std::endl;
333  return(1);
334  }
335  std::ostringstream Ostr;
336  if(dim == 0)
337  Ostr << "Vertex-" << ++nvertex;
338  if(dim == 1)
339  Ostr << "Curve-" << ++ncurve;
340  if(dim == 2)
341  Ostr << "Surface-" << ++nsurf;
342  if(dim == 3)
343  Ostr << "Region-" << ++nregion;
344  Mesh::IndexType index = nvertex+ncurve+nsurf+nregion-1;
345  t3d_geovec[index] = *gli;
346  geometries[index].first = Ostr.str();
347  geometries[index].second = dim;
348  geometries[index]._collections.resize((unsigned int)dim+1);
349  gli++;
350  }
351  geomlist.resize(0);
352  std::cerr << "Geometries created, populating..." << std::endl;
353  gi = nodal_geoms.begin();
354  Mesh::IndexType nid = 1;
355  std::cerr << " populating nodal geoms..." << std::endl;
356  PropOut.open("t3d_properties");
357  std::list<int>::iterator pi = nodal_properties.begin();
358  while(gi != nodal_geoms.end()){
359  Mesh::IndexType gindex = std::find(t3d_geovec.begin(),t3d_geovec.end(),*gi++) -
360  t3d_geovec.begin();
361  geometries[gindex]._collections[0].push_back(nid++);
362  PropOut << gindex+1 << "\t" << *pi++ << std::endl;
363  }
364  // blank line between nodal and elemental properties
365  PropOut << std::endl;
366  gi = elem_geoms.begin();
367  pi = elemental_properties.begin();
368  std::cerr << " populating elemental geoms..." << std::endl;
369  Mesh::IndexType eid = 1;
370  while(gi != elem_geoms.end()){
371  Mesh::IndexType gindex = std::find(t3d_geovec.begin(),t3d_geovec.end(),*gi++) -
372  t3d_geovec.begin();
373  geometries[gindex]._collections[geometries[gindex].second].push_back(eid++);
374  PropOut << gindex+1 << "\t" << *pi++ << std::endl;
375  }
376  PropOut.close();
377  std::ofstream GeomOut;
378  GeomOut.open("t3d_geometries");
379  GeomOut << geometries.size() << std::endl;
380  for(Mesh::IndexType gindex = 0;gindex < geometries.size();gindex++)
381  GeomOut << geometries[gindex].first << "\t"
382  << (uint)geometries[gindex].second << std::endl;
383  for(Mesh::IndexType gindex = 0;gindex < geometries.size();gindex++)
384  GeomOut << geometries[gindex] << std::endl;
385  GeomOut.close();
386  return(0);
387 }
bool operator<(const T3DGeomEnt &ge) const
Definition: t3d2mesh.C:30
T3DGeomEnt()
Definition: t3d2mesh.C:17
void int int REAL REAL * y
Definition: read.cpp:74
Mesh Stuff.
void int int int REAL REAL REAL * z
Definition: write.cpp:76
const T3DGeomEnt & operator=(const T3DGeomEnt &ing)
Definition: t3d2mesh.C:24
blockLoc i
Definition: read.cpp:79
void int int REAL * x
Definition: read.cpp:74
int main(int argc, char *argv[])
Definition: blastest.C:94
j indices j
Definition: Indexing.h:6
const double pi
bool operator==(const T3DGeomEnt &ge) const
Definition: t3d2mesh.C:32
IRAD::Primitive::IndexType IndexType
Definition: Mesh.H:57
T3DGeomEnt(const T3DGeomEnt &ing)
Definition: t3d2mesh.C:19