Commit 1b8034a6 authored by prudhomm's avatar prudhomm
Browse files

added netgen read/writer(solution) by Mauro

parent 3691cffb
......@@ -15,6 +15,7 @@ GF :: G. Fourestey
JFG :: J.F. Gerbeau
LF :: L. Formaggia
MF :: M. Fernandez
MM :: M. Meneghin
MP :: M. Prosi
SD :: S. Deparis
VM :: V. Martin
......
......@@ -31,6 +31,9 @@ New in 0.5.0:
* Factorized and merged FSI::operFS (GF)
* Added NavierStokesWithFlux solver to impose fluxes (CV, CP)
[lifefilters]
* Added netgen reader/writer(solution) (MM)
[testsuite]
* Added test for 3D/0D coupling (AM)
* Moved files of test_darcy in lifesolver (VM, CP)
......@@ -186,4 +189,4 @@ New in 0.1.0pre3:
* fixed some bugs and miscompilation in a few places
-- Christophe Prud'homme <prudhomm@debian.org>, Thu Oct 28 19:18:06 2004
-- Christophe Prud'homme <prudhomm@debian.org>, Mon Nov 8 14:55:09 2004
......@@ -633,5 +633,368 @@ readINRIAMeshFileHead( std::ifstream & mystream,
UInt & numVolumes,
ReferenceShapes & shape );
//meneghin: added support for reading mesh3d meshes generated by netgen
template<typename GeoShape, typename MC>
bool
readNetgenMesh(RegionMesh3D<GeoShape,MC> & mesh,
const std::string & filename,
EntityFlag regionFlag);
/*
2004 Mauro Meneghin
added support for using meshes
generated by netgen
I referred a bit to:
LifeV/.../mesh/readMesh3D.h/.cc
netgen/libsrc/meshing/meshclass.cpp
#include "markers.h"
#include "regionMesh3D.h"
#include "util_string.h"
#include "mesh_util.h"
*/
template<typename GeoShape, typename MC>
bool
readNetgenMesh(RegionMesh3D<GeoShape,MC> & mesh,
const std::string & filename,
EntityFlag regionFlag)
{
std::string line;
UInt nVe(0), nBVe(0),nFa(0),nBFa(0),nPo(0),nBPo(0);
UInt nVo(0),nEd(0),nBEd(0);
UInt i;
UInt fake;
UInt surfnr,bcnr,domin,domout;
UInt surf1,surf2,matnr,np,p1,p2,p3,p4,t;
Real x,y,z;
std::vector<bool> bpoints;
// std::vector<bool> bedges;
BareItemsHandler<BareEdge> bihBedges;
std::vector<EntityFlag> bcnsurf,bcnpoints;
UInt flag;
typename MC::PointMarker PMarker;
typename MC::EdgeMarker EMarker;
std::ifstream hstream(filename.c_str()); //.c_str());
if (hstream.fail()) {
std::cerr<<"Error in readNetgenMesh: File not found or locked"<<std::endl;
abort();
}
std::cout<< "Reading netgen mesh file: "<<filename<<std::endl;
hstream>>line;
if(line!="mesh3d"){
std::cerr<<"Error in readNetgenMesh: mesh file is not in mesh3d format (netgen)"\
<<std::endl;
abort();
}
/* I assume as I tested that faces stored are only boundary faces
and edges stored are only a part of boundary ones with
inside a face with a common marker, instead points can be
inside the domain too, but this format doesn't say me
which, so I'll find them from myself
*/
flag=1|2|4|8;
while(!hstream.eof())
{
hstream>>line;
if(line=="points" && flag&1){
hstream>>nVe;
bcnpoints.reserve(nVe+1);
bpoints.reserve(nVe+1);
for(i=0;i<nVe+1;i++){
bpoints[i]=false;
bcnpoints[i]=NULLFLAG;
}
flag&=~1;
break;
}
}
hstream.seekg(0,std::ios_base::beg);
while(!hstream.eof())
{
hstream>>line;
if(line=="surfaceelementsgi" && flag&8){
hstream>>nBFa;
bcnsurf.reserve(nBFa+1);
bcnsurf[0]=NULLFLAG;
for(i=0;i<nBFa;i++){
hstream>>surfnr>>bcnr>>domin>>domout>>np>>p1>>p2>>p3\
>>t>>t>>t;
bcnsurf[i+1]=bcnr;
if(np!=3){
std::cerr<<"Error in readNetgenMesh: " \
"only triangular surfaces supported"<<std::endl;
abort();
}
//assume p1!=p2!=p3!=p1
nBVe+=bpoints[p1]?0:1;
nBVe+=bpoints[p2]?0:1;
nBVe+=bpoints[p3]?0:1;
bpoints[p1]=bpoints[p2]=bpoints[p3]=true;
/* here I set the boundart points marker
note: this works only with my patch (meneghin)
of strongerFlag */
bcnpoints[p1]=PMarker.setStrongerMarker(bcnpoints[p1],bcnr);
bcnpoints[p2]=PMarker.setStrongerMarker(bcnpoints[p2],bcnr);
bcnpoints[p3]=PMarker.setStrongerMarker(bcnpoints[p3],bcnr);
/* now I have the surface and points, so I can calculate
the num of edges on the boundary, useful in case this is
a quadratic tetra, I calculate the bcnr too using
BareItemHandler<BareEdge> to store results
*/
//I've found this silly but easy way
BareEdge bed=setBareEdge(p1,p2);
bihBedges.addIfNotThere(bed,(ID)NULLFLAG);
bihBedges[bed]=(ID)EMarker.setStrongerMarker(
bihBedges[bed],
bcnr);
bed=setBareEdge(p2,p3);
bihBedges.addIfNotThere(bed,(ID)NULLFLAG);
bihBedges[bed]=(ID)EMarker.setStrongerMarker(
bihBedges[bed],
bcnr);
bed=setBareEdge(p3,p1);
bihBedges.addIfNotThere(bed,(ID)NULLFLAG);
bihBedges[bed]=(ID)EMarker.setStrongerMarker(
bihBedges[bed],
bcnr);
}
flag&=~8;
break;
}
}
nBEd=bihBedges.howMany();
/* I assume as I tested that edges stored are only a part of
boundary ones, so really they are meaningless to me */
hstream.seekg(0,std::ios_base::beg);
while(!hstream.eof())
{
hstream>>line;
if(line=="edgesegmentsgi2" && flag&2){
hstream>>fake; //this is not really the nBEd
for(i=0;i<fake;i++){
hstream>>surf1>>surf2>>p1>>p2>>t>>t>>t>>t>>t>>t>>t>>t;
}
flag&=~2;
}
else if(line=="volumeelements" && flag&4){
hstream>>nVo;
for(i=0;i<nVo;i++){
hstream>>matnr>>np>>p1>>p2>>p3>>p4;
if(np!=4){
std::cerr<<"Error in readNetgenMesh: " \
"only tetrahedra elements supported"<<std::endl;
abort();
}
}
flag&=~4;
}
if(!flag)break;
}
if(flag!=0){
std::cerr<<"Error in readNetgenMesh: " \
"the mesh file does not have all the required sections" \
<<std::endl;
abort();
}
// Euler formulas
nFa=2*nVo+(nBFa/2);
nEd=nVo+nVe+(3*nBFa-2*nBVe)/4;
// Be a little verbose
if (GeoShape::numPoints > 4 ){
std::cout << "Quadratic Tetra Mesh (from Linear geometry)" <<std::endl;
nPo=nVe+nEd;
nBPo=nBVe+nBEd; // I calculated the real nBEd before
} else {
std::cout << "Linear Tetra Mesh" <<std::endl;
nPo=nVe;
nBPo=nBVe;
}
//points can be only vertices or on edges too
std::cout<< "#Vertices= "<<nVe;
std::cout<< " #BVertices= "<<nBVe<<std::endl;
std::cout<< "#Faces= "<<nFa;
std::cout<< " #Boundary Faces= "<<nBFa<<std::endl;
std::cout<< "#Edges= "<<nEd;
std::cout<< " #Boundary Edges= "<<nBEd<<std::endl;
std::cout<< "#Points= "<<nPo;
std::cout<< " #Boundary Points= "<<nBPo<<std::endl;
std::cout<< "#Volumes= "<<nVo<<std::endl;
// Set all basic data structure
// I store all Points
mesh.setMaxNumPoints(nPo,true);
mesh.setNumBPoints(nBPo);
mesh.numVertices()=nVe;
mesh.numBVertices()=nBVe;
// Only Boundary Edges (in a next version I will allow for different choices)
mesh.setMaxNumEdges(nBEd);
mesh.numEdges()=nEd; // Here the REAL number of edges (all of them)
mesh.setNumBEdges(nBEd); /////////????????????
// Only Boundary Faces
mesh.setMaxNumFaces(nBFa);
mesh.numFaces()=nFa; // Here the REAL number of edges (all of them)
mesh.setNumBFaces(nBFa);
mesh.setMaxNumVolumes(nVo,true);
mesh.setMarker(regionFlag); // Mark the region ????????what if more then one<<<<<<<
typename RegionMesh3D<GeoShape,MC>::PointType * pp=0;
typename RegionMesh3D<GeoShape,MC>::EdgeType * pe=0;
typename RegionMesh3D<GeoShape,MC>::FaceType * pf=0;
typename RegionMesh3D<GeoShape,MC>::VolumeType * pv=0;
// addPoint()/Face()/Edge() returns a reference to the last stored point
// I use that information to set all point info, by using a pointer.
hstream.seekg(0,std::ios_base::beg);
flag=1|2|4|8;
while(!hstream.eof())
{
hstream>>line;
if(line=="points" && flag&1){
hstream>>nVe;
for(i=0;i<nVe;i++){
hstream>>x>>y>>z;
pp=&mesh.addPoint(bpoints[i+1]); //true if boundary point
pp->setMarker(bcnpoints[i+1]);
pp->x()=x;
pp->y()=y;
pp->z()=z;
}
flag&=~1;
break;
}
}
hstream.seekg(0,std::ios_base::beg);
while(!hstream.eof())
{
hstream>>line;
if(line=="edgesegmentsgi2" && flag&2){
hstream>>fake; //this is not really the nBEd
for(i=0;i<fake;i++){
hstream>>surf1>>surf2>>p1>>p2>>t>>t>>t>>t>>t>>t>>t>>t;
}
/* here I set the real boundary edges that I stored
in bihBedges
*/
BareItemsHandler<BareEdge>::const_iterator bedge=bihBedges.begin();
for(i=0;i<nBEd;i++){
pe = &mesh.addEdge( true ); // Only boundary edges.
pe->setMarker( EntityFlag(bedge->second));
p1=bedge->first.first;
p2=bedge->first.second;
pe->setPoint( 1, mesh.point( p1 ) ); // set edge conn.
pe->setPoint( 2, mesh.point( p2 ) ); // set edge conn.
bedge++;
}
flag&=~2;
}
else if(line=="volumeelements" && flag&4){
hstream>>nVo;
for(i=0;i<nVo;i++){
hstream>>matnr>>np>>p1>>p2>>p3>>p4;
pv=&mesh.addVolume();
pv->id()=i+1;
pv->setPoint(1, mesh.point(p1) );
pv->setPoint(2, mesh.point(p2) );
pv->setPoint(3, mesh.point(p3) );
pv->setPoint(4, mesh.point(p4) );
}
flag&=~4;
}
else if(line=="surfaceelementsgi" && flag&8){
hstream>>nFa;
for(i=0;i<nFa;i++){
hstream>>surfnr>>bcnr>>domin>>domout>>np>>p1>>p2>>p3\
>>t>>t>>t;
pf=&mesh.addFace(true); // Only boundary faces
pf->setMarker(EntityFlag(bcnr));
pf->setPoint(1,mesh.point(p1)); // set face conn.
pf->setPoint(2,mesh.point(p2)); // set face conn.
pf->setPoint(3,mesh.point(p3)); // set face conn.
}
flag&=~8;
}
}
// This part is to build a P2 mesh from a P1 geometry
if (GeoShape::numPoints > 4 )p1top2(mesh);
hstream.close();
// Test mesh
Switch sw;
///// CORRECTION JFG
//if (mesh.check(1, true,true))done=0;
if(!checkMesh3D(mesh, sw, true,true,std::cout,std::cout,std::cout)) abort(); // CORRECTION JFG
Real vols[3];
getVolumeFromFaces(mesh, vols,std::cout);
std::cout<< " VOLUME ENCLOSED BY THE MESH COMPUTED BY INTEGRATION ON"<<
" BOUNDARY FACES"<<std::endl;
std::cout << "INT(X) INT(Y) INT(Z) <- they should be equal and equal to"<<std::endl<<
" the voulume enclosed by the mesh "<<std::endl;
std::cout<<vols[0]<<" "<<vols[1]<<" "<<vols[2]<<std::endl;
std::cout<< " BOUNDARY FACES ARE DEFINING A CLOSED SURFACE IF "<<testClosedDomain(mesh,std::cout)<< std::endl<<
" IS (ALMOST) ZERO"<<std::endl;
return true;
}
/*template <typename TheMesh,typename TheDof,typename VectorType>
void saveNetgenSolution(std::string filename,VectorType U,
TODO!!!!!!! order U output with same order of mesh points
*/
#include <iostream>
#include <fstream>
/*
/usr/local/src/ng431/libsrc/interface/importsolution.cpp
*/
template <typename VectorType>
void saveNetgenSolution(std::string filename,VectorType U)
{
std::ofstream of(filename.c_str());
std::string fctname=std::string("f");
of<<"solution "<<fctname<<" -size="<<U.size()
<<" -components=1 -type=nodal"<<std::endl;
for(UInt i=0;i<U.size();i++)
of<<U(i)<<std::endl;
of.close();
}
}
#endif
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment