ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-15 10:28:34 +05:00
Code Issues Packages Projects Releases Wiki Activity
f1ec19b35b
netgen/libsrc/interface/writeOpenFOAM15x.cpp
Joachim Schoeberl e1f7935f0b fixing warnings
2023-08-05 12:01:01 +02:00

762 lines
26 KiB
C++
Raw Blame History

/*! \file writeOpenFOAM15x.cpp
* \brief Export Netgen Mesh in the OpenFOAM 1.5+ File format
* \author Philippose Rajan
* \date 25 October 2009
*
* This function extends the export capabilities of
* Netgen to include the OpenFOAM 1.5+ File Format.
*
* The OpenFOAM 1.5+ mesh format consists of a set of 5 files
* which together define the mesh points, faces, cells and
* boundary conditions.
*
* The files are:
* 1. points -> A list of the point co-ordinates
* 2. faces -> A list of the faces with format <n>(pnt_ind1 pnt_ind2 .... pnt_ind<n>)
* 3. owner -> The owner cell of each face
* 4. neighbour -> The neighbour cell of each face
* 5. boundary -> The set of boundaries with name, start face, and num. of faces
*
* For a detailed description of the format, refer to the following link:
* http://openfoamwiki.net/index.php/Write_OpenFOAM_meshes
*
*/
#include <mystdlib.h>
#include <myadt.hpp>
#include <linalg.hpp>
#include <csg.hpp>
#include <meshing.hpp>
#include <sys/stat.h>
#include "../general/gzstream.h"
namespace netgen
{
#include "writeuser.hpp"
extern MeshingParameters mparam;
// Global arrays used to maintain the owner, neighbour and face lists
// so that they are accessible across functions
static NgArray<int> owner_facelist;
static NgArray<int> owner_celllist;
static NgArray<int> neighbour_celllist;
static NgArray<int> surfelem_bclist;
static NgArray<INDEX_2> surfelem_lists;
static void WriteOpenFOAM15xBanner(ostream * outfile)
{
static char FOAMversion[4] = "1.5";
static char spaces[40];
memset(spaces, ' ', 40);
spaces[38 - strlen(FOAMversion)] = '\0';
*outfile <<
"/*--------------------------------*- C++ -*----------------------------------*\\\n";
*outfile <<
"| ========= | |\n"
"| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\n"
"| \\\\ / O peration | Version: " << FOAMversion << spaces << "|\n"
"| \\\\ / A nd | Web: http://www.OpenFOAM.org |\n"
"| \\\\/ M anipulation | |\n"
"\\*---------------------------------------------------------------------------*/\n";
}
static void WriteOpenFOAM15xDividerStart(ostream * outfile)
{
*outfile <<
"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n";
}
static void WriteOpenFOAM15xDividerEnd(ostream * outfile)
{
*outfile <<
"// ************************************************************************* //\n";
}
static bool BuildOwnerNeighbourLists (const Mesh & mesh)
{
// Clear all the arrays
owner_facelist.DeleteAll();
owner_celllist.DeleteAll();
neighbour_celllist.DeleteAll();
surfelem_bclist.DeleteAll();
surfelem_lists.DeleteAll();
const MeshTopology& meshtopo = mesh.GetTopology();
// Update the mesh topology structures
const_cast<MeshTopology&> (meshtopo).SetBuildEdges(true);
const_cast<MeshTopology&> (meshtopo).SetBuildFaces(true);
const_cast<MeshTopology&> (meshtopo).Update();
// Extract important mesh metrics
int ne = mesh.GetNE();
int nse = mesh.GetNSE();
int totfaces = meshtopo.GetNFaces();
// Preset the size of the arrays to speed up future operations
// Number of internal faces = total faces - num. of surface faces
owner_facelist.SetSize(totfaces - nse);
owner_celllist.SetSize(totfaces - nse);
neighbour_celllist.SetSize(totfaces - nse);
surfelem_bclist.SetSize(nse);
surfelem_lists.SetSize(nse);
// Initialise arrays to zero if required
neighbour_celllist = 0;
// NgArray used to keep track of Faces which have already been
// processed and added to the Owner list... In addition, also the
// location where the face appears in the Owner list is also stored
// to speed up creation of the Neighbour list
NgArray<int> ownerfaces(totfaces);
ownerfaces = 0;
// NgArray to hold the set of local faces of each volume element
// while running through the set of volume elements
// NOTE: The size is set automatically by the Netgen topology function
NgArray<int> locfaces;
// Secondary indices used to independently advance the owner
// and boundary condition arrays within the main loop
int owner_ind = 1;
int bc_ind = 1;
// Loop through all the volume elements
for(int elind = 1; elind <= ne; elind++)
{
// Extract the current volume element
// const Element & el = mesh.VolumeElement(elind);
// Get the face numbers of the faces of the current volume element
// The values returned are given a sign depending on the orientation
// of the faces. This is used while writing the faces file, to
// determine whether or not to invert the face triangle before writing
// it to file
meshtopo.GetElementFaces(elind,locfaces,true);
// Loop through the faces
for(int i = 1; i <= locfaces.Size(); i++)
{
// The absolute value of a face number (because the faces
// returned by the GetElementFaces function prepend it
// with a sign depending on the face orientation)
int absfacenr = abs(locfaces.Elem(i));
// If the face already exists in the owner list, add
// the current cell into the neighbour list, in the
// same location where the face appears in the owner list
int owner_face = ownerfaces.Elem(absfacenr);
if(owner_face)
{
neighbour_celllist.Elem(owner_face) = elind;
// From this point on, the code within this "if" block
// basically sorts the order of the Neighbour cells (along
// with the faces list) in ascending order.
// The approach used is..... to traverse the owner and neighbour cell lists
// up and down, and sort the neighbour cells of a given owner cell
// as the list evolves.
// NOTE: A value of "zero" in the neighbour list implies that
// the neighbour has not been found yet, so the "zero" locations need
// to be skipped while sorting in ascending order
int curr_owner = owner_celllist.Elem(owner_face);
int peek_loc = owner_face - 1;
int new_loc = owner_face;
// Traversing upwards in the list
while((owner_celllist.Elem(peek_loc) == curr_owner) && (peek_loc >= 1))
{
if((neighbour_celllist.Elem(peek_loc) != 0)
&& (neighbour_celllist.Elem(new_loc) < neighbour_celllist.Elem(peek_loc)))
{
Swap(neighbour_celllist.Elem(new_loc),neighbour_celllist.Elem(peek_loc));
Swap(owner_facelist.Elem(new_loc),owner_facelist.Elem(peek_loc));
new_loc = peek_loc;
}
peek_loc--;
}
peek_loc = owner_face + 1;
// Traversing downwards in the list
while((owner_celllist.Elem(peek_loc) == curr_owner) && (peek_loc <= owner_ind))
{
if((neighbour_celllist.Elem(peek_loc) != 0)
&& (neighbour_celllist.Elem(new_loc) > neighbour_celllist.Elem(peek_loc)))
{
Swap(neighbour_celllist.Elem(new_loc),neighbour_celllist.Elem(peek_loc));
Swap(owner_facelist.Elem(new_loc),owner_facelist.Elem(peek_loc));
new_loc = peek_loc;
}
peek_loc++;
}
continue;
}
// Check if the face is a surface element (boundary face)
// if not, add the current volume element and the corresponding face into
// the owner list
int surfelem = meshtopo.GetFace2SurfaceElement(absfacenr);
if(!surfelem)
{
// If it is a new face which has not been listed before,
// add the current cell into the owner list, and save
// the index location to be used later by the neighbour list
owner_celllist.Elem(owner_ind) = elind;
owner_facelist.Elem(owner_ind) = locfaces.Elem(i);
// Update the array to indicate that the face is already processed
ownerfaces.Elem(absfacenr) = owner_ind;
owner_ind++;
}
// If the face is a boundary face, extract the boundary condition number of the
// face, and append that along with the face number and the current cell
// into the various surface elements lists
else
{
Element2d sel = mesh.SurfaceElement(surfelem);
surfelem_bclist.Elem(bc_ind) = mesh.GetFaceDescriptor(sel.GetIndex()).BCProperty();
surfelem_lists.Elem(bc_ind) = INDEX_2(locfaces.Elem(i),elind);
bc_ind++;
}
}
}
// This correction is required in cases where the mesh has been "uniform refined".... for
// some reason, the number of faces reported by Netgen is higher than the actual number
// of faces in the mesh
owner_facelist.SetSize(owner_ind-1);
owner_celllist.SetSize(owner_ind-1);
neighbour_celllist.SetSize(owner_ind-1);
// Sort the list of surface elements in ascending order of boundary condition number
// also sort the cell list in the same manner
QuickSort(surfelem_bclist,surfelem_lists);
/*
// Debugging output to a file
ofstream dbg("OpenFOAMDebug.log");
dbg << " ------- Boundary List -------- \n";
for(int i = 1; i <= surfelem_bclist.Size(); i++)
{
dbg << "bc = " << surfelem_bclist.Elem(i)
<< " : face = " << surfelem_lists.Elem(i).I1()
<< " : cell = " << surfelem_lists.Elem(i).I2() << "\n";
}
dbg << "\n ------- Owner / Face / Neighbour List ------- \n";
for(int i = 1; i <= owner_celllist.Size(); i++)
{
dbg << "Ind:" << i << " :: ("
<< owner_celllist.Elem(i) << " "
<< owner_facelist.Elem(i) << " "
<< neighbour_celllist.Elem(i) << ")\n";
}
dbg.close();
*/
return(false);
}
static void WriteNeighbourFile (ostream * outfile)
{
// Write the OpenFOAM standard banner and dividers, etc...
WriteOpenFOAM15xBanner(outfile);
*outfile << "FoamFile \n"
<< "{ \n"
<< " version 2.0; \n"
<< " format ascii; \n"
<< " class labelList; \n"
<< " note \"Mesh generated and converted using NETGEN-" << PACKAGE_VERSION << "\"; \n"
<< " location \"constant\\polyMesh\"; \n"
<< " object neighbour; \n"
<< "} \n";
WriteOpenFOAM15xDividerStart(outfile);
*outfile << "\n\n";
int nneighbours = neighbour_celllist.Size();
*outfile << nneighbours << "\n";
*outfile << "(\n";
// Write the neighbour cells to file
for(int i = 1; i <= neighbour_celllist.Size(); i++)
{
*outfile << neighbour_celllist.Elem(i) - 1 << "\n";
}
*outfile << ")\n\n";
WriteOpenFOAM15xDividerEnd(outfile);
}
static void WriteOwnerFile (ostream * outfile)
{
// Write the OpenFOAM standard banner and dividers, etc...
WriteOpenFOAM15xBanner(outfile);
*outfile << "FoamFile \n"
<< "{ \n"
<< " version 2.0; \n"
<< " format ascii; \n"
<< " class labelList; \n"
<< " note \"Mesh generated and converted using NETGEN-" << PACKAGE_VERSION << "\"; \n"
<< " location \"constant\\polyMesh\"; \n"
<< " object owner; \n"
<< "} \n";
WriteOpenFOAM15xDividerStart(outfile);
*outfile << "\n\n";
int nowners = owner_celllist.Size() + surfelem_lists.Size();
*outfile << nowners << "\n";
*outfile << "(\n";
// Write the owners of the internal cells to file
for(int i = 1; i <= owner_celllist.Size(); i++)
{
*outfile << owner_celllist.Elem(i) - 1 << "\n";
}
// Write the owners of the boundary cells to file
// (Written in order of ascending boundary condition numbers)
for(int i = 1; i <= surfelem_lists.Size(); i++)
{
*outfile << surfelem_lists.Elem(i).I2() - 1 << "\n";
}
*outfile << ")\n\n";
WriteOpenFOAM15xDividerEnd(outfile);
}
static void WriteFacesFile (ostream * outfile, const Mesh & mesh)
{
const MeshTopology& meshtopo = mesh.GetTopology();
// Write the OpenFOAM standard banner and dividers, etc...
WriteOpenFOAM15xBanner(outfile);
*outfile << "FoamFile \n"
<< "{ \n"
<< " version 2.0; \n"
<< " format ascii; \n"
<< " class faceList; \n"
<< " note \"Mesh generated and converted using NETGEN-" << PACKAGE_VERSION << "\"; \n"
<< " location \"constant\\polyMesh\"; \n"
<< " object faces; \n"
<< "} \n";
WriteOpenFOAM15xDividerStart(outfile);
*outfile << "\n\n";
int nfaces = owner_facelist.Size() + surfelem_lists.Size();
*outfile << nfaces << "\n";
*outfile << "(\n";
// NgArray to hold the indices of the points of each face to
// flip if required
NgArray<int> facepnts;
// Write the faces in the order specified in the owners lists of the
// internal cells and the boundary cells
for(int i = 1; i <= owner_facelist.Size(); i++)
{
int face_w_orientation = owner_facelist.Elem(i);
int facenr = abs(face_w_orientation);
meshtopo.GetFaceVertices(facenr,facepnts);
// Get the orientation of the face, and invert it if required
// Since the faces already have the orientation "embedded" into
// them by means of the prepended sign, only this needs to be
// checked for...
if(face_w_orientation > 0)
{
int tmppnts = 0;
if(facepnts.Size() == 4)
{
tmppnts = facepnts.Elem(1);
facepnts.Elem(1) = facepnts.Elem(2);
facepnts.Elem(2) = tmppnts;
tmppnts = facepnts.Elem(3);
facepnts.Elem(3) = facepnts.Elem(4);
facepnts.Elem(4) = tmppnts;
}
else if(facepnts.Size() == 3)
{
tmppnts = facepnts.Elem(1);
facepnts.Elem(1) = facepnts.Elem(3);
facepnts.Elem(3) = tmppnts;
}
}
*outfile << facepnts.Size();
*outfile << "(";
for(int j = 1; j <= facepnts.Size(); j++)
{
*outfile << facepnts.Elem(j)-1;
if(j != facepnts.Size()) *outfile << " ";
}
*outfile << ")\n";
}
// Now append the faces of the surface elements (written in
// ascending order of boundary condition number) also into
// the faces file
for(int i = 1; i <= surfelem_lists.Size(); i++)
{
int face_w_orientation = surfelem_lists.Elem(i).I1();
int facenr = abs(face_w_orientation);
meshtopo.GetFaceVertices(facenr,facepnts);
// Get the orientation of the face, and invert it if required
if(face_w_orientation > 0)
{
int tmppnts = 0;
if(facepnts.Size() == 4)
{
tmppnts = facepnts.Elem(1);
facepnts.Elem(1) = facepnts.Elem(2);
facepnts.Elem(2) = tmppnts;
tmppnts = facepnts.Elem(3);
facepnts.Elem(3) = facepnts.Elem(4);
facepnts.Elem(4) = tmppnts;
}
else if(facepnts.Size() == 3)
{
tmppnts = facepnts.Elem(1);
facepnts.Elem(1) = facepnts.Elem(3);
facepnts.Elem(3) = tmppnts;
}
}
*outfile << facepnts.Size();
*outfile << "(";
for(int j = 1; j <= facepnts.Size(); j++)
{
*outfile << facepnts.Elem(j)-1;
if(j != facepnts.Size()) *outfile << " ";
}
*outfile << ")\n";
}
*outfile << ")\n\n";
WriteOpenFOAM15xDividerEnd(outfile);
}
static void WritePointsFile (ostream * outfile, const Mesh & mesh)
{
int np = mesh.GetNP();
// Write the OpenFOAM standard banner and dividers, etc...
WriteOpenFOAM15xBanner(outfile);
*outfile << "FoamFile \n"
<< "{ \n"
<< " version 2.0; \n"
<< " format ascii; \n"
<< " class vectorField; \n"
<< " note \"Mesh generated and converted using NETGEN-" << PACKAGE_VERSION << "\"; \n"
<< " location \"constant\\polyMesh\"; \n"
<< " object points; \n"
<< "} \n";
WriteOpenFOAM15xDividerStart(outfile);
*outfile << "\n\n";
// Number of points in the following list
*outfile << np << "\n";
outfile->precision(6);
outfile->setf (ios::fixed, ios::floatfield);
outfile->setf (ios::showpoint);
// Coordinate list starts here
*outfile << "(\n";
for(int i = 1; i <= np; i++)
{
const Point3d & p = mesh.Point(i);
// Write coordinates to file
*outfile << "(";
*outfile << p.X() << " ";
*outfile << p.Y() << " ";
*outfile << p.Z();
*outfile << ")\n";
}
*outfile << ")\n\n";
WriteOpenFOAM15xDividerEnd(outfile);
}
static void WriteBoundaryFile (ostream * outfile)
{
// Write the OpenFOAM standard banner and dividers, etc...
WriteOpenFOAM15xBanner(outfile);
*outfile << "FoamFile \n"
<< "{ \n"
<< " version 2.0; \n"
<< " format ascii; \n"
<< " class polyBoundaryMesh; \n"
<< " note \"Mesh generated and converted using NETGEN-" << PACKAGE_VERSION << "\"; \n"
<< " location \"constant\\polyMesh\"; \n"
<< " object boundary; \n"
<< "} \n";
WriteOpenFOAM15xDividerStart(outfile);
*outfile << "\n";
NgArray<INDEX_3> bcarray;
int ind = 1;
// Since the boundary conditions are already sorted in ascending
// order, the last element will give the maximum number of possible
// boundary condition entries
int bcmax = surfelem_bclist.Elem(surfelem_bclist.Size());
bcarray.SetSize(bcmax+1);
bcarray.Elem(ind) = INDEX_3(surfelem_bclist.Elem(1),1,0);
for(int i = 2; i <= surfelem_bclist.Size(); i++)
{
if(surfelem_bclist.Elem(i) == bcarray.Elem(ind).I1())
{
bcarray.Elem(ind).I2() = bcarray.Elem(ind).I2()+1;
}
else
{
ind++;
bcarray.Elem(ind) = INDEX_3(surfelem_bclist.Elem(i),1,i-1);
}
}
bcarray.SetSize(ind);
*outfile << bcarray.Size() << "\n";
*outfile << "(\n";
int startface = 0;
for(int i = 1; i <= bcarray.Size(); i++)
{
startface = owner_celllist.Size() + bcarray.Elem(i).I3();
*outfile << " patch" << bcarray.Elem(i).I1() << "\n"
<< " {\n"
<< " type patch;\n"
<< " physicalType patch;\n"
<< " nFaces " << bcarray.Elem(i).I2() << ";\n"
<< " startFace " << startface << ";\n"
<< " }\n";
}
*outfile << ")\n\n";
WriteOpenFOAM15xDividerEnd(outfile);
}
void WriteOpenFOAM15xFormat (const Mesh & mesh, const filesystem::path & dirname, const bool compressed)
{
bool error = false;
// char casefiles[256];
// Make sure that the mesh data has been updated
const_cast<Mesh&> (mesh).Compress();
const_cast<Mesh&> (mesh).CalcSurfacesOfNode();
const_cast<Mesh&> (mesh).RebuildSurfaceElementLists();
const_cast<Mesh&> (mesh).BuildElementSearchTree();
int np = mesh.GetNP();
int nse = mesh.GetNSE();
int ne = mesh.GetNE();
cout << "Write OpenFOAM 1.5+ Mesh Files....\n";
// Abort if there are no points, surface elements or volume elements
if((np <= 0) || (ne <= 0) || (nse <= 0))
{
cout << "Export Error: Invalid mesh.... Aborting!\n";
return;
}
// OpenFOAM only supports linear meshes!
if(mparam.secondorder || mesh.GetCurvedElements().IsHighOrder())
{
cout << "Export Error: OpenFOAM 1.5+ does not support non-linear elements.... Aborting!\n";
return;
}
if(( (mesh.SurfaceElement(nse/2).GetType() != TRIG)
&& (mesh.SurfaceElement(nse/2).GetType() != QUAD) )
|| (mesh.VolumeElement(ne/2).GetType() == TET10)
|| (mesh.VolumeElement(ne/2).GetType() == PRISM12))
{
cout << "Export Error: OpenFOAM 1.5+ does not support non-linear elements.... Aborting!\n";
return;
}
cout << "Writing OpenFOAM 1.5+ Mesh files to case: " << dirname.string() << "\n";
// Create the case directory if it does not already exist
// NOTE: This needs to be improved for the Linux variant....!!!
auto mesh_dir = filesystem::path(dirname).append("constant").append("polyMesh");
filesystem::create_directories(mesh_dir);
// Open handles to the five required mesh files
// points
// faces
// owner
// neighbour
// boundary
auto get_name = [compressed, &mesh_dir]( string s ) {
auto p = filesystem::path(mesh_dir).append(s);
if(compressed)
p.concat(".gz");
return p;
};
auto outfile_pnts = make_unique<ofstream>(get_name("points"));
auto outfile_faces = make_unique<ofstream>(get_name("faces"));
auto outfile_own = make_unique<ofstream>(get_name("owner"));
auto outfile_nei = make_unique<ofstream>(get_name("neighbor"));
// Note... the boundary file is not compressed
auto outfile_bnd = make_unique<ofstream>(mesh_dir.append("boundary"));
ResetTime();
// Build the owner, neighbour, faces and boundary lists
// from the Netgen mesh
cout << "\nBuilding Owner, Neighbour and Face Lists: ";
error = BuildOwnerNeighbourLists(mesh);
cout << "Done! (Time Elapsed = " << GetTime() << " sec)\n";
// Write the "owner" file
if(outfile_own->good() && !error)
{
cout << "Writing the owner file: ";
WriteOwnerFile(outfile_own.get());
cout << "Done! (Time Elapsed = " << GetTime() << " sec)\n";
}
else
{
cout << "Export Error: Error creating file: owner.... Aborting\n";
error = true;
}
// Write the "neighbour" file
if(outfile_nei->good() && !error)
{
cout << "Writing the neighbour file: ";
WriteNeighbourFile(outfile_nei.get());
cout << "Done! (Time Elapsed = " << GetTime() << " sec)\n";
}
else
{
cout << "Export Error: Error creating file: neighbour.... Aborting\n";
error = true;
}
// Write the "faces" file
if(outfile_faces->good() && !error)
{
cout << "Writing the faces file: ";
WriteFacesFile(outfile_faces.get(), mesh);
cout << "Done! (Time Elapsed = " << GetTime() << " sec)\n";
}
else
{
cout << "Export Error: Error creating file: faces.... Aborting\n";
error = true;
}
// Write the "points" file
if(outfile_pnts->good() && !error)
{
cout << "Writing the points file: ";
WritePointsFile(outfile_pnts.get(),mesh);
cout << "Done! (Time Elapsed = " << GetTime() << " sec)\n";
}
else
{
cout << "Export Error: Error creating file: points.... Aborting\n";
error = true;
}
// Write the "boundary" file
if(outfile_bnd->good() && !error)
{
cout << "Writing the boundary file: ";
WriteBoundaryFile(outfile_bnd.get());
cout << "Done! (Time Elapsed = " << GetTime() << " sec)\n";
}
else
{
cout << "Export Error: Error creating file: boundary.... Aborting\n";
error = true;
}
if(!error)
{
cout << "OpenFOAM 1.5+ Export successfully completed (Time elapsed = " << GetTime() << " sec) !\n";
}
else
{
cout << "Error in OpenFOAM 1.5+ Export.... Aborted!\n";
}
}
}
Reference in New Issue View Git Blame Copy Permalink