netgen/libsrc/interface/wuchemnitz.cpp

// Write Chemnitz file format


#include <mystdlib.h>

#include <myadt.hpp>

#include <linalg.hpp>
#include <csg.hpp>
#include <meshing.hpp>

namespace netgen
{

  class POINT3D
  {
  public:
    POINT3D () { };
    double x, y, z;
  };

  class VOLELEMENT
  {
  public:
    int domnr, p1, p2, p3, p4;
    int faces[4];

    VOLELEMENT () 
    { for (int i = 0; i < 4; i++) faces[i] = 0; }
  };
  
  class SURFELEMENT
  {
  public:
    SURFELEMENT () { };
    int snr, p1, p2, p3;
  };
  

  class FACE
  {
  public:
    int p1, p2, p3;
    int edges[3];

    FACE () 
    { for (int i = 0; i < 3; i++) edges[i] = 0; }
  };

  class EDGE
  {
  public:
    EDGE () { };
    int p1, p2;
  };

  static Array<POINT3D> points;
  static Array<VOLELEMENT> volelements;
  static Array<SURFELEMENT> surfelements;

  static Array<FACE> faces;
  static Array<EDGE> edges;


  void ReadFile (char * filename)
  {
    int i, n;
    ifstream infile(filename);
    char reco[100];
  
  
    infile >> reco;  // file format recognition
  
    infile >> n;   // number of surface elements
    cout << n << " Surface elements" << endl;
  
    for (i = 1; i <= n; i++)
      {
        SURFELEMENT sel;
        infile >> sel.snr >> sel.p1 >> sel.p2 >> sel.p3;
        surfelements.Append (sel);
      }
    
    infile >> n;   // number of volume elements
    cout << n << " Volume elements" << endl;
  
    for (i = 1; i <= n; i++)
      {
        VOLELEMENT el;
        infile >> el.p1 >> el.p2 >> el.p3 >> el.p4;
        volelements.Append (el);
      }
    
    infile >> n;   // number of points 
    cout << n << " Points" << endl;
  
    for (i = 1; i <= n; i++)
      {
        POINT3D p;
        infile >> p.x >> p.y >> p.z;
        points.Append (p);
      }
  }
  
  

  void ReadFileMesh (const Mesh & mesh)
  {
    int i, n;
  
    n = mesh.GetNSE();   // number of surface elements
    cout << n << " Surface elements" << endl;
  
    for (i = 1; i <= n; i++)
      {
        SURFELEMENT sel;
        const Element2d & el = mesh.SurfaceElement(i);
        sel.snr = el.GetIndex();
        sel.p1 = el.PNum(1);
        sel.p2 = el.PNum(2);
        sel.p3 = el.PNum(3);
        surfelements.Append (sel);
      }
    
    n = mesh.GetNE();   // number of volume elements
    cout << n << " Volume elements" << endl;
  
    for (i = 1; i <= n; i++)
      {
        VOLELEMENT el;
        const Element & nel = mesh.VolumeElement(i);
        el.p1 = nel.PNum(1);
        el.p2 = nel.PNum(2);
        el.p3 = nel.PNum(3);
        el.p4 = nel.PNum(4);
        //      infile >> el.p1 >> el.p2 >> el.p3 >> el.p4;
        volelements.Append (el);
      }
    
    n = mesh.GetNP();   // number of points 
    cout << n << " Points" << endl;
  
    for (i = 1; i <= n; i++)
      {
        POINT3D p;
        Point3d mp = mesh.Point(i);
        p.x = mp.X();
        p.y = mp.Y();
        p.z = mp.Z();
        //      infile >> p.x >> p.y >> p.z;
        points.Append (p);
      }
  }
  



