ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-11 16:49:16 +05:00
Code Issues Packages Projects Releases Wiki Activity
9f32a5c3ad
netgen/libsrc/interface/writejcm.cpp
Matthias Hochsteger cb87362f64 Rename Array to NgArray
2019-07-09 10:39:16 +02:00

431 lines
13 KiB
C++
Raw Blame History

//
// Write JCMwave file
// 07.07.2005, Sven Burger, ZIB Berlin
//
#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <csg.hpp>
#include <meshing.hpp>
#include <sys/stat.h>
namespace netgen
{
#include "writeuser.hpp"
void WriteJCMFormat (const Mesh & mesh,
const NetgenGeometry & geom,
const string & filename)
{
if (mesh.GetDimension() != 3)
{
cout <<"\n Error: Dimension 3 only supported by this output format!"<<endl;
return;
}
int bc_at_infinity = 0;
int i, j, jj, ct(0), counter;
double dx1, dx2, dx3, dy1, dy2, dy3, dz1, dz2, dz3, vol;
// number of points
int np = mesh.GetNP();
// Identic points
NgArray<int,PointIndex::BASE> identmap1, identmap2, identmap3;
mesh.GetIdentifications().GetMap(1, identmap1);
mesh.GetIdentifications().GetMap(2, identmap2);
mesh.GetIdentifications().GetMap(3, identmap3);
// number of volume elements
int ne = mesh.GetNE();
int ntets = 0;
int nprisms = 0;
for (i = 1; i <= ne; i++)
{
Element el = mesh.VolumeElement(i);
if (el.GetNP() == 4)
{
ntets++;
// Check that no two points on a tetrahedron are identified with each other
for (j = 1; j <= 4; j++)
for (jj = 1; jj <=4; jj++)
{
if (identmap1.Elem(el.PNum(j)) == el.PNum(jj))
{
cout << "\n Error: two points on a tetrahedron identified (1) with each other"
<< "\n REFINE MESH !" << endl;
return;
}
if (identmap2.Elem(el.PNum(j)) == el.PNum(jj))
{
cout << "\n Error: two points on a tetrahedron identified (2) with each other"
<< "\n REFINE MESH !" << endl;
return;
}
if (identmap3.Elem(el.PNum(j)) == el.PNum(jj))
{
cout << "\n Error: two points on a tetrahedron identified (3) with each other"
<< "\n REFINE MESH !" << endl;
return;
}
}
}
else if (el.GetNP() == 6)
nprisms++;
}
if ( ne != (ntets+nprisms))
{
cout<< "\n Error in determining number of volume elements!\n"
<< "\n Prisms and tetrahedra only implemented in the JCMwave format!\n"<<endl;
return;
}
if (nprisms > 0)
cout << " Please note: Boundaries at infinity have to carry the bc-attribute '-bc="
<< bc_at_infinity <<"'."<<endl;
// number of surface elements
int nse = mesh.GetNSE();
// number of boundary triangles
int nbtri = 0;
// number of boundary quadrilaterals
int nbquad = 0;
// array with 1 if point on any tetra, 0 else
// this is needed in order to arrange the prism points in the right order
NgArray<int,1> pointsOnTetras;
pointsOnTetras.SetSize (mesh.GetNP());
pointsOnTetras = 0;
for (i = 1; i <= ne; i++)
{
Element el = mesh.VolumeElement(i);
if (el.GetNP() == 4)
{
for (j = 1; j <= 4; j++)
pointsOnTetras.Set(int (el.PNum(j)),1);
}
}
// number of boundary triangles and boundary quadrilaterals
for (i = 1; i <= nse; i++)
{
Element2d el = mesh.SurfaceElement(i);
if (el.GetNP() == 3 &&
( mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0 ||
mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0 ) )
nbtri++;
else if (el.GetNP() == 4 &&
( mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0 ||
mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0 ) )
nbquad++;
}
ofstream outfile (filename.c_str());
outfile.precision(6);
outfile.setf (ios::fixed, ios::floatfield);
outfile.setf (ios::showpoint);
outfile << "/* <BLOBHead>\n";
outfile << "__BLOBTYPE__=Grid\n";
outfile << "__OWNER__=JCMwave\n";
outfile << "<I>SpaceDim=3\n";
outfile << "<I>ManifoldDim=3\n";
outfile << "<I>NRefinementSteps=0\n";
outfile << "<I>NPoints="<<np<<"\n";
outfile << "<I>NTetrahedra="<<ntets<<"\n";
outfile << "<I>NPrisms="<<nprisms<<"\n";
outfile << "<I>NBoundaryTriangles="<<nbtri<<"\n";
outfile << "<I>NBoundaryQuadrilaterals="<<nbquad<<"\n";
outfile << "*/\n";
outfile << "\n";
outfile << "# output from Netgen\n\n";
int nDomains=mesh.GetNDomains();
for (i=1; i<=nDomains; i++)
{
if (mesh.GetMaterialPtr(i))
outfile << "#" << mesh.GetMaterial(i)
<< ": Material ID = "
<< i << "\n";
}
outfile << "# Points\n";
cout << " Please note: The unit of length in the .geo file is assumed to be 'microns'."<<endl;
for (i = 1; i <= np; i++)
{
const Point<3> & p = mesh.Point(i);
outfile << i << "\n";
outfile << p(0) << "e-6\n";
outfile << p(1) << "e-6\n";
outfile << p(2) << "e-6\n\n";
}
outfile << "\n";
outfile << "# Tetrahedra\n";
counter = 0;
for (i = 1; i <= ne; i++)
{
Element el = mesh.VolumeElement(i);
if (el.GetNP() == 4)
{
counter++;
dx1 = mesh.Point(el.PNum(2))(0) - mesh.Point(el.PNum(1))(0);
dx2 = mesh.Point(el.PNum(3))(0) - mesh.Point(el.PNum(1))(0);
dx3 = mesh.Point(el.PNum(4))(0) - mesh.Point(el.PNum(1))(0);
dy1 = mesh.Point(el.PNum(2))(1) - mesh.Point(el.PNum(1))(1);
dy2 = mesh.Point(el.PNum(3))(1) - mesh.Point(el.PNum(1))(1);
dy3 = mesh.Point(el.PNum(4))(1) - mesh.Point(el.PNum(1))(1);
dz1 = mesh.Point(el.PNum(2))(2) - mesh.Point(el.PNum(1))(2);
dz2 = mesh.Point(el.PNum(3))(2) - mesh.Point(el.PNum(1))(2);
dz3 = mesh.Point(el.PNum(4))(2) - mesh.Point(el.PNum(1))(2);
vol = (dy1*dz2-dz1*dy2)*dx3 + (dz1*dx2-dx1*dz2)*dy3 + (dx1*dy2-dy1*dx2)*dz3;
if ( vol > 0 )
for (j = 1; j <= 4; j++)
outfile << el.PNum(j)<<"\n";
else
{
for (j = 2; j >= 1; j--)
outfile << el.PNum(j)<<"\n";
for (j = 3; j <= 4; j++)
outfile << el.PNum(j)<<"\n";
}
outfile << el.GetIndex() << "\n\n";
}
}
if ( counter != ntets)
{
cout<< "\n Error in determining number of tetras!\n"<<endl;
return;
}
outfile << "\n";
outfile << "# Prisms\n";
counter = 0;
for (i = 1; i <= ne; i++)
{
Element el = mesh.VolumeElement(i);
if (el.GetNP() == 6)
{
counter++;
dx1 = mesh.Point(el.PNum(2))(0) - mesh.Point(el.PNum(1))(0);
dx2 = mesh.Point(el.PNum(3))(0) - mesh.Point(el.PNum(1))(0);
dx3 = mesh.Point(el.PNum(4))(0) - mesh.Point(el.PNum(1))(0);
dy1 = mesh.Point(el.PNum(2))(1) - mesh.Point(el.PNum(1))(1);
dy2 = mesh.Point(el.PNum(3))(1) - mesh.Point(el.PNum(1))(1);
dy3 = mesh.Point(el.PNum(4))(1) - mesh.Point(el.PNum(1))(1);
dz1 = mesh.Point(el.PNum(2))(2) - mesh.Point(el.PNum(1))(2);
dz2 = mesh.Point(el.PNum(3))(2) - mesh.Point(el.PNum(1))(2);
dz3 = mesh.Point(el.PNum(4))(2) - mesh.Point(el.PNum(1))(2);
vol = (dy1*dz2-dz1*dy2)*dx3 + (dz1*dx2-dx1*dz2)*dy3 + (dx1*dy2-dy1*dx2)*dz3;
if (pointsOnTetras.Get(el.PNum(1)) &&
pointsOnTetras.Get(el.PNum(2)) &&
pointsOnTetras.Get(el.PNum(3)))
{
if (vol > 0)
for (j = 1; j <= 6; j++)
outfile << el.PNum(j)<<"\n";
else
{
for (j = 3; j >= 1; j--)
outfile << el.PNum(j)<<"\n";
for (j = 6; j >= 4; j--)
outfile << el.PNum(j)<<"\n";
}
}
else if ( pointsOnTetras.Get(el.PNum(4)) &&
pointsOnTetras.Get(el.PNum(5)) &&
pointsOnTetras.Get(el.PNum(6)) )
{
if ( vol < 0 )
{
for (j = 4; j <= 6; j++)
outfile << el.PNum(j)<<"\n";
for (j = 1; j <= 3; j++)
outfile << el.PNum(j)<<"\n";
}
else
{
for (j = 6; j >= 4; j--)
outfile << el.PNum(j)<<"\n";
for (j = 3; j >= 1; j--)
outfile << el.PNum(j)<<"\n";
}
}
else
{
cout << "\n Error in determining prism point numbering!\n"<<endl;
return;
}
outfile << el.GetIndex() << "\n\n";
}
}
if ( counter != nprisms)
{
cout<< "\n Error in determining number of prisms!\n"<<endl;
return;
}
int npid1 = 0;
int npid2 = 0;
int npid3 = 0;
for (i=1; i<=np; i++)
{
if (identmap1.Elem(i))
npid1++;
if (identmap2.Elem(i))
npid2++;
if (identmap3.Elem(i))
npid3++;
}
outfile << "\n";
outfile << "# Boundary triangles\n";
outfile << "# Number of identified points in 1-direction: " << npid1 << "\n";
outfile << "# Number of identified points in 2-direction: " << npid2 << "\n";
outfile << "# Number of identified points in 3-direction: " << npid3 << "\n";
for (i = 1; i <= nse; i++)
{
Element2d el = mesh.SurfaceElement(i);
if (el.GetNP() == 3
&& (mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0
|| mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0))
{
outfile <<"# T\n";
for (j = 1; j <= 3; j++)
outfile << el.PNum(j)<<"\n";
if (mesh.GetFaceDescriptor (el.GetIndex()).BCProperty()==bc_at_infinity)
outfile << 1000 << "\n";
else
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << "\n";
if (mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() == bc_at_infinity)
outfile << "-2\n\n";
else if (identmap1.Elem(el.PNum(1))
&&identmap1.Elem(el.PNum(2))
&&identmap1.Elem(el.PNum(3)))
{
outfile << "-1\n";
for (j = 1; j <= 3; j++)
outfile << identmap1.Elem(el.PNum(j))<<"\n";
outfile << "\n";
}
else if (identmap2.Elem(el.PNum(1))
&&identmap2.Elem(el.PNum(2))
&&identmap2.Elem(el.PNum(3)))
{
outfile << "-1\n";
for (j = 1; j <= 3; j++)
outfile << identmap2.Elem(el.PNum(j))<<"\n";
outfile << "\n";
}
else if (identmap3.Elem(el.PNum(1))
&&identmap3.Elem(el.PNum(2))
&&identmap3.Elem(el.PNum(3)))
{
outfile << "-1\n";
for (j = 1; j <= 3; j++)
outfile << identmap3.Elem(el.PNum(j))<<"\n";
outfile << "\n";
}
else
outfile << "1\n\n";
}
}
outfile << "\n";
outfile << "# Boundary quadrilaterals\n";
for (i = 1; i <= nse; i++)
{
Element2d el = mesh.SurfaceElement(i);
if (el.GetNP() == 4
&& (mesh.GetFaceDescriptor (el.GetIndex()).DomainIn()==0
|| mesh.GetFaceDescriptor (el.GetIndex()).DomainOut()==0))
{
if (pointsOnTetras.Get(el.PNum(1)) &&
pointsOnTetras.Get(el.PNum(2)))
ct = 0;
else if (pointsOnTetras.Get(el.PNum(2)) &&
pointsOnTetras.Get(el.PNum(3)))
ct = 1;
else if (pointsOnTetras.Get(el.PNum(3)) &&
pointsOnTetras.Get(el.PNum(4)))
ct = 2;
else if (pointsOnTetras.Get(el.PNum(4)) &&
pointsOnTetras.Get(el.PNum(1)))
ct = 3;
else
cout << "\nWarning: Quadrilateral with inconsistent points found!"<<endl;
for (j = 1; j <= 4; j++)
{
jj = j + ct;
if ( jj >= 5 )
jj = jj - 4;
outfile << el.PNum(jj)<<"\n";
}
outfile << mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() << "\n";
if (mesh.GetFaceDescriptor (el.GetIndex()).BCProperty() == bc_at_infinity)
{
outfile << "-2\n\n";
cout << "\nWarning: Quadrilateral at infinity found (this should not occur)!"<<endl;
}
else if ( identmap1.Elem(el.PNum(1)) &&
identmap1.Elem(el.PNum(2)) &&
identmap1.Elem(el.PNum(3)) &&
identmap1.Elem(el.PNum(4)) )
{
outfile << "-1\n";
for (j = 1; j <= 4; j++)
{
jj = j + ct;
if ( jj >= 5 )
jj = jj - 4;
outfile << identmap1.Elem(el.PNum(jj))<<"\n";
}
outfile << "\n";
}
else if ( identmap2.Elem(el.PNum(1)) &&
identmap2.Elem(el.PNum(2)) &&
identmap2.Elem(el.PNum(3)) &&
identmap2.Elem(el.PNum(4)) )
{
outfile << "-1\n";
for (j = 1; j <= 4; j++)
{
jj = j + ct;
if ( jj >= 5 )
jj = jj - 4;
outfile << identmap2.Elem(el.PNum(jj))<<"\n";
}
outfile << "\n";
}
else if ( identmap3.Elem(el.PNum(1)) &&
identmap3.Elem(el.PNum(2)) &&
identmap3.Elem(el.PNum(3)) &&
identmap3.Elem(el.PNum(4)) )
{
outfile << "-1\n";
for (j = 1; j <= 4; j++)
{
jj = j + ct;
if ( jj >= 5 )
jj = jj - 4;
outfile << identmap3.Elem(el.PNum(jj))<<"\n";
}
outfile << "\n";
}
else
outfile << "1\n\n";
}
}
cout << " JCMwave grid file written." << endl;
}
}
Reference in New Issue View Git Blame Copy Permalink