netgen/libsrc/interface/writeelmer.cpp

209 lines
5.1 KiB
C++
Raw Normal View History

2009-01-13 04:40:13 +05:00
//
// Write Elmer file
//
//
#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <csg.hpp>
#include <meshing.hpp>
#include <sys/stat.h>
2016-07-10 21:07:36 +05:00
#include "writeuser.hpp"
2009-01-13 04:40:13 +05:00
namespace netgen
{
extern MeshingParameters mparam;
2009-01-13 04:40:13 +05:00
void WriteElmerFormat (const Mesh &mesh,
2014-08-30 06:15:59 +06:00
const string &filename)
2009-01-13 04:40:13 +05:00
{
cout << "write elmer mesh files" << endl;
std::map<ELEMENT_TYPE, int> tmap;
tmap[TRIG] = 303;
tmap[TRIG6] = 306;
tmap[QUAD] = 404;
tmap[QUAD8] = 408;
tmap[TET] = 504;
tmap[TET10] = 510;
tmap[PYRAMID] = 605;
tmap[PYRAMID13] = 613;
tmap[PRISM] = 706;
tmap[PRISM15] = 715;
tmap[HEX] = 808;
tmap[HEX20] = 820;
std::map<int, Array<int,int>> pmap;
pmap[TRIG] = {1,2,3};
pmap[TRIG6] = {1,2,3, 6,4,5};
pmap[QUAD] = {1,2,3,4};
pmap[QUAD8] = {1,2,3,4, 5,8,6,7};
pmap[TET] = {1,2,3,4};
pmap[TET10] = {1,2,3,4, 5,8,6,7,9,10};
pmap[PYRAMID]={1,2,3,4,5};
pmap[PYRAMID13]= {1,2,3,4,5,6,7,8,9,10,11,12,13};
pmap[PRISM] = {1,2,3,4,5,6};
pmap[PRISM15] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
pmap[HEX] = {1,2,3,4,5,6,7,8};
pmap[HEX20] = {1,2,3,4,5,8,6,7,8, 9,12,10,11, 17,20,19,18, 13,16,14,15};
2009-01-13 04:40:13 +05:00
int np = mesh.GetNP();
int ne = mesh.GetNE();
int nse = mesh.GetNSE();
int i, j;
char str[200];
int inverttets = mparam.inverttets;
int invertsurf = mparam.inverttrigs;
#ifdef WIN32
char a[256];
sprintf( a, "mkdir %s", filename.c_str() );
system( a );
#else
2009-04-20 03:15:26 +06:00
// int rc =
mkdir(filename.c_str(), S_IRWXU|S_IRWXG);
2009-01-13 04:40:13 +05:00
#endif
sprintf( str, "%s/mesh.header", filename.c_str() );
ofstream outfile_h(str);
sprintf( str, "%s/mesh.nodes", filename.c_str() );
ofstream outfile_n(str);
sprintf( str, "%s/mesh.elements", filename.c_str() );
ofstream outfile_e(str);
sprintf( str, "%s/mesh.boundary", filename.c_str() );
ofstream outfile_b(str);
sprintf( str, "%s/mesh.names", filename.c_str() );
ofstream outfile_names(str);
for( auto codim : IntRange(0, mesh.GetDimension()-1) )
{
auto & names = const_cast<Mesh&>(mesh).GetRegionNamesCD(codim);
for (auto i0 : Range(names) )
{
if(names[i0] == nullptr)
continue;
string name = *names[i0];
if(name == "" || name == "default")
continue;
outfile_names << "$" << name << "=" << i0+1 << "\n";
}
}
auto get3FacePoints = [](const Element2d & el)
{
INDEX_3 i3;
INDEX_4 i4;
auto eltype = el.GetType();
switch (eltype)
{
case TRIG:
case TRIG6:
i3 = {el[0], el[1], el[2]};
i3.Sort();
break;
case QUAD:
case QUAD8:
i4 = {el[0], el[1], el[2], el[3]};
i4.Sort();
i3 = {i4[0], i4[1], i4[2]};
break;
default:
throw Exception("Got invalid type (no face)");
}
return i3;
};
2009-01-13 04:40:13 +05:00
// fill hashtable
// use lowest three point numbers of lowest-order face to index faces
2009-01-13 04:40:13 +05:00
INDEX_3_HASHTABLE<int> face2volelement(ne);
for (int i = 1; i <= ne; i++)
2009-01-13 04:40:13 +05:00
{
const Element & el = mesh.VolumeElement(i);
// getface not working for second order elements -> reconstruct linear element here
Element linear_el = el;
linear_el.SetNP(el.GetNV()); // GetNV returns 8 for HEX20 for instance
for (auto j : Range(1,el.GetNFaces()+1))
2009-01-13 04:40:13 +05:00
{
Element2d face;
linear_el.GetFace(j, face);
face2volelement.Set (get3FacePoints(face), i);
cout << "set " << get3FacePoints(face) << "\tto " << i << endl;
2009-01-13 04:40:13 +05:00
}
}
// outfile.precision(6);
// outfile.setf (ios::fixed, ios::floatfield);
// outfile.setf (ios::showpoint);
std::map<ELEMENT_TYPE, size_t> elcount;
2009-01-13 04:40:13 +05:00
for (i = 1; i <= np; i++)
{
const Point3d & p = mesh.Point(i);
outfile_n << i << " -1 ";
outfile_n << p.X() << " ";
outfile_n << p.Y() << " ";
outfile_n << p.Z() << "\n";
}
for (i = 1; i <= ne; i++)
{
Element el = mesh.VolumeElement(i);
if (inverttets) el.Invert();
auto eltype = el.GetType();
elcount[eltype]++;
outfile_e << i << " " << el.GetIndex() << " " << tmap[eltype] << " ";
auto & map = pmap[eltype];
2009-01-13 04:40:13 +05:00
for (j = 1; j <= el.GetNP(); j++)
{
outfile_e << " ";
outfile_e << el.PNum(map[j-1]);
2009-01-13 04:40:13 +05:00
}
outfile_e << "\n";
}
for (i = 1; i <= nse; i++)
{
Element2d el = mesh.SurfaceElement(i);
if (invertsurf) el.Invert();
auto eltype = el.GetType();
elcount[eltype]++;
int elind = face2volelement.Get(get3FacePoints(el));
cout << "get " << get3FacePoints(el) << "\t " << elind << endl;
outfile_b << i << " " << mesh.GetFaceDescriptor(el.GetIndex()).BCProperty() <<
" " << elind << " 0 " << tmap[eltype] << " ";
auto & map = pmap[el.GetType()];
2009-01-13 04:40:13 +05:00
for (j = 1; j <= el.GetNP(); j++)
{
outfile_b << " ";
outfile_b << el.PNum(map[j-1]);
2009-01-13 04:40:13 +05:00
}
outfile_b << "\n";
}
outfile_h << np << " " << ne << " " << nse << "\n";
outfile_h << "2" << "\n";
for( auto & [eltype,count] : elcount )
outfile_h << tmap[eltype] << " " << count << "\n";
2009-01-13 04:40:13 +05:00
}
}