1 //
2 // Write user dependent output file
3 //
4
5 #include <mystdlib.h>
6
7 #include <myadt.hpp>
8 #include <linalg.hpp>
9 #include <csg.hpp>
10 #include <geometry2d.hpp>
11 #include <meshing.hpp>
12
13 #include "writeuser.hpp"
14
15
16 namespace netgen
17 {
18 extern MeshingParameters mparam;
19
20
RegisterUserFormats(NgArray<const char * > & names,NgArray<const char * > & extensions)21 void RegisterUserFormats (NgArray<const char*> & names,
22 NgArray<const char*> & extensions)
23
24 {
25 const char *types[] =
26 {
27 "Neutral Format", ".mesh",
28 "Surface Mesh Format", ".mesh" ,
29 "DIFFPACK Format", ".mesh",
30 "TecPlot Format", ".mesh",
31 "Tochnog Format", ".mesh",
32 "Abaqus Format", ".mesh",
33 "Fluent Format", ".mesh",
34 "Permas Format", ".mesh",
35 "FEAP Format", ".mesh",
36 "Elmer Format", "*",
37 "STL Format", ".stl",
38 "STL Extended Format", ".stl",
39 "VRML Format", ".*",
40 "Gmsh Format", ".gmsh",
41 "Gmsh2 Format", ".gmsh2",
42 "OpenFOAM 1.5+ Format", "*",
43 "OpenFOAM 1.5+ Compressed", "*",
44 "JCMwave Format", ".jcm",
45 "TET Format", ".tet",
46 "CGNS Format", ".cgns",
47 // { "Chemnitz Format" },
48 0
49 };
50
51 for (int i = 0; types[2*i]; i++)
52 {
53 names.Append (types[2*i]);
54 extensions.Append (types[2*i+1]);
55 }
56 }
57
58
59
WriteUserFormat(const string & format,const Mesh & mesh,const string & filename)60 bool WriteUserFormat (const string & format,
61 const Mesh & mesh,
62 const string & filename)
63 {
64 // cout << "write user &hgeom = " << &hgeom << endl;
65 // const CSGeometry & geom = *dynamic_cast<const CSGeometry*> (&hgeom);
66 const CSGeometry & geom = *dynamic_pointer_cast<CSGeometry> (mesh.GetGeometry());
67
68 PrintMessage (1, "Export mesh to file ", filename,
69 ", format is ", format);
70
71 if (format == "Neutral Format")
72 WriteNeutralFormat (mesh, geom, filename);
73
74 else if (format == "Surface Mesh Format")
75 WriteSurfaceFormat (mesh, filename);
76
77 else if (format == "DIFFPACK Format")
78 WriteDiffPackFormat (mesh, geom, filename);
79
80 else if (format == "Tochnog Format")
81 WriteTochnogFormat (mesh, filename);
82
83 else if (format == "TecPlot Format")
84 cerr << "ERROR: TecPlot format currently out of order" << endl;
85 // WriteTecPlotFormat (mesh, geom, filename);
86
87 else if (format == "Abaqus Format")
88 WriteAbaqusFormat (mesh, filename);
89
90 else if (format == "Fluent Format")
91 WriteFluentFormat (mesh, filename);
92
93 else if (format == "Permas Format")
94 WritePermasFormat (mesh, filename);
95
96 else if (format == "FEAP Format")
97 WriteFEAPFormat (mesh, filename);
98
99 else if (format == "Elmer Format")
100 WriteElmerFormat (mesh, filename);
101
102 else if (format == "STL Format")
103 WriteSTLFormat (mesh, filename);
104
105 // Philippose - 16 August 2010
106 // Added additional STL Export in which
107 // each face of the geometry is treated
108 // as a separate "solid" entity
109 else if (format == "STL Extended Format")
110 WriteSTLExtFormat (mesh, filename);
111
112 else if (format == "VRML Format")
113 WriteVRMLFormat (mesh, 1, filename);
114
115 else if (format == "Fepp Format")
116 WriteFEPPFormat (mesh, geom, filename);
117
118 else if (format == "EdgeElement Format")
119 WriteEdgeElementFormat (mesh, geom, filename);
120
121 else if (format == "Chemnitz Format")
122 WriteUserChemnitz (mesh, filename);
123
124 else if (format == "Gmsh Format")
125 WriteGmshFormat (mesh, geom, filename);
126
127 // Philippose - 29/01/2009
128 // Added Gmsh v2.xx Mesh export capability
129 else if (format == "Gmsh2 Format")
130 WriteGmsh2Format (mesh, geom, filename);
131
132 // Philippose - 25/10/2009
133 // Added OpenFOAM 1.5+ Mesh export capability
134 else if (format == "OpenFOAM 1.5+ Format")
135 WriteOpenFOAM15xFormat (mesh, filename, false);
136
137 else if (format == "OpenFOAM 1.5+ Compressed")
138 WriteOpenFOAM15xFormat (mesh, filename, true);
139
140 else if (format == "JCMwave Format")
141 WriteJCMFormat (mesh, geom, filename);
142
143 #ifdef OLIVER
144 else if (format == "TET Format")
145 WriteTETFormat( mesh, filename);//, "High Frequency" );
146 #endif
147
148 else if (format == "CGNS Format")
149 WriteCGNSMesh( mesh, filename);
150
151 else
152 {
153 return 1;
154 }
155
156 return 0;
157 }
158
159
160
161
162 /*
163 * Neutral mesh format
164 * points, elements, surface elements
165 */
166
WriteNeutralFormat(const Mesh & mesh,const NetgenGeometry & geom,const string & filename)167 void WriteNeutralFormat (const Mesh & mesh,
168 const NetgenGeometry & geom,
169 const string & filename)
170 {
171 cout << "write neutral, new" << endl;
172 int np = mesh.GetNP();
173 int ne = mesh.GetNE();
174 int nse = mesh.GetNSE();
175 int nseg = mesh.GetNSeg();
176 int i, j;
177
178 int inverttets = mparam.inverttets;
179 int invertsurf = mparam.inverttrigs;
180
181 ofstream outfile (filename.c_str());
182
183 outfile.precision(6);
184 outfile.setf (ios::fixed, ios::floatfield);
185 outfile.setf (ios::showpoint);
186
187 outfile << np << "\n";
188
189 for (i = 1; i <= np; i++)
190 {
191 const Point3d & p = mesh.Point(i);
192
193 outfile.width(10);
194 outfile << p.X() << " ";
195 outfile.width(9);
196 outfile << p.Y() << " ";
197 if (mesh.GetDimension() == 3)
198 {
199 outfile.width(9);
200 outfile << p.Z();
201 }
202 outfile << "\n";
203 }
204
205 if (mesh.GetDimension() == 3)
206 {
207 outfile << ne << "\n";
208 for (i = 1; i <= ne; i++)
209 {
210 Element el = mesh.VolumeElement(i);
211 if (inverttets)
212 el.Invert();
213 outfile.width(4);
214 outfile << el.GetIndex() << " ";
215 for (j = 1; j <= el.GetNP(); j++)
216 {
217 outfile << " ";
218 outfile.width(8);
219 outfile << el.PNum(j);
220 }
221 outfile << "\n";
222 }
223 }
224
225 outfile << nse << "\n";
226 for (i = 1; i <= nse; i++)
227 {
228 Element2d el = mesh.SurfaceElement(i);
229 if (invertsurf)
230 el.Invert();
231 outfile.width(4);
232 outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
233 for (j = 1; j <= el.GetNP(); j++)
234 {
235 outfile << " ";
236 outfile.width(8);
237 outfile << el.PNum(j);
238 }
239 outfile << "\n";
240 }
241
242
243 if (mesh.GetDimension() == 2)
244 {
245 outfile << nseg << "\n";
246 for (int i = 1; i <= nseg; i++)
247 {
248 const Segment & seg = mesh.LineSegment(i);
249 outfile.width(4);
250 outfile << seg.si << " ";
251
252 for (int j = 0; j < seg.GetNP(); j++)
253 {
254 outfile << " ";
255 outfile.width(8);
256 outfile << seg[j];
257 }
258 /*
259 outfile << " ";
260 outfile.width(8);
261 outfile << seg[0];
262 outfile << " ";
263 outfile.width(8);
264 outfile << seg[1];
265 if (seg[2] != -1)
266 {
267 outfile.width(8);
268 outfile << seg[2];
269 }
270 */
271 outfile << "\n";
272 }
273 }
274 }
275
276
277
278
279
280
281
282
283
WriteSurfaceFormat(const Mesh & mesh,const string & filename)284 void WriteSurfaceFormat (const Mesh & mesh,
285 const string & filename)
286 {
287 // surface mesh
288 int i, j;
289
290 cout << "Write Surface Mesh" << endl;
291
292 ofstream outfile (filename.c_str());
293
294 outfile << "surfacemesh" << endl;
295
296 outfile << mesh.GetNP() << endl;
297 for (i = 1; i <= mesh.GetNP(); i++)
298 {
299 for (j = 0; j < 3; j++)
300 {
301 outfile.width(10);
302 outfile << mesh.Point(i)(j) << " ";
303 }
304 outfile << endl;
305 }
306 outfile << mesh.GetNSE() << endl;
307 for (i = 1; i <= mesh.GetNSE(); i++)
308 {
309 for (j = 1; j <= 3; j++)
310 {
311 outfile.width(8);
312 outfile << mesh.SurfaceElement(i).PNum(j);
313 }
314 outfile << endl;
315 }
316 }
317
318
319
320
321
322 /*
323 * save surface mesh as STL file
324 */
325
WriteSTLFormat(const Mesh & mesh,const string & filename)326 void WriteSTLFormat (const Mesh & mesh,
327 const string & filename)
328 {
329 cout << "\nWrite STL Surface Mesh" << endl;
330
331 ostream *outfile;
332
333 if(filename.substr(filename.length()-3,3) == ".gz")
334 outfile = new ogzstream(filename.c_str());
335 else
336 outfile = new ofstream(filename.c_str());
337
338 int i;
339
340 outfile->precision(10);
341
342 *outfile << "solid" << endl;
343
344 for (i = 1; i <= mesh.GetNSE(); i++)
345 {
346 *outfile << "facet normal ";
347 const Point3d& p1 = mesh.Point(mesh.SurfaceElement(i).PNum(1));
348 const Point3d& p2 = mesh.Point(mesh.SurfaceElement(i).PNum(2));
349 const Point3d& p3 = mesh.Point(mesh.SurfaceElement(i).PNum(3));
350
351 Vec3d normal = Cross(p2-p1,p3-p1);
352 if (normal.Length() != 0)
353 {
354 normal /= (normal.Length());
355 }
356
357 *outfile << normal.X() << " " << normal.Y() << " " << normal.Z() << "\n";
358 *outfile << "outer loop\n";
359
360 *outfile << "vertex " << p1.X() << " " << p1.Y() << " " << p1.Z() << "\n";
361 *outfile << "vertex " << p2.X() << " " << p2.Y() << " " << p2.Z() << "\n";
362 *outfile << "vertex " << p3.X() << " " << p3.Y() << " " << p3.Z() << "\n";
363
364 *outfile << "endloop\n";
365 *outfile << "endfacet\n";
366 }
367 *outfile << "endsolid" << endl;
368 }
369
370
371
372
373
374 /*
375 * Philippose - 16 August 2010
376 * Save surface mesh as STL file
377 * with a separate solid definition
378 * for each face
379 * - This helps in splitting up the
380 * STL into named boundary faces
381 * when using a third-party mesher
382 */
WriteSTLExtFormat(const Mesh & mesh,const string & filename)383 void WriteSTLExtFormat (const Mesh & mesh,
384 const string & filename)
385 {
386 cout << "\nWrite STL Surface Mesh (with separated boundary faces)" << endl;
387
388 ostream *outfile;
389
390 if(filename.substr(filename.length()-3,3) == ".gz")
391 outfile = new ogzstream(filename.c_str());
392 else
393 outfile = new ofstream(filename.c_str());
394
395 outfile->precision(10);
396
397 int numBCs = 0;
398
399 NgArray<int> faceBCs;
400 TABLE<int> faceBCMapping;
401
402 faceBCs.SetSize(mesh.GetNFD());
403 faceBCMapping.SetSize(mesh.GetNFD());
404
405 faceBCs = -1;
406
407 // Collect the BC numbers used in the mesh
408 for(int faceNr = 1; faceNr <= mesh.GetNFD(); faceNr++)
409 {
410 int bcNum = mesh.GetFaceDescriptor(faceNr).BCProperty();
411
412 if(faceBCs.Pos(bcNum) < 0)
413 {
414 numBCs++;
415 faceBCs.Set(numBCs,bcNum);
416 faceBCMapping.Add1(numBCs,faceNr);
417 }
418 else
419 {
420 faceBCMapping.Add1(faceBCs.Pos(bcNum)+1,faceNr);
421 }
422 }
423
424 faceBCs.SetSize(numBCs);
425 faceBCMapping.ChangeSize(numBCs);
426
427 // Now actually write the data to file
428 for(int bcInd = 1; bcInd <= faceBCs.Size(); bcInd++)
429 {
430 *outfile << "solid Boundary_" << faceBCs.Elem(bcInd) << "\n";
431
432 for(int faceNr = 1;faceNr <= faceBCMapping.EntrySize(bcInd); faceNr++)
433 {
434 Array<SurfaceElementIndex> faceSei;
435 mesh.GetSurfaceElementsOfFace(faceBCMapping.Get(bcInd,faceNr),faceSei);
436
437 for (int i = 0; i < faceSei.Size(); i++)
438 {
439 *outfile << "facet normal ";
440 const Point3d& p1 = mesh.Point(mesh[faceSei[i]].PNum(1));
441 const Point3d& p2 = mesh.Point(mesh[faceSei[i]].PNum(2));
442 const Point3d& p3 = mesh.Point(mesh[faceSei[i]].PNum(3));
443
444 Vec3d normal = Cross(p2-p1,p3-p1);
445 if (normal.Length() != 0)
446 {
447 normal /= (normal.Length());
448 }
449
450 *outfile << normal.X() << " " << normal.Y() << " " << normal.Z() << "\n";
451 *outfile << "outer loop\n";
452
453 *outfile << "vertex " << p1.X() << " " << p1.Y() << " " << p1.Z() << "\n";
454 *outfile << "vertex " << p2.X() << " " << p2.Y() << " " << p2.Z() << "\n";
455 *outfile << "vertex " << p3.X() << " " << p3.Y() << " " << p3.Z() << "\n";
456
457 *outfile << "endloop\n";
458 *outfile << "endfacet\n";
459 }
460 }
461 *outfile << "endsolid Boundary_" << faceBCs.Elem(bcInd) << "\n";
462 }
463 }
464
465
466
467
468 /*
469 *
470 * write surface mesh as VRML file
471 *
472 */
473
WriteVRMLFormat(const Mesh & mesh,bool faces,const string & filename)474 void WriteVRMLFormat (const Mesh & mesh,
475 bool faces,
476 const string & filename)
477 {
478
479 if (faces)
480
481 {
482 // Output in VRML, IndexedFaceSet is used
483 // Bartosz Sawicki <sawickib@ee.pw.edu.pl>
484
485 int np = mesh.GetNP();
486 int nse = mesh.GetNSE();
487 int i, j;
488
489 ofstream outfile (filename.c_str());
490
491 outfile.precision(6);
492 outfile.setf (ios::fixed, ios::floatfield);
493 outfile.setf (ios::showpoint);
494
495 outfile << "#VRML V2.0 utf8 \n"
496 "Background {\n"
497 " skyColor [1 1 1]\n"
498 " groundColor [1 1 1]\n"
499 "}\n"
500 "Group{ children [\n"
501 "Shape{ \n"
502 "appearance Appearance { material Material { }} \n"
503 "geometry IndexedFaceSet { \n"
504 "coord Coordinate { point [ \n";
505
506
507 for (i = 1; i <= np; i++)
508 {
509 const Point3d & p = mesh.Point(i);
510 outfile.width(10);
511 outfile << p.X() << " ";
512 outfile << p.Y() << " ";
513 outfile << p.Z() << " \n";
514 }
515
516 outfile << " ] } \n"
517 "coordIndex [ \n";
518
519 for (i = 1; i <= nse; i++)
520 {
521 const Element2d & el = mesh.SurfaceElement(i);
522
523 for (j = 1; j <= 3; j++)
524 {
525 outfile.width(8);
526 outfile << el.PNum(j)-1;
527 }
528 outfile << " -1 \n";
529 }
530
531 outfile << " ] \n";
532
533 //define number and RGB definitions of colors
534 outfile << "color Color { color [1 0 0, 0 1 0, 0 0 1, 1 1 0]} \n"
535 "colorIndex [\n";
536
537 for (i = 1; i <= nse; i++)
538 {
539 outfile << mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex ()).BCProperty();
540 outfile << endl;
541 }
542
543 outfile << " ] \n"
544 "colorPerVertex FALSE \n"
545 "creaseAngle 0 \n"
546 "solid FALSE \n"
547 "ccw FALSE \n"
548 "convex TRUE \n"
549 "} } # end of Shape\n"
550 "] }\n";
551
552 } /* end of VRMLFACES */
553
554
555 else
556
557 {
558 // Output in VRML, IndexedLineSet is used
559 // Bartosz Sawicki <sawickib@ee.pw.edu.pl>
560
561 int np = mesh.GetNP();
562 int nse = mesh.GetNSE();
563 int i, j;
564
565 ofstream outfile (filename.c_str());
566
567 outfile.precision(6);
568 outfile.setf (ios::fixed, ios::floatfield);
569 outfile.setf (ios::showpoint);
570
571 outfile << "#VRML V2.0 utf8 \n"
572 "Background {\n"
573 " skyColor [1 1 1]\n"
574 " groundColor [1 1 1]\n"
575 "}\n"
576 "Group{ children [\n"
577 "Shape{ \n"
578 "appearance Appearance { material Material { }} \n"
579 "geometry IndexedLineSet { \n"
580 "coord Coordinate { point [ \n";
581
582
583 for (i = 1; i <= np; i++)
584 {
585 const Point3d & p = mesh.Point(i);
586 outfile.width(10);
587 outfile << p.X() << " ";
588 outfile << p.Y() << " ";
589 outfile << p.Z() << " \n";
590 }
591
592 outfile << " ] } \n"
593 "coordIndex [ \n";
594
595 for (i = 1; i <= nse; i++)
596 {
597 const Element2d & el = mesh.SurfaceElement(i);
598
599 for (j = 1; j <= 3; j++)
600 {
601 outfile.width(8);
602 outfile << el.PNum(j)-1;
603 }
604 outfile.width(8);
605 outfile << el.PNum(1)-1;
606 outfile << " -1 \n";
607 }
608
609 outfile << " ] \n";
610
611 /* Uncomment if you want color mesh
612 outfile << "color Color { color [1 1 1, 0 1 0, 0 0 1, 1 1 0]} \n"
613 "colorIndex [\n";
614
615 for (i = 1; i <= nse; i++)
616 {
617 outfile << mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex ()).BCProperty();
618 outfile << endl;
619 }
620
621 outfile << " ] \n"
622 */
623 outfile << "colorPerVertex FALSE \n"
624 "} } #end of Shape\n"
625 "] } \n";
626
627 }
628
629 }
630
631
632
633
634
635
636 /*
637 * FEPP .. a finite element package developed at University Linz, Austria
638 */
WriteFEPPFormat(const Mesh & mesh,const NetgenGeometry & geom,const string & filename)639 void WriteFEPPFormat (const Mesh & mesh,
640 const NetgenGeometry & geom,
641 const string & filename)
642 {
643
644 ofstream outfile (filename.c_str());
645
646 if (mesh.GetDimension() == 3)
647
648 {
649
650 // output for FEPP
651
652 int np = mesh.GetNP();
653 int ne = mesh.GetNE();
654 int nse = mesh.GetNSE();
655 // int ns = mesh.GetNFD();
656 int i, j;
657
658 outfile.precision(5);
659 outfile.setf (ios::fixed, ios::floatfield);
660 outfile.setf (ios::showpoint);
661
662 outfile << "volumemesh4" << endl;
663 outfile << nse << endl;
664 for (i = 1; i <= nse; i++)
665 {
666 const Element2d & el = mesh.SurfaceElement(i);
667
668 // int facenr = mesh.facedecoding.Get(el.GetIndex()).surfnr;
669 outfile.width(4);
670 outfile << el.GetIndex() << " ";
671 outfile.width(4);
672 // outfile << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << " ";
673 outfile << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() << " ";
674 outfile.width(4);
675 outfile << el.GetNP() << " ";
676 for (j = 1; j <= el.GetNP(); j++)
677 {
678 outfile.width(8);
679 outfile << el.PNum(j);
680 }
681 outfile << "\n";
682 }
683
684
685 outfile << ne << "\n";
686 for (i = 1; i <= ne; i++)
687 {
688 const Element & el = mesh.VolumeElement(i);
689 outfile.width(4);
690 outfile << el.GetIndex() << " ";
691 outfile.width(4);
692 outfile << el.GetNP() << " ";
693 for (j = 1; j <= el.GetNP(); j++)
694 {
695 outfile.width(8);
696 outfile << el.PNum(j);
697 }
698 outfile << "\n";
699 }
700
701 outfile << np << "\n";
702 for (i = 1; i <= np; i++)
703 {
704 const Point3d & p = mesh.Point(i);
705
706 outfile.width(10);
707 outfile << p.X() << " ";
708 outfile.width(9);
709 outfile << p.Y() << " ";
710 outfile.width(9);
711 outfile << p.Z() << "\n";
712 }
713
714 /*
715 if (typ == WRITE_FEPPML)
716 {
717 int nbn = mesh.mlbetweennodes.Size();
718 outfile << nbn << "\n";
719 for (i = 1; i <= nbn; i++)
720 outfile << mesh.mlbetweennodes.Get(i).I1() << " "
721 << mesh.mlbetweennodes.Get(i).I2() << "\n";
722
723
724 // int ncon = mesh.connectedtonode.Size();
725 // outfile << ncon << "\n";
726 // for (i = 1; i <= ncon; i++)
727 // outfile << i << " " << mesh.connectedtonode.Get(i) << endl;
728 }
729 */
730
731 /*
732 // write CSG surfaces
733 if (&geom && geom.GetNSurf() >= ns)
734 {
735 outfile << ns << endl;
736 for (i = 1; i <= ns; i++)
737 geom.GetSurface(mesh.GetFaceDescriptor(i).SurfNr())->Print(outfile);
738 }
739 else
740 */
741 outfile << "0" << endl;
742 }
743
744
745 else
746
747 { // 2D fepp format
748
749 ;
750 /*
751 extern SplineGeometry2d * geometry2d;
752 if (geometry2d)
753 Save2DMesh (mesh, &geometry2d->GetSplines(), outfile);
754 else
755 Save2DMesh (mesh, 0, outfile);
756 */
757 }
758 }
759
760
761
762
763
764
765 /*
766 * Edge element mesh format
767 * points, elements, edges
768 */
769
WriteEdgeElementFormat(const Mesh & mesh,const NetgenGeometry & geom,const string & filename)770 void WriteEdgeElementFormat (const Mesh & mesh,
771 const NetgenGeometry & geom,
772 const string & filename)
773 {
774 cout << "write edge element format" << endl;
775
776 const MeshTopology * top = &mesh.GetTopology();
777 int npoints = mesh.GetNP();
778 int nelements = mesh.GetNE();
779 int nsurfelem = mesh.GetNSE();
780 int nedges = top->GetNEdges();
781 int i, j;
782
783 int inverttets = mparam.inverttets;
784 int invertsurf = mparam.inverttrigs;
785 NgArray<int> edges;
786
787 ofstream outfile (filename.c_str());
788
789 outfile.precision(6);
790 outfile.setf (ios::fixed, ios::floatfield);
791 outfile.setf (ios::showpoint);
792
793
794 // vertices with coordinates
795 outfile << npoints << "\n";
796 for (i = 1; i <= npoints; i++)
797 {
798 const Point3d & p = mesh.Point(i);
799
800 outfile.width(10);
801 outfile << p.X() << " ";
802 outfile.width(9);
803 outfile << p.Y() << " ";
804 outfile.width(9);
805 outfile << p.Z() << "\n";
806 }
807
808 // element - edge - list
809 outfile << nelements << " " << nedges << "\n";
810 for (i = 1; i <= nelements; i++)
811 {
812 Element el = mesh.VolumeElement(i);
813 if (inverttets)
814 el.Invert();
815 outfile.width(4);
816 outfile << el.GetIndex() << " ";
817 outfile.width(8);
818 outfile << el.GetNP();
819 for (j = 1; j <= el.GetNP(); j++)
820 {
821 outfile << " ";
822 outfile.width(8);
823 outfile << el.PNum(j);
824 }
825
826 top->GetElementEdges(i,edges);
827 outfile << endl << " ";
828 outfile.width(8);
829 outfile << edges.Size();
830 for (j=1; j <= edges.Size(); j++)
831 {
832 outfile << " ";
833 outfile.width(8);
834 outfile << edges[j-1];
835 }
836 outfile << "\n";
837
838 // orientation:
839 top->GetElementEdgeOrientations(i,edges);
840 outfile << " ";
841 for (j=1; j <= edges.Size(); j++)
842 {
843 outfile << " ";
844 outfile.width(8);
845 outfile << edges[j-1];
846 }
847 outfile << "\n";
848 }
849
850 // surface element - edge - list (with boundary conditions)
851 outfile << nsurfelem << "\n";
852 for (i = 1; i <= nsurfelem; i++)
853 {
854 Element2d el = mesh.SurfaceElement(i);
855 if (invertsurf)
856 el.Invert();
857 outfile.width(4);
858 outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << " ";
859 outfile.width(8);
860 outfile << el.GetNP();
861 for (j = 1; j <= el.GetNP(); j++)
862 {
863 outfile << " ";
864 outfile.width(8);
865 outfile << el.PNum(j);
866 }
867
868 top->GetSurfaceElementEdges(i,edges);
869 outfile << endl << " ";
870 outfile.width(8);
871 outfile << edges.Size();
872 for (j=1; j <= edges.Size(); j++)
873 {
874 outfile << " ";
875 outfile.width(8);
876 outfile << edges[j-1];
877 }
878 outfile << "\n";
879 }
880
881
882 int v1, v2;
883 // edge - vertex - list
884 outfile << nedges << "\n";
885 for (i=1; i <= nedges; i++)
886 {
887 top->GetEdgeVertices(i,v1,v2);
888 outfile.width(4);
889 outfile << v1;
890 outfile << " ";
891 outfile.width(8);
892 outfile << v2 << endl;
893 }
894 }
895
896
897
898
899
900
901
902
903
904 #ifdef OLDSTYLE_WRITE
905
906
WriteFile(int typ,const Mesh & mesh,const CSGeometry & geom,const char * filename,const char * geomfile,double h)907 void WriteFile (int typ,
908 const Mesh & mesh,
909 const CSGeometry & geom,
910 const char * filename,
911 const char * geomfile,
912 double h)
913 {
914
915
916 int inverttets = mparam.inverttets;
917 int invertsurf = mparam.inverttrigs;
918
919
920
921
922
923
924
925
926 if (typ == WRITE_EDGEELEMENT)
927 {
928 // write edge element file
929 // Peter Harscher, ETHZ
930
931 cout << "Write Edge-Element Format" << endl;
932
933 ofstream outfile (filename);
934
935 int i, j;
936 int ned;
937
938 // hash table representing edges;
939 INDEX_2_HASHTABLE<int> edgeht(mesh.GetNP());
940
941 // list of edges
942 NgArray<INDEX_2> edgelist;
943
944 // edge (point) on boundary ?
945 NgBitArray bedge, bpoint(mesh.GetNP());
946
947 static int eledges[6][2] = { { 1, 2 } , { 1, 3 } , { 1, 4 },
948 { 2, 3 } , { 2, 4 } , { 3, 4 } };
949
950 // fill hashtable (point1, point2) ----> edgenr
951 for (i = 1; i <= mesh.GetNE(); i++)
952 {
953 const Element & el = mesh.VolumeElement (i);
954 INDEX_2 edge;
955 for (j = 1; j <= 6; j++)
956 {
957 edge.I1() = el.PNum (eledges[j-1][0]);
958 edge.I2() = el.PNum (eledges[j-1][1]);
959 edge.Sort();
960
961 if (!edgeht.Used (edge))
962 {
963 edgelist.Append (edge);
964 edgeht.Set (edge, edgelist.Size());
965 }
966 }
967 }
968
969
970 // set bedges, bpoints
971 bedge.SetSize (edgelist.Size());
972 bedge.Clear();
973 bpoint.Clear();
974
975 for (i = 1; i <= mesh.GetNSE(); i++)
976 {
977 const Element2d & sel = mesh.SurfaceElement(i);
978 for (j = 1; j <= 3; j++)
979 {
980 bpoint.Set (sel.PNum(j));
981
982 INDEX_2 edge;
983 edge.I1() = sel.PNum(j);
984 edge.I2() = sel.PNum(j%3+1);
985 edge.Sort();
986
987 bedge.Set (edgeht.Get (edge));
988 }
989 }
990
991
992
993 outfile << mesh.GetNE() << endl;
994 // write element ---> point
995 for (i = 1; i <= mesh.GetNE(); i++)
996 {
997 const Element & el = mesh.VolumeElement(i);
998
999 outfile.width(8);
1000 outfile << i;
1001 for (j = 1; j <= 4; j++)
1002 {
1003 outfile.width(8);
1004 outfile << el.PNum(j);
1005 }
1006 outfile << endl;
1007 }
1008
1009 // write element ---> edge
1010 for (i = 1; i <= mesh.GetNE(); i++)
1011 {
1012 const Element & el = mesh.VolumeElement (i);
1013 INDEX_2 edge;
1014 for (j = 1; j <= 6; j++)
1015 {
1016 edge.I1() = el.PNum (eledges[j-1][0]);
1017 edge.I2() = el.PNum (eledges[j-1][1]);
1018 edge.Sort();
1019
1020 outfile.width(8);
1021 outfile << edgeht.Get (edge);
1022 }
1023 outfile << endl;
1024 }
1025
1026 // write points
1027 outfile << mesh.GetNP() << endl;
1028 outfile.precision (6);
1029 for (i = 1; i <= mesh.GetNP(); i++)
1030 {
1031 const Point3d & p = mesh.Point(i);
1032
1033 for (j = 1; j <= 3; j++)
1034 {
1035 outfile.width(8);
1036 outfile << p.X(j);
1037 }
1038 outfile << " "
1039 << (bpoint.Test(i) ? "1" : 0) << endl;
1040 }
1041
1042 // write edges
1043 outfile << edgelist.Size() << endl;
1044 for (i = 1; i <= edgelist.Size(); i++)
1045 {
1046 outfile.width(8);
1047 outfile << edgelist.Get(i).I1();
1048 outfile.width(8);
1049 outfile << edgelist.Get(i).I2();
1050 outfile << " "
1051 << (bedge.Test(i) ? "1" : "0") << endl;
1052 }
1053 }
1054
1055
1056
1057
1058 }
1059 #endif
1060 }
1061
1062