  void Convert ()
  {
    int i, j, facei, edgei;
    INDEX_3 i3;
    INDEX_2 i2;

    INDEX_3_HASHTABLE<int> faceindex(volelements.Size()/5 + 1);
    INDEX_2_HASHTABLE<int> edgeindex(volelements.Size()/5 + 1);
  
    for (i = 1; i <= volelements.Size(); i++)
      {
        for (j = 1; j <= 4; j++)
          {
            switch (j)
              {
              case 1:
                i3.I1() = volelements.Get(i).p2;
                i3.I2() = volelements.Get(i).p3;
                i3.I3() = volelements.Get(i).p4;
                break;
              case 2:
                i3.I1() = volelements.Get(i).p1;
                i3.I2() = volelements.Get(i).p3;
                i3.I3() = volelements.Get(i).p4;
                break;
              case 3:
                i3.I1() = volelements.Get(i).p1;
                i3.I2() = volelements.Get(i).p2;
                i3.I3() = volelements.Get(i).p4;
                break;
              case 4:
                i3.I1() = volelements.Get(i).p1;
                i3.I2() = volelements.Get(i).p2;
                i3.I3() = volelements.Get(i).p3;
                break;
              default:
                i3.I1()=i3.I2()=i3.I3()=0;
              }
            i3.Sort();
            if (faceindex.Used (i3)) 
              facei = faceindex.Get(i3);
            else
              {
                FACE fa;
                fa.p1 = i3.I1();
                fa.p2 = i3.I2();
                fa.p3 = i3.I3();
                faces.Append (fa);
                facei = faces.Size();
                faceindex.Set (i3, facei);
              } 
        
            volelements.Elem(i).faces[j-1] = facei;  
          }    
    
      } 
 

    for (i = 1; i <= faces.Size(); i++)
      {
        for (j = 1; j <= 3; j++)
          {
            switch (j)
              {
              case 1:
                i2.I1() = faces.Get(i).p2;
                i2.I2() = faces.Get(i).p3;
                break;
              case 2:
                i2.I1() = faces.Get(i).p1;
                i2.I2() = faces.Get(i).p3;
                break;
              case 3:
                i2.I1() = faces.Get(i).p1;
                i2.I2() = faces.Get(i).p2;
                break;
              default:
                i2.I1()=i2.I2()=0;
              }
            if (i2.I1() > i2.I2()) swap (i2.I1(), i2.I2());
            if (edgeindex.Used (i2)) 
              edgei = edgeindex.Get(i2);
            else
              {
                EDGE ed;
                ed.p1 = i2.I1();
                ed.p2 = i2.I2();
                edges.Append (ed);
                edgei = edges.Size();
                edgeindex.Set (i2, edgei);
              } 
        
            faces.Elem(i).edges[j-1] = edgei;  
          }    
    
      }  
 
  }  
  
  
  void WriteFile (ostream & outfile)
  {
    int i;
  
    outfile 
      << "#VERSION: 1.0" << endl
      << "#PROGRAM: NETGEN" << endl
      << "#EQN_TYPE: POISSON" << endl
      << "#DIMENSION: 3D" << endl
      << "#DEG_OF_FREE: 1" << endl
      << "#DESCRIPTION: I don't know" << endl
      << "##RENUM: not done" << endl
      << "#USER: Kleinzen" << endl
      << "DATE: 10.06.1996" << endl;
  
    outfile << "#HEADER:   8" << endl
            << points.Size() << "  " << edges.Size() << "  " 
            << faces.Size() << "  " << volelements.Size() << "  0  0  0  0" << endl;
  
    outfile << "#VERTEX:   " << points.Size() << endl;
    for (i = 1; i <= points.Size(); i++)
      outfile << "  " << i << "  " << points.Get(i).x << "  " << points.Get(i).y 
              << "  " << points.Get(i).z << endl;
    	
    outfile << "#EDGE:  " << edges.Size() << endl;
    for (i = 1; i <= edges.Size(); i++)
      outfile << "  " << i << "  1  " 
              << edges.Get(i).p1 << "  " 
              << edges.Get(i).p2 
              << "  0" << endl;
    
    outfile << "#FACE:  " << faces.Size() << endl;  
    for (i = 1; i <= faces.Size(); i++)
      outfile << "  " << i << "  1  3  " 
              << faces.Get(i).edges[0] << "  " 
              << faces.Get(i).edges[1] << "  " 
              << faces.Get(i).edges[2] << endl;
    	
    outfile << "#SOLID:  " << volelements.Size() << endl;
    for (i = 1; i <= volelements.Size(); i++)
      outfile << "  " << i << "  1  4  " 
              << volelements.Get(i).faces[0] << "  "
              << volelements.Get(i).faces[1] << "  "
              << volelements.Get(i).faces[2] << "  "
              << volelements.Get(i).faces[3] << endl;
    	
    outfile << "#END_OF_DATA" << endl;
  }
    

  void WriteUserChemnitz (const Mesh & mesh,
                          const string & filename)
  {
    ofstream outfile (filename.c_str());

    ReadFileMesh (mesh);
    Convert ();
  
    WriteFile (outfile);
    cout << "Wrote Chemnitz standard file" << endl;
  }
}