mirror of
https://github.com/NGSolve/netgen.git
synced 2024-11-14 10:08:32 +05:00
326 lines
8.6 KiB
C++
326 lines
8.6 KiB
C++
//
|
|
// Write diffpack file
|
|
//
|
|
// by
|
|
// Bartosz Sawicki <sawickib@ee.pw.edu.pl>
|
|
// extended by
|
|
// Jacques Lechelle <jacques.lechelle@wanadoo.fr>
|
|
//
|
|
#include <mystdlib.h>
|
|
|
|
#include <myadt.hpp>
|
|
#include <linalg.hpp>
|
|
#include <csg.hpp>
|
|
#include <meshing.hpp>
|
|
|
|
|
|
namespace netgen
|
|
{
|
|
#include "writeuser.hpp"
|
|
|
|
|
|
void WriteDiffPackFormat (const Mesh & mesh,
|
|
const CSGeometry & geom,
|
|
const string & filename)
|
|
{
|
|
// double scale = globflags.GetNumFlag ("scale", 1);
|
|
double scale = 1;
|
|
|
|
ofstream outfile(filename.c_str());
|
|
outfile.precision(14);
|
|
|
|
|
|
if (mesh.GetDimension() == 3)
|
|
|
|
{
|
|
// Output compatible to Diffpack grid format
|
|
// Bartosz Sawicki <sawickib@ee.pw.edu.pl>
|
|
|
|
int np = mesh.GetNP();
|
|
int ne = mesh.GetNE();
|
|
int nse = mesh.GetNSE();
|
|
Array <int> BIname;
|
|
Array <int> BCsinpoint;
|
|
int i, j, k, l;
|
|
|
|
|
|
outfile.precision(6);
|
|
outfile.setf (ios::fixed, ios::floatfield);
|
|
outfile.setf (ios::showpoint);
|
|
|
|
const Element & eldummy = mesh.VolumeElement((int)1);
|
|
outfile << "\n\n"
|
|
"Finite element mesh (GridFE):\n\n"
|
|
" Number of space dim. = 3\n"
|
|
" Number of elements = " << ne << "\n"
|
|
" Number of nodes = " << np << "\n\n"
|
|
" All elements are of the same type : dpTRUE\n"
|
|
" Max number of nodes in an element: "<< eldummy.GetNP() << "\n"
|
|
" Only one subdomain : dpFALSE\n"
|
|
" Lattice data ? 0\n\n\n\n";
|
|
|
|
for (i = 1; i <= nse; i++)
|
|
{
|
|
int BI=mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex()).BCProperty();
|
|
int nbi=BIname.Size();
|
|
int found=0;
|
|
for (j = 1; j <= nbi; j++)
|
|
if(BI == BIname.Get(j)) found = 1;
|
|
if( ! found ) BIname.Append(BI);
|
|
}
|
|
|
|
outfile << " " << BIname.Size() << " Boundary indicators: ";
|
|
for (i =1 ; i <= BIname.Size(); i++)
|
|
outfile << BIname.Get(i) << " ";
|
|
outfile << "\n\n\n";
|
|
|
|
outfile << " Nodal coordinates and nodal boundary indicators,\n"
|
|
" the columns contain:\n"
|
|
" - node number\n"
|
|
" - coordinates\n"
|
|
" - no of boundary indicators that are set (ON)\n"
|
|
" - the boundary indicators that are set (ON) if any.\n"
|
|
"#\n";
|
|
|
|
|
|
// setup point-to-surfaceelement table
|
|
TABLE<SurfaceElementIndex, PointIndex::BASE> point2sel(np);
|
|
for (SurfaceElementIndex sei = 0; sei < nse; sei++)
|
|
{
|
|
const Element2d & el = mesh[sei];
|
|
for (int j = 0; j < el.GetNP(); j++)
|
|
point2sel.Add (el[j], sei);
|
|
}
|
|
|
|
|
|
for (i = 1; i <= np; i++)
|
|
{
|
|
const Point3d & p = mesh.Point(i);
|
|
|
|
outfile.width(12);
|
|
outfile << i << " (";
|
|
outfile.width(16);
|
|
outfile << p.X()/scale << ", ";
|
|
outfile.width(16);
|
|
outfile << p.Y()/scale << ", ";
|
|
outfile.width(16);
|
|
outfile << p.Z()/scale << ") ";
|
|
|
|
if(mesh[PointIndex(i)].Type() != INNERPOINT)
|
|
{
|
|
BCsinpoint.DeleteAll();
|
|
/*
|
|
for (j = 1; j <= nse; j++)
|
|
*/
|
|
FlatArray<SurfaceElementIndex> sels = point2sel[i];
|
|
for (int jj = 0; jj < sels.Size(); jj++)
|
|
{
|
|
for (k = 1; k <= mesh[sels[jj]].GetNP(); k++)
|
|
{
|
|
if(mesh[sels[jj]].PNum(k)==i)
|
|
{
|
|
int BC=mesh.GetFaceDescriptor(mesh[sels[jj]].GetIndex()).BCProperty();
|
|
int nbcsp=BCsinpoint.Size();
|
|
int found = 0;
|
|
for (l = 1; l <= nbcsp; l++)
|
|
if(BC == BCsinpoint.Get(l)) found = 1;
|
|
if( ! found ) BCsinpoint.Append(BC);
|
|
}
|
|
}
|
|
}
|
|
int nbcsp = BCsinpoint.Size();
|
|
outfile << "[" << nbcsp << "] ";
|
|
for (j = 1; j <= nbcsp; j++)
|
|
outfile << BCsinpoint.Get(j) << " ";
|
|
outfile << "\n";
|
|
}
|
|
else outfile << "[0]\n";
|
|
|
|
}
|
|
|
|
outfile << "\n"
|
|
" Element types and connectivity\n"
|
|
" the columns contain:\n"
|
|
" - element number\n"
|
|
" - element type\n"
|
|
" - subdomain number\n"
|
|
" - the global node numbers of the nodes in the element.\n"
|
|
"#\n";
|
|
|
|
for (i = 1; i <= ne; i++)
|
|
{
|
|
const Element & el = mesh.VolumeElement(i);
|
|
outfile.width(5);
|
|
if(el.GetNP()==4)
|
|
outfile << i << " ElmT4n3D ";
|
|
else
|
|
outfile << i << " ElmT10n3D ";
|
|
outfile.width(4);
|
|
outfile << el.GetIndex() << " ";
|
|
if(el.GetNP()==10)
|
|
{
|
|
outfile.width(8);
|
|
outfile << el.PNum(1);
|
|
outfile.width(8);
|
|
outfile << el.PNum(3);
|
|
outfile.width(8);
|
|
outfile << el.PNum(2);
|
|
outfile.width(8);
|
|
outfile << el.PNum(4);
|
|
outfile.width(8);
|
|
outfile << el.PNum(6);
|
|
outfile.width(8);
|
|
outfile << el.PNum(8);
|
|
outfile.width(8);
|
|
outfile << el.PNum(5);
|
|
outfile.width(8);
|
|
outfile << el.PNum(7);
|
|
outfile.width(8);
|
|
outfile << el.PNum(10);
|
|
outfile.width(8);
|
|
outfile << el.PNum(9);
|
|
}
|
|
else
|
|
{
|
|
outfile.width(8);
|
|
outfile << el.PNum(1);
|
|
outfile.width(8);
|
|
outfile << el.PNum(3);
|
|
outfile.width(8);
|
|
outfile << el.PNum(2);
|
|
outfile.width(8);
|
|
outfile << el.PNum(4);
|
|
}
|
|
outfile << "\n";
|
|
}
|
|
} /* Diffpack */
|
|
|
|
else
|
|
|
|
{
|
|
// Output compatible to Diffpack grid format 2D
|
|
|
|
int np = mesh.GetNP();
|
|
//int ne = mesh.GetNE();
|
|
int nse = mesh.GetNSE();
|
|
Array <int> BIname;
|
|
Array <int> BCsinpoint;
|
|
int i, j, k, l;
|
|
|
|
|
|
outfile.precision(6);
|
|
outfile.setf (ios::fixed, ios::floatfield);
|
|
outfile.setf (ios::showpoint);
|
|
const Element2d & eldummy = mesh.SurfaceElement((int)1);
|
|
outfile << "\n\n"
|
|
"Finite element mesh (GridFE):\n\n"
|
|
" Number of space dim. = 2\n"
|
|
" Number of elements = " << nse << "\n"
|
|
" Number of nodes = " << np << "\n\n"
|
|
" All elements are of the same type : dpTRUE\n"
|
|
" Max number of nodes in an element: "<<eldummy.GetNP()<<"\n"
|
|
" Only one subdomain : dpFALSE\n"
|
|
" Lattice data ? 0\n\n\n\n";
|
|
|
|
for (i = 1; i <= nse; i++)
|
|
{
|
|
int BI=mesh.GetFaceDescriptor(mesh.SurfaceElement(i).GetIndex()).BCProperty();
|
|
int nbi=BIname.Size();
|
|
int found=0;
|
|
for (j = 1; j <= nbi; j++)
|
|
if(BI == BIname.Get(j)) found = 1;
|
|
if( ! found ) BIname.Append(BI);
|
|
}
|
|
|
|
outfile << " " << BIname.Size() << " Boundary indicators: ";
|
|
for (i =1 ; i <= BIname.Size(); i++)
|
|
outfile << BIname.Get(i) << " ";
|
|
outfile << "\n\n\n";
|
|
|
|
outfile << " Nodal coordinates and nodal boundary indicators,\n"
|
|
" the columns contain:\n"
|
|
" - node number\n"
|
|
" - coordinates\n"
|
|
" - no of boundary indicators that are set (ON)\n"
|
|
" - the boundary indicators that are set (ON) if any.\n"
|
|
"#\n";
|
|
|
|
for (i = 1; i <= np; i++)
|
|
{
|
|
const Point3d & p = mesh.Point(i);
|
|
|
|
outfile.width(12);
|
|
outfile << i << " (";
|
|
outfile.width(16);
|
|
outfile << p.X()/scale << ", ";
|
|
outfile.width(16);
|
|
outfile << p.Y()/scale << ", ";
|
|
|
|
if(mesh[PointIndex(i)].Type() != INNERPOINT)
|
|
{
|
|
BCsinpoint.DeleteAll();
|
|
for (j = 1; j <= nse; j++)
|
|
{
|
|
for (k = 1; k <= 2; k++)
|
|
{
|
|
if(mesh.SurfaceElement(j).PNum(k)==i)
|
|
{
|
|
int BC=mesh.GetFaceDescriptor(mesh.SurfaceElement(j).GetIndex()).BCProperty();
|
|
int nbcsp=BCsinpoint.Size();
|
|
int found = 0;
|
|
for (l = 1; l <= nbcsp; l++)
|
|
if(BC == BCsinpoint.Get(l)) found = 1;
|
|
if( ! found ) BCsinpoint.Append(BC);
|
|
}
|
|
}
|
|
}
|
|
int nbcsp = BCsinpoint.Size();
|
|
outfile << "[" << nbcsp << "] ";
|
|
for (j = 1; j <= nbcsp; j++)
|
|
outfile << BCsinpoint.Get(j) << " ";
|
|
outfile << "\n";
|
|
}
|
|
else outfile << "[0]\n";
|
|
|
|
}
|
|
|
|
outfile << "\n"
|
|
" Element types and connectivity\n"
|
|
" the columns contain:\n"
|
|
" - element number\n"
|
|
" - element type\n"
|
|
" - subdomain number\n"
|
|
" - the global node numbers of the nodes in the element.\n"
|
|
"#\n";
|
|
|
|
for (i = 1; i <= nse; i++)
|
|
{
|
|
const Element2d & el = mesh.SurfaceElement(i);
|
|
outfile.width(12);
|
|
if(eldummy.GetNP()==3)
|
|
outfile << i << " ElmT3n2D ";
|
|
else
|
|
outfile << i << " ElmT6n2D ";
|
|
outfile.width(12);
|
|
outfile << el.GetIndex() << " ";
|
|
outfile.width(16);
|
|
outfile << el.PNum(1);
|
|
outfile.width(16);
|
|
outfile << el.PNum(2);
|
|
outfile.width(16);
|
|
outfile << el.PNum(3);
|
|
if(eldummy.GetNP()==6)
|
|
{
|
|
outfile.width(16);
|
|
outfile << el.PNum(6);
|
|
outfile.width(16);
|
|
outfile << el.PNum(4);
|
|
outfile.width(16);
|
|
outfile << el.PNum(5);
|
|
}
|
|
outfile << "\n";
|
|
}
|
|
}
|
|
}
|
|
}
|