1 //
2 //  Read user dependent output file
3 //
4 
5 
6 #include <mystdlib.h>
7 #include <myadt.hpp>
8 #include <linalg.hpp>
9 #include <csg.hpp>
10 #include <stlgeom.hpp>
11 #include <meshing.hpp>
12 
13 #include "writeuser.hpp"
14 
15 namespace netgen
16 {
ReadFile(Mesh & mesh,const string & hfilename)17   void ReadFile (Mesh & mesh,
18                  const string & hfilename)
19   {
20     PrintMessage(3, "Read User File");
21 
22     const char * filename = hfilename.c_str();
23 
24     char reco[100];
25     int np, nbe;
26 
27 
28 
29     // ".surf" - mesh
30 
31     if ( (strlen (filename) > 5) &&
32          strcmp (&filename[strlen (filename)-5], ".surf") == 0 )
33 
34       {
35         cout << "Surface file" << endl;
36 
37         ifstream in (filename);
38 
39         in >> reco;
40         in >> np;
41         for (int i = 1; i <= np; i++)
42           {
43             Point3d p;
44             in >> p.X() >> p.Y() >> p.Z();
45             mesh.AddPoint (p);
46           }
47 
48         mesh.ClearFaceDescriptors();
49         mesh.AddFaceDescriptor (FaceDescriptor(1,1,0,0));
50 
51         in >> nbe;
52         //      int invert = globflags.GetDefineFlag ("invertsurfacemesh");
53         for (int i = 1; i <= nbe; i++)
54           {
55             Element2d el;
56             el.SetIndex(1);
57 
58             for (int j = 1; j <= 3; j++)
59               {
60                 in >> el.PNum(j);
61                 // el.PNum(j)++;
62                 if (el.PNum(j) < PointIndex(1) ||
63                     el.PNum(j) > PointIndex(np))
64                   {
65                     cerr << "Point Number " << el.PNum(j) << " out of range 1..."
66                          << np << endl;
67                     return;
68                   }
69               }
70             /*
71               if (invert)
72               swap (el.PNum(2), el.PNum(3));
73             */
74 
75             mesh.AddSurfaceElement (el);
76           }
77 
78 
79         cout << "points: " << np << " faces: " << nbe << endl;
80       }
81 
82 
83 
84 
85 
86     if ( (strlen (filename) > 4) &&
87          strcmp (&filename[strlen (filename)-4], ".unv") == 0 )
88       {
89         char reco[100];
90         // int invert;
91 
92         ifstream in(filename);
93 
94         mesh.ClearFaceDescriptors();
95         mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0));
96         mesh.GetFaceDescriptor(1).SetBCProperty (1);
97         // map from unv element nr to our element number + an index if it is vol (0), bnd(1), ...
98         std::map<size_t, std::tuple<size_t, int>> element_map;
99 	int dim = 3;
100 	int bccounter = 0;
101 
102         NgArray<Segment> tmp_segments;
103         while (in.good())
104           {
105             in >> reco;
106             if (strcmp(reco, "-1") == 0)
107               continue;
108 
109             else if (strcmp (reco, "2411") == 0)
110               {
111                 cout << "nodes found" << endl;
112 
113                 while (1)
114                   {
115                     int pi, hi;
116                     Point<3> p;
117 
118                     in >> pi;
119                     if (pi == -1)
120                       break;
121 
122                     in >> hi >> hi >> hi;
123                     in >> p(0) >> p(1) >> p(2);
124 
125                     mesh.AddPoint (p);
126                   }
127 		cout << "read " << mesh.GetNP() << " points" << endl;
128                 Point3d pmin, pmax;
129                 mesh.GetBox (pmin, pmax);
130                 if(fabs(pmin.Z() - pmax.Z()) < 1e-10 * Dist(pmin, pmax))
131                 {
132                        cout << "Set Dimension to 2." << endl;
133                        mesh.SetDimension(2);
134                        dim = 2 ;
135 		}
136 
137               }
138 
139             else if (strcmp (reco, "2412") == 0)
140               {
141                 cout << "elements found" << endl;
142 
143                 while (1)
144                   {
145 		    int label, fe_id, phys_prop, mat_prop, color, nnodes;
146 		    int nodes[100];
147 		    int hi;
148 
149 		    in >> label;
150 		    if (label == -1) break;
151 		    in >> fe_id >> phys_prop >> mat_prop >> color >> nnodes;
152 
153 		    if (fe_id >= 11 && fe_id <= 32)
154 		      in >> hi >> hi >> hi;
155 
156 
157 		    for (int j = 0; j < nnodes; j++)
158 		      in >> nodes[j];
159 
160 		    switch (fe_id)
161 		      {
162                       case 11: // (Rod) SEGM
163                         {
164                           Segment el;
165                           el[0] = nodes[0];
166                           el[1] = nodes[1];
167                           el[2] = -1;
168 
169                           if(dim == 3){
170                             auto nr = tmp_segments.Size();
171                             tmp_segments.Append(el);
172                             element_map[label] = std::make_tuple(nr+1, 2);
173                           }
174                           else if(dim == 2){
175 		            el.si = -1; // add label to segment, will be changed later when BC's are assigned
176                             auto nr = mesh.AddSegment(el);
177                             element_map[label] = std::make_tuple(nr+1, 2);
178                           }
179                           break;
180                         }
181 
182                       case 22: // (Tapered beam) SEGM
183                         {
184                           Segment el;
185                           el[0] = nodes[0];
186                           el[1] = nodes[2];
187                           el[2] = nodes[1];
188 
189                           if(dim == 3){
190                             auto nr = tmp_segments.Size();
191                             tmp_segments.Append(el);
192                             element_map[label] = std::make_tuple(nr+1, 2);
193                           }
194                           else if(dim == 2){
195 		            el.si = -1; // add label to segment, will be changed later when BC's are assigned
196                             auto nr = mesh.AddSegment(el);
197                             element_map[label] = std::make_tuple(nr+1, 2);
198                           }
199 
200                           break;
201                         }
202 		      case 41: // TRIG
203 			{
204 			  Element2d el (TRIG);
205 			  el.SetIndex (1);
206 			  for (int j = 0; j < nnodes; j++)
207 			    el[j] = nodes[j];
208 			  auto nr = mesh.AddSurfaceElement (el);
209                           element_map[label] = std::make_tuple(nr+1, 1);
210 			  break;
211 			}
212                       case 42: // TRIG6
213                         {
214                           Element2d el(TRIG6);
215                           el.SetIndex(1);
216                           int jj = 0;
217                           for(auto j : {0,2,4,3,5,1})
218                               el[jj++] = nodes[j];
219                           auto nr = mesh.AddSurfaceElement(el);
220                           element_map[label] = std::make_tuple(nr+1, 1);
221                           break;
222                         }
223 		      case 111: // TET
224 			{
225 			  Element el (TET);
226 			  el.SetIndex (1);
227 			  for (int j = 0; j < nnodes; j++)
228 			    el[j] = nodes[j];
229 			  auto nr = mesh.AddVolumeElement (el);
230 			  element_map[label] = std::make_tuple(nr+1, 0);
231 			  break;
232 			}
233                       case 118: // TET10
234                         {
235                           Element el(TET10);
236                           el.SetIndex(1);
237                           int jj = 0;
238                           for(auto j : {0,2,4,9,1,5,6,3,7,8})
239                             el[jj++] = nodes[j];
240                           auto nr = mesh.AddVolumeElement(el);
241                           element_map[label] = std::make_tuple(nr+1, 0);
242                           break;
243                         }
244                       default:
245                         cout << "Do not know fe_id = " << fe_id << ", skipping it." << endl;
246                         break;
247 		      }
248                   }
249                 cout << mesh.GetNE() << " elements found" << endl;
250                 cout << mesh.GetNSE() << " surface elements found" << endl;
251 
252               }
253             else if(strcmp (reco, "2467") == 0)
254               {
255                 int matnr = 1;
256                 cout << "Groups found" << endl;
257                 while(in.good())
258                   {
259                     int len;
260                     string name;
261                     in >> len;
262                     if(len == -1)
263                       break;
264                     for(int i=0; i < 7; i++)
265                       in >> len;
266                     in >> name;
267                     cout << len << " element are in group " << name << endl;
268                     int hi, index;
269                     int fdnr, ednr;
270 
271                     in >> hi >> index >> hi >> hi;
272                     int codim = get<1>(element_map[index]);
273                     // use first element to determine if boundary or volume
274 
275                     switch (codim)
276                       {
277                       case 0:
278                         {
279                           mesh.SetMaterial(++matnr, name);
280                           mesh.VolumeElement(get<0>(element_map[index])).SetIndex(matnr);
281                           break;
282                         }
283                       case 1:
284                         {
285                           if(dim == 3)
286                           {
287                             int bcpr = mesh.GetNFD();
288                             fdnr = mesh.AddFaceDescriptor(FaceDescriptor(bcpr, 0,0,0));
289                             mesh.GetFaceDescriptor(fdnr).SetBCProperty(bcpr+1);
290                             mesh.SetBCName(bcpr, name);
291                             mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(fdnr);
292                             bccounter++;
293                           }
294                           else if(dim == 2)
295                           {
296                             mesh.SetMaterial(matnr, name);
297                             fdnr = mesh.AddFaceDescriptor(FaceDescriptor(matnr, 0,0,0));
298                             mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(matnr);
299                             mesh.GetFaceDescriptor(fdnr).SetBCProperty(matnr);
300 			    matnr++;
301                           }
302                           break;
303 
304                         }
305                       case 2:
306                         {
307                          if(dim == 3)
308                           {
309                             int bcpr = mesh.GetNCD2Names()+1;
310                             auto ed = EdgeDescriptor();
311                             ed.SetSurfNr(0,bcpr);//?
312                             ednr = mesh.AddEdgeDescriptor(ed);
313                             mesh.SetCD2Name(bcpr, name);
314                             auto nr = mesh.AddSegment(tmp_segments[get<0>(element_map[index])-1]);
315                             mesh[nr].edgenr = ednr+1;
316                           }
317                           else if(dim == 2)
318                           {
319                             Segment & seg = mesh.LineSegment(get<0>(element_map[index]));
320 			    seg.si = bccounter + 1;
321 			    mesh.SetBCName(bccounter, name);
322 		            bccounter++;
323                           }
324                           break;
325 
326                         }
327                       default:
328                         {
329                           cout << "Codim " << codim << " not implemented yet!" << endl;
330                         }
331                       }
332 
333                     for(int i=0; i<len-1; i++)
334                       {
335                         in >> hi >> index >> hi >> hi;
336                         switch (codim)
337                           {
338                           case 0:
339                             mesh.VolumeElement(get<0>(element_map[index])).SetIndex(matnr);
340                             break;
341                           case 1:
342 			    if(dim == 3) mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(fdnr);
343 			    else if (dim == 2){
344                                     mesh.SurfaceElement(get<0>(element_map[index])).SetIndex(matnr-1);
345 				    mesh.GetFaceDescriptor(fdnr).SetBCProperty(matnr);
346 			    }
347                             break;
348                           case 2:
349 	   		    if(dim == 3)
350                             {
351                               auto nr = mesh.AddSegment(tmp_segments[get<0>(element_map[index])-1]);
352                               mesh[nr].edgenr = ednr+1;
353                             }
354 			    else if(dim == 2)
355 			    {
356 	 			    Segment & seg = mesh.LineSegment(get<0>(element_map[index]));
357 			            seg.si = bccounter;
358 			    }
359                             break;
360                           default:
361                             break;
362                           }
363                       }
364                   }
365               }
366             else
367               {
368                 cout << "Do not know data field type " << reco << ", skipping it" << endl;
369                 while(in.good())
370                   {
371                     in >> reco;
372                     if(strcmp(reco, "-1") == 0)
373                       break;
374                   }
375               }
376           }
377 
378 	if(dim == 2){
379 		// loop through segments to assign default BC to unmarked edges
380 		int bccounter_tmp = bccounter;
381 		for(int index=1; index <= mesh.GetNSeg(); index++){
382                 	Segment & seg = mesh.LineSegment(index);
383 			if(seg.si == -1){
384 			  seg.si = bccounter + 1;
385 			  if(bccounter_tmp == bccounter) mesh.SetBCName(bccounter, "default"); // could be more efficient
386 			  bccounter_tmp++;
387 			}
388 		}
389 		if(bccounter_tmp > bccounter) bccounter++;
390 	}
391 
392 
393         Point3d pmin, pmax;
394         mesh.ComputeNVertices();
395         mesh.RebuildSurfaceElementLists();
396         mesh.GetBox (pmin, pmax);
397         mesh.UpdateTopology();
398         if(dim == 3) bccounter++;
399         cout << "bounding-box = " << pmin << "-" << pmax << endl;
400 	cout << "Created " << bccounter << " boundaries." << endl;
401 	for(int i=0; i<bccounter; i++){
402 		cout << mesh.GetBCName(i) << endl;
403 	}
404       }
405 
406 
407 
408     // fepp format2d:
409 
410     if ( (strlen (filename) > 7) &&
411          strcmp (&filename[strlen (filename)-7], ".mesh2d") == 0 )
412       {
413         cout << "Reading FEPP2D Mesh" << endl;
414 
415         char buf[100];
416         int np, ne, nseg, i, j;
417 
418         ifstream in (filename);
419 
420         in >> buf;
421 
422         in >> nseg;
423         for (i = 1; i <= nseg; i++)
424           {
425             int bound, p1, p2;
426             in >> bound >> p1 >> p2;
427             // forget them
428           }
429 
430         in >> ne;
431         for (i = 1; i <= ne; i++)
432           {
433             int mat, nelp;
434             in >> mat >> nelp;
435             Element2d el (nelp == 3 ? TRIG : QUAD);
436             el.SetIndex (mat);
437             for (j = 1; j <= nelp; j++)
438               in >> el.PNum(j);
439             mesh.AddSurfaceElement (el);
440           }
441 
442         in >> np;
443         for (i = 1; i <= np; i++)
444           {
445             Point3d p(0,0,0);
446             in >> p.X() >> p.Y();
447             mesh.AddPoint (p);
448           }
449       }
450 
451 
452     else if ( (strlen (filename) > 5) &&
453               strcmp (&filename[strlen (filename)-5], ".mesh") == 0 )
454       {
455         cout << "Reading Neutral Format" << endl;
456 
457         int np, ne, nse, i, j;
458 
459         ifstream in (filename);
460 
461         in >> np;
462 
463         if (in.good())
464           {
465             // file starts with an integer
466 
467             for (i = 1; i <= np; i++)
468               {
469                 Point3d p(0,0,0);
470                 in >> p.X() >> p.Y() >> p.Z();
471                 mesh.AddPoint (p);
472               }
473 
474             in >> ne;
475             for (i = 1; i <= ne; i++)
476               {
477                 int mat;
478                 in >> mat;
479                 Element el (4);
480                 el.SetIndex (mat);
481                 for (j = 1; j <= 4; j++)
482                   in >> el.PNum(j);
483                 mesh.AddVolumeElement (el);
484               }
485 
486             mesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 0));
487             int nfd = 1;
488 
489             in >> nse;
490             for (i = 1; i <= nse; i++)
491               {
492                 int mat; // , nelp;
493                 in >> mat;
494                 Element2d el (TRIG);
495                 el.SetIndex (mat);
496                 while(nfd<mat)
497                   {
498                     ++nfd;
499                     mesh.AddFaceDescriptor(FaceDescriptor(nfd,nfd,0,0));
500                   }
501                 for (j = 1; j <= 3; j++)
502                   in >> el.PNum(j);
503                 mesh.AddSurfaceElement (el);
504               }
505           }
506         else
507           {
508             char buf[100];
509             in.clear();
510             do
511               {
512                 in >> buf;
513                 cout << "buf = " << buf << endl;
514                 if (strcmp (buf, "points") == 0)
515                   {
516                     in >> np;
517                     cout << "np = " << np << endl;
518                   }
519               }
520             while (in.good());
521           }
522       }
523 
524 
525     if ( (strlen (filename) > 4) &&
526          strcmp (&filename[strlen (filename)-4], ".emt") == 0 )
527       {
528         ifstream inemt (filename);
529 
530         string pktfile = filename;
531         int len = strlen (filename);
532         pktfile[len-3] = 'p';
533         pktfile[len-2] = 'k';
534         pktfile[len-1] = 't';
535         cout << "pktfile = " << pktfile << endl;
536 
537         int np, nse, i;
538         int bcprop;
539         ifstream inpkt (pktfile.c_str());
540         inpkt >> np;
541         NgArray<double> values(np);
542         for (i = 1; i <= np; i++)
543           {
544             Point3d p(0,0,0);
545             inpkt >> p.X() >> p.Y() >> p.Z()
546                   >> bcprop >> values.Elem(i);
547             mesh.AddPoint (p);
548           }
549 
550         mesh.ClearFaceDescriptors();
551         mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0));
552         mesh.GetFaceDescriptor(1).SetBCProperty (1);
553         mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0));
554         mesh.GetFaceDescriptor(2).SetBCProperty (2);
555         mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0));
556         mesh.GetFaceDescriptor(3).SetBCProperty (3);
557         mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0));
558         mesh.GetFaceDescriptor(4).SetBCProperty (4);
559         mesh.AddFaceDescriptor (FaceDescriptor(0,1,0,0));
560         mesh.GetFaceDescriptor(5).SetBCProperty (5);
561 
562         int p1, p2, p3;
563         double value;
564         inemt >> nse;
565         for (i = 1; i <= nse; i++)
566           {
567             inemt >> p1 >> p2 >> p3 >> bcprop >> value;
568 
569             if (bcprop < 1 || bcprop > 4)
570               cerr << "bcprop out of range, bcprop = " << bcprop << endl;
571             p1++;
572             p2++;
573             p3++;
574             if (p1 < 1 || p1 > np || p2 < 1 || p2 > np || p3 < 1 || p3 > np)
575               {
576                 cout << "p1 = " << p1 << " p2 = " << p2 << " p3 = " << p3 << endl;
577               }
578 
579             if (i > 110354) Swap (p2, p3);
580             if (mesh.Point(p1)(0) < 0.25)
581               Swap (p2,p3);
582 
583             Element2d el(TRIG);
584 
585             if (bcprop == 1)
586               {
587                 if (values.Get(p1) < -69999)
588                   el.SetIndex(1);
589                 else
590                   el.SetIndex(2);
591               }
592             else
593               el.SetIndex(3);
594 
595 
596             el.PNum(1) = p1;
597             el.PNum(2) = p2;
598             el.PNum(3) = p3;
599             mesh.AddSurfaceElement (el);
600           }
601 
602 
603         ifstream incyl ("ngusers/guenter/cylinder.surf");
604         int npcyl, nsecyl;
605         incyl >> npcyl;
606         cout << "npcyl = " << npcyl << endl;
607         for (i = 1; i <= npcyl; i++)
608           {
609             Point3d p(0,0,0);
610             incyl >> p.X() >> p.Y() >> p.Z();
611             mesh.AddPoint (p);
612           }
613         incyl >> nsecyl;
614         cout << "nsecyl = " << nsecyl << endl;
615         for (i = 1; i <= nsecyl; i++)
616           {
617             incyl >> p1 >> p2 >> p3;
618             p1 += np;
619             p2 += np;
620             p3 += np;
621             Element2d el(TRIG);
622             el.SetIndex(5);
623             el.PNum(1) = p1;
624             el.PNum(2) = p2;
625             el.PNum(3) = p3;
626             mesh.AddSurfaceElement (el);
627           }
628       }
629 
630 
631     // .tet mesh
632     if ( (strlen (filename) > 4) &&
633          strcmp (&filename[strlen (filename)-4], ".tet") == 0 )
634       {
635         ReadTETFormat (mesh, filename);
636       }
637 
638 
639     // .fnf mesh (FNF - PE neutral format)
640     if ( (strlen (filename) > 4) &&
641          strcmp (&filename[strlen (filename)-4], ".fnf") == 0 )
642       {
643         ReadFNFFormat (mesh, filename);
644       }
645 
646     // .cgns file - CFD General Notation System
647     if ( (strlen (filename) > 5) &&
648          strcmp (&filename[strlen (filename)-5], ".cgns") == 0 )
649       {
650         ReadCGNSMesh (mesh, filename);
651       }
652 
653     if ( ( (strlen (filename) > 4) && strcmp (&filename[strlen (filename)-4], ".stl") == 0 ) ||
654          ( (strlen (filename) > 5) && strcmp (&filename[strlen (filename)-5], ".stlb") == 0 ) )
655       {
656         ifstream ist{string{filename}};
657         auto geom = shared_ptr<STLGeometry>(STLGeometry::Load(ist));
658 
659         mesh.SetDimension (3);
660 
661         auto & points = geom->GetPoints();
662 
663         for (auto & p : points)
664           mesh.AddPoint(MeshPoint(p));
665 
666         mesh.AddFaceDescriptor (FaceDescriptor (1, 1, 0, 1));
667 
668         for (auto ti : IntRange(geom->GetNT()))
669         {
670           Element2d el(TRIG);
671           for (auto i : IntRange(3))
672             el[i] = int((*geom)[STLTrigId(ti+IndexBASE<netgen::STLTrigId>())][i]);
673 
674           el.SetIndex(1);
675 
676           mesh.AddSurfaceElement(el);
677         }
678       }
679   }
680 
681 }
682 
683