-// Copyright (C) 2004-2008 CEA/DEN, EDF R&D
+// Copyright (C) 2004-2010 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
+
//=============================================================================
// File : GHS3DPlugin_GHS3D.cxx
// Created :
// Author : Edward AGAPOV, modified by Lioka RAZAFINDRAZAKA (CEA) 09/02/2007
// Project : SALOME
-// $Header$
//=============================================================================
//
-using namespace std;
-
#include "GHS3DPlugin_GHS3D.hxx"
#include "GHS3DPlugin_Hypothesis.hxx"
+
+#include <Basics_Utils.hxx>
+
#include "SMESH_Gen.hxx"
#include "SMESH_Mesh.hxx"
#include "SMESH_Comment.hxx"
#include "SMDS_FaceOfNodes.hxx"
#include "SMDS_VolumeOfNodes.hxx"
+#include <StdMeshers_QuadToTriaAdaptor.hxx>
+
#include <BRepAdaptor_Surface.hxx>
#include <BRepBndLib.hxx>
#include <BRepClass3d_SolidClassifier.hxx>
#include <OSD_File.hxx>
#include <Precision.hxx>
#include <Quantity_Parameter.hxx>
+#include <Standard_ProgramError.hxx>
#include <Standard_ErrorHandler.hxx>
#include <Standard_Failure.hxx>
#include <TopExp.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopoDS.hxx>
//#include <BRepClass_FaceClassifier.hxx>
-//#include <BRepGProp.hxx>
-//#include <GProp_GProps.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <BRepGProp.hxx>
+#include <GProp_GProps.hxx>
#include "utilities.h"
-#ifndef WIN32
+#ifdef WIN32
+#include <io.h>
+#else
#include <sys/sysinfo.h>
#endif
+using namespace std;
+
//#include <Standard_Stream.hxx>
#ifdef _DEBUG_
#define DUMP(txt) \
-// cout << txt
+// std::cout << txt
#else
#define DUMP(txt)
#endif
#define HOLE_ID -1
+typedef const list<const SMDS_MeshFace*> TTriaList;
+
+static void removeFile( const TCollection_AsciiString& fileName )
+{
+ try {
+ OSD_File( fileName ).Remove();
+ }
+ catch ( Standard_ProgramError ) {
+ MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
+ }
+}
+
//=============================================================================
/*!
*
aPnt = p;
break;
}
- aPnt += p;
+ aPnt += p / nbNode;
}
BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
static char* readMapIntLine(char* ptr, int tab[]) {
long int intVal;
- cout << endl;
+ std::cout << std::endl;
for ( int i=0; i<17; i++ ) {
intVal = strtol(ptr, &ptr, 10);
//purpose :
//=======================================================================
-static bool writeFaces (ofstream & theFile,
- SMESHDS_Mesh * theMesh,
- const map <int,int> & theSmdsToGhs3dIdMap)
+static bool writeFaces (ofstream & theFile,
+ SMESHDS_Mesh* theMesh,
+ const TopoDS_Shape& theShape,
+ vector<StdMeshers_QuadToTriaAdaptor>& theQuad2Trias,
+ const map <int,int> & theSmdsToGhs3dIdMap)
{
// record structure:
//
// Loop from 1 to NB_ELEMS
// NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
- int nbShape = countShape( theMesh, TopAbs_FACE );
-
- int *tabID; tabID = new int[nbShape];
- TopoDS_Shape *tabShape; tabShape = new TopoDS_Shape[nbShape];
TopoDS_Shape aShape;
SMESHDS_SubMesh* theSubMesh;
const SMDS_MeshElement* aFace;
const int dummyint = 0;
map<int,int>::const_iterator itOnMap;
SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
- int shapeID, nbNodes, aSmdsID;
- bool idFound;
+ int nbNodes, aSmdsID;
- cout << " " << theMesh->NbFaces() << " shapes of 2D dimension" << endl;
- cout << endl;
+ // count triangles bound to geometry
+ int nbTriangles = 0;
- theFile << space << theMesh->NbFaces() << space << dummyint << endl;
+ TopTools_IndexedMapOfShape facesMap;
+ TopExp::MapShapes( theShape, TopAbs_FACE, facesMap );
- TopExp_Explorer expface( theMesh->ShapeToMesh(), TopAbs_FACE );
- for ( int i = 0; expface.More(); expface.Next(), i++ ) {
- tabID[i] = 0;
- aShape = expface.Current();
- shapeID = theMesh->ShapeToIndex( aShape );
- idFound = false;
- for ( int j=0; j<=i; j++) {
- if ( shapeID == tabID[j] ) {
- idFound = true;
- break;
- }
- }
- if ( ! idFound ) {
- tabID[i] = shapeID;
- tabShape[i] = aShape;
+ // 2 adaptors for each face in facesMap, as a face can belong to 2 solids
+ typedef vector< StdMeshers_QuadToTriaAdaptor* > TwoAdaptors;
+ vector< TwoAdaptors > qttaByFace;
+ if ( theQuad2Trias.empty() )
+ {
+ // case w/o quadrangles
+ for ( int i = 1; i <= facesMap.Extent(); ++i )
+ if (( theSubMesh = theMesh->MeshElements( facesMap(i))))
+ nbTriangles += theSubMesh->NbElements();
+ }
+ else
+ {
+ // case with quadrangles
+ qttaByFace.resize( facesMap.Extent() );
+ for ( unsigned i = 0; i < theQuad2Trias.size(); ++i )
+ {
+ TopoDS_Shape solid = theQuad2Trias[i].GetShape();
+ TopExp_Explorer expface( solid, TopAbs_FACE );
+ for ( ; expface.More(); expface.Next() )
+ if (( theSubMesh = theMesh->MeshElements( expface.Current()) ))
+ {
+ const int faceIndex = facesMap.Add( expface.Current() );
+ TwoAdaptors& aTwoAdaptors = qttaByFace[ faceIndex-1 ];
+ const bool newFaceEncounters = aTwoAdaptors.empty();
+ aTwoAdaptors.push_back( & theQuad2Trias[i] );
+
+ // on a shared face encountered for the second time
+ // we count only triangles of pyramids
+ const int countTrias = int( newFaceEncounters );
+ itOnSubMesh = theSubMesh->GetElements();
+ while ( itOnSubMesh->more() )
+ {
+ aFace = itOnSubMesh->next();
+ if ( aFace->NbCornerNodes() != 4 )
+ nbTriangles += countTrias;
+ else if ( TTriaList* trias = theQuad2Trias[i].GetTriangles( aFace ))
+ nbTriangles += trias->size();
+ else
+ nbTriangles += countTrias;
+ }
+ }
}
}
- for ( int i =0; i < nbShape; i++ ) {
- if ( tabID[i] != 0 ) {
- aShape = tabShape[i];
- shapeID = tabID[i];
- theSubMesh = theMesh->MeshElements( aShape );
- if ( !theSubMesh ) continue;
- itOnSubMesh = theSubMesh->GetElements();
- while ( itOnSubMesh->more() ) {
- aFace = itOnSubMesh->next();
+
+ std::cout << " " << facesMap.Extent() << " shapes of 2D dimension" << std::endl;
+ std::cout << std::endl;
+
+ theFile << space << nbTriangles << space << dummyint << std::endl;
+
+ vector< const SMDS_MeshElement* > trias;
+ trias.resize( 8 ); // 4 triangles from 2 pyramids basing on one quadranle
+
+ for ( int i = 1; i <= facesMap.Extent(); i++ )
+ {
+ aShape = facesMap(i);
+ theSubMesh = theMesh->MeshElements( aShape );
+ if ( !theSubMesh ) continue;
+ TwoAdaptors& aTwoAdaptors = qttaByFace[ i-1 ];
+ itOnSubMesh = theSubMesh->GetElements();
+ while ( itOnSubMesh->more() )
+ {
+ aFace = itOnSubMesh->next();
+ if ( aFace->NbCornerNodes() == 4 )
+ {
+ nbTriangles = 0;
+ for ( unsigned j = 0; j < aTwoAdaptors.size(); ++j )
+ if ( TTriaList* t = aTwoAdaptors[j]->GetTriangles( aFace ))
+ for ( TTriaList::const_iterator tIt = t->begin(); tIt != t->end(); ++tIt)
+ trias[nbTriangles++] = *tIt;
+ if ( nbTriangles == 0 )
+ {
+ nbTriangles = 1;
+ trias[0] = aFace;
+ }
+ }
+ else
+ {
+ nbTriangles = 1;
+ trias[0] = aFace;
+ }
+ for ( int j = 0; j < nbTriangles; ++j )
+ {
+ aFace = trias[j];
nbNodes = aFace->NbNodes();
theFile << space << nbNodes;
// find GHS3D ID
aSmdsID = itOnSubFace->next()->GetID();
itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
+// if ( itOnMap == theSmdsToGhs3dIdMap.end() ) {
+// cout << "not found node: " << aSmdsID << endl;
+// return false;
+// }
ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
theFile << space << (*itOnMap).second;
for ( int j=0; j<=nbNodes; j++)
theFile << space << dummyint;
- theFile << endl;
+ theFile << std::endl;
}
}
}
-
- delete [] tabID;
- delete [] tabShape;
-
+
return true;
}
//purpose : Write Faces in case if generate 3D mesh w/o geometry
//=======================================================================
-static bool writeFaces (ofstream & theFile,
- SMESHDS_Mesh * theMesh,
+static bool writeFaces (ofstream & theFile,
+ SMESHDS_Mesh * theMesh,
+ StdMeshers_QuadToTriaAdaptor& theQuad2Trias,
vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
{
// record structure:
// Loop from 1 to NB_ELEMS
// NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
-
- int nbFaces = 0;
- list< const SMDS_MeshElement* > faces;
- list< const SMDS_MeshElement* >::iterator f;
+ int nbNodes, nbTriangles = 0;
map< const SMDS_MeshNode*,int >::iterator it;
SMDS_ElemIteratorPtr nodeIt;
const SMDS_MeshElement* elem;
- int nbNodes;
const char* space = " ";
const int dummyint = 0;
- //get all faces from mesh
- SMDS_FaceIteratorPtr eIt = theMesh->facesIterator();
- while ( eIt->more() ) {
- const SMDS_MeshElement* elem = eIt->next();
- if ( !elem )
- return false;
- faces.push_back( elem );
- nbFaces++;
- }
+ // count faces
- if ( nbFaces == 0 )
+ nbTriangles = theQuad2Trias.TotalNbOfTriangles();
+ if ( nbTriangles == 0 )
+ // theQuad2Trias not computed as there are no quadrangles in mesh
+ nbTriangles = theMesh->NbFaces();
+
+ if ( nbTriangles == 0 )
return false;
-
- cout << "The initial 2D mesh contains " << nbFaces << " faces and ";
+
+ std::cout << "The initial 2D mesh contains " << nbTriangles << " faces and ";
// NB_ELEMS DUMMY_INT
- theFile << space << nbFaces << space << dummyint << endl;
+ theFile << space << nbTriangles << space << dummyint << std::endl;
- // Loop from 1 to NB_ELEMS
map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
- f = faces.begin();
- for ( ; f != faces.end(); ++f )
+
+ vector< const SMDS_MeshElement* > trias;
+ trias.resize( 8 ); // 8 == 4 triangles from 2 pyramids basing on one quadranle
+
+ // Loop from 1 to NB_ELEMS
+
+ SMDS_FaceIteratorPtr eIt = theMesh->facesIterator();
+ while ( eIt->more() )
{
- // NB_NODES PER FACE
- elem = *f;
- nbNodes = elem->NbNodes();
- theFile << space << nbNodes;
-
- // NODE_NB_1 NODE_NB_2 ...
- nodeIt = elem->nodesIterator();
- while ( nodeIt->more() )
+ elem = eIt->next();
+ // get triangles
+ if ( elem->NbCornerNodes() == 4 )
{
- // find GHS3D ID
- const SMDS_MeshNode* node = castToNode( nodeIt->next() );
- int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
- it = aNodeToGhs3dIdMap.insert( make_pair( node, newId )).first;
- theFile << space << it->second;
+ nbTriangles = 0;
+ if ( TTriaList* t = theQuad2Trias.GetTriangles( elem ))
+ for ( TTriaList::const_iterator tIt = t->begin(); tIt != t->end(); ++tIt)
+ trias[nbTriangles++] = *tIt;
+ if ( nbTriangles == 0 )
+ {
+ nbTriangles = 1;
+ trias[0] = elem;
+ }
}
+ else
+ {
+ nbTriangles = 1;
+ trias[0] = elem;
+ }
+ // write triangles
+ for ( int j = 0; j < nbTriangles; ++j )
+ {
+ // NB_NODES PER FACE
+ elem = trias[j];
+ nbNodes = elem->NbNodes();
+ theFile << space << nbNodes;
+
+ // NODE_NB_1 NODE_NB_2 ...
+ nodeIt = elem->nodesIterator();
+ while ( nodeIt->more() )
+ {
+ // find GHS3D ID
+ const SMDS_MeshNode* node = castToNode( nodeIt->next() );
+ int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
+ it = aNodeToGhs3dIdMap.insert( make_pair( node, newId )).first;
+ theFile << space << it->second;
+ }
- // (NB_NODES + 1) times: DUMMY_INT
- for ( int i=0; i<=nbNodes; i++)
- theFile << space << dummyint;
-
- theFile << endl;
+ // (NB_NODES + 1) times: DUMMY_INT
+ for ( int i=0; i<=nbNodes; i++)
+ theFile << space << dummyint;
+ theFile << std::endl;
+ }
}
// put nodes to theNodeByGhs3dId vector
//=======================================================================
static bool writePoints (ofstream & theFile,
- SMESHDS_Mesh * theMesh,
+ SMESH_MesherHelper& theHelper,
map <int,int> & theSmdsToGhs3dIdMap,
- map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap)
+ map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap,
+ map<vector<double>,double> & theEnforcedVertices)
{
// record structure:
//
// Loop from 1 to NB_NODES
// X Y Z DUMMY_INT
+ SMESHDS_Mesh * theMesh = theHelper.GetMeshDS();
int nbNodes = theMesh->NbNodes();
if ( nbNodes == 0 )
return false;
+ int nbEnforcedVertices = theEnforcedVertices.size();
+ // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
+ // The problem is in nodes on degenerated edges, we need to skip them
+ if ( theHelper.HasDegeneratedEdges() )
+ {
+ // here we decrease total nb of nodes by nb of nodes on degenerated edges
+ set<int> checkedSM;
+ for (TopExp_Explorer e(theMesh->ShapeToMesh(), TopAbs_EDGE ); e.More(); e.Next())
+ {
+ SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh( e.Current() );
+ if ( checkedSM.insert( sm->GetId() ).second && theHelper.IsDegenShape(sm->GetId() )) {
+ if ( sm->GetSubMeshDS() )
+ nbNodes -= sm->GetSubMeshDS()->NbNodes();
+ }
+ }
+ }
const char* space = " ";
const int dummyint = 0;
const SMDS_MeshNode* node;
// NB_NODES
- theFile << space << nbNodes << endl;
- cout << endl;
- cout << "The initial 2D mesh contains :" << endl;
- cout << " " << nbNodes << " vertices" << endl;
+ std::cout << std::endl;
+ std::cout << "The initial 2D mesh contains :" << std::endl;
+ std::cout << " " << nbNodes << " nodes" << std::endl;
+ if (nbEnforcedVertices > 0) {
+ std::cout << " " << nbEnforcedVertices << " enforced vertices" << std::endl;
+ }
+ std::cout << std::endl;
+ std::cout << "Start writing in 'points' file ..." << std::endl;
+ theFile << space << nbNodes << std::endl;
// Loop from 1 to NB_NODES
while ( it->more() )
{
node = it->next();
+ if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
+ theHelper.IsDegenShape( node->GetPosition()->GetShapeId() )) // Issue 020674
+ continue;
+
theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node ));
aGhs3dID++;
// X Y Z DUMMY_INT
theFile
- << space << node->X()
- << space << node->Y()
- << space << node->Z()
- << space << dummyint;
+ << space << node->X()
+ << space << node->Y()
+ << space << node->Z()
+ << space << dummyint;
+
+ theFile << std::endl;
- theFile << endl;
}
+
+ // Iterate over the enforced vertices
+ GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
+ const TopoDS_Shape shapeToMesh = theMesh->ShapeToMesh();
+ for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
+ double x = vertexIt->first[0];
+ double y = vertexIt->first[1];
+ double z = vertexIt->first[2];
+ // Test if point is inside shape to mesh
+ gp_Pnt myPoint(x,y,z);
+ BRepClass3d_SolidClassifier scl(shapeToMesh);
+ scl.Perform(myPoint, 1e-7);
+ TopAbs_State result = scl.State();
+ if ( result == TopAbs_IN ) {
+ MESSAGE("Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second);
+ // X Y Z PHY_SIZE DUMMY_INT
+ theFile
+ << space << x
+ << space << y
+ << space << z
+ << space << vertexIt->second
+ << space << dummyint;
+
+ theFile << std::endl;
+ }
+ else
+ MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added ");
+ }
+
+
+ std::cout << std::endl;
+ std::cout << "End writing in 'points' file." << std::endl;
return true;
}
//=======================================================================
static bool writePoints (ofstream & theFile,
- SMESHDS_Mesh * theMesh,
- const vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
+ SMESH_Mesh * theMesh,
+ const vector <const SMDS_MeshNode*> & theNodeByGhs3dId,
+ map<vector<double>,double> & theEnforcedVertices)
{
// record structure:
//
// NB_NODES
// Loop from 1 to NB_NODES
// X Y Z DUMMY_INT
-
- //int nbNodes = theMesh->NbNodes();
+
int nbNodes = theNodeByGhs3dId.size();
if ( nbNodes == 0 )
return false;
+ int nbEnforcedVertices = theEnforcedVertices.size();
+
const char* space = " ";
const int dummyint = 0;
const SMDS_MeshNode* node;
// NB_NODES
- theFile << space << nbNodes << endl;
- cout << nbNodes << " nodes" << endl;
+ std::cout << std::endl;
+ std::cout << "The initial 2D mesh contains :" << std::endl;
+ std::cout << " " << nbNodes << " nodes" << std::endl;
+ std::cout << " " << nbEnforcedVertices << " enforced vertices" << std::endl;
+ std::cout << std::endl;
+ std::cout << "Start writing in 'points' file ..." << std::endl;
+ theFile << space << nbNodes << std::endl;
// Loop from 1 to NB_NODES
// X Y Z DUMMY_INT
theFile
- << space << node->X()
- << space << node->Y()
- << space << node->Z()
- << space << dummyint;
+ << space << node->X()
+ << space << node->Y()
+ << space << node->Z()
+ << space << dummyint;
+
+ theFile << std::endl;
- theFile << endl;
}
+
+ // Iterate over the enforced vertices
+ GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt;
+ auto_ptr< SMESH_ElementSearcher > pntCls ( SMESH_MeshEditor( theMesh ).GetElementSearcher());
+ for(vertexIt = theEnforcedVertices.begin() ; vertexIt != theEnforcedVertices.end() ; ++vertexIt) {
+ double x = vertexIt->first[0];
+ double y = vertexIt->first[1];
+ double z = vertexIt->first[2];
+ // Test if point is inside shape to mesh
+ gp_Pnt myPoint(x,y,z);
+ TopAbs_State result = pntCls->GetPointState( myPoint );
+ if ( result == TopAbs_IN ) {
+ std::cout << "Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second << std::endl;
+
+ // X Y Z PHY_SIZE DUMMY_INT
+ theFile
+ << space << x
+ << space << y
+ << space << z
+ << space << vertexIt->second
+ << space << dummyint;
+
+ theFile << std::endl;
+ }
+ }
+ std::cout << std::endl;
+ std::cout << "End writing in 'points' file." << std::endl;
return true;
}
const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
if ( !face )
return invalidID;
-
+#ifdef _DEBUG_
+ std::cout << "bnd face " << face->GetID() << " - ";
+#endif
// geom face the face assigned to
SMESH_MeshEditor editor(&mesh);
int geomFaceID = editor.FindShape( face );
if ( meshNormal.SquareMagnitude() < DBL_MIN )
return invalidID;
- // find UV of node1 on geomFace
+ // get normale to geomFace at any node
+ bool geomNormalOK = false;
+ gp_Vec geomNormal;
+ const SMDS_MeshNode* nodes[3] = { node1, node2, node3 };
SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace );
- const SMDS_MeshNode* nNotOnSeamEdge = 0;
- if ( helper.IsSeamShape( node1->GetPosition()->GetShapeId() ))
- if ( helper.IsSeamShape( node2->GetPosition()->GetShapeId() ))
- nNotOnSeamEdge = node3;
- else
- nNotOnSeamEdge = node2;
- gp_XY uv = helper.GetNodeUV( geomFace, node1, nNotOnSeamEdge );
- // check that uv is correct
- BRepAdaptor_Surface surface( geomFace );
- double tol = BRep_Tool::Tolerance( geomFace );
- if ( node1Pnt.Distance( surface.Value( uv.X(), uv.Y() )) > 2 * tol ) {
- GeomAPI_ProjectPointOnSurf projector( node1Pnt, surface.Surface().Surface(), 2 * tol );
- if ( !projector.IsDone() || projector.NbPoints() < 1 || projector.LowerDistance() > 2 * tol)
- return invalidID;
- Quantity_Parameter U,V;
- projector.LowerDistanceParameters(U,V);
- if ( node1Pnt.Distance( surface.Value( U, V )) > 2 * tol )
- return invalidID;
- uv.SetCoord( U,V );
+ for ( int i = 0; !geomNormalOK && i < 3; ++i )
+ {
+ // find UV of i-th node on geomFace
+ const SMDS_MeshNode* nNotOnSeamEdge = 0;
+ if ( helper.IsSeamShape( nodes[i]->GetPosition()->GetShapeId() ))
+ if ( helper.IsSeamShape( nodes[(i+1)%3]->GetPosition()->GetShapeId() ))
+ nNotOnSeamEdge = nodes[(i+2)%3];
+ else
+ nNotOnSeamEdge = nodes[(i+1)%3];
+ bool uvOK;
+ gp_XY uv = helper.GetNodeUV( geomFace, nodes[i], nNotOnSeamEdge, &uvOK );
+ // check that uv is correct
+ if (uvOK) {
+ double tol = 1e-6;
+ TopoDS_Shape nodeShape = helper.GetSubShapeByNode( nodes[i], meshDS );
+ if ( !nodeShape.IsNull() )
+ switch ( nodeShape.ShapeType() )
+ {
+ case TopAbs_FACE: tol = BRep_Tool::Tolerance( TopoDS::Face( nodeShape )); break;
+ case TopAbs_EDGE: tol = BRep_Tool::Tolerance( TopoDS::Edge( nodeShape )); break;
+ case TopAbs_VERTEX: tol = BRep_Tool::Tolerance( TopoDS::Vertex( nodeShape )); break;
+ default:;
+ }
+ gp_Pnt nodePnt ( nodes[i]->X(), nodes[i]->Y(), nodes[i]->Z() );
+ BRepAdaptor_Surface surface( geomFace );
+ uvOK = ( nodePnt.Distance( surface.Value( uv.X(), uv.Y() )) < 2 * tol );
+ if ( uvOK ) {
+ // normale to geomFace at UV
+ gp_Vec du, dv;
+ surface.D1( uv.X(), uv.Y(), nodePnt, du, dv );
+ geomNormal = du ^ dv;
+ if ( geomFace.Orientation() == TopAbs_REVERSED )
+ geomNormal.Reverse();
+ geomNormalOK = ( geomNormal.SquareMagnitude() > DBL_MIN * 1e3 );
+ }
+ }
}
- // normale to geomFace at UV
- gp_Vec du, dv;
- surface.D1( uv.X(), uv.Y(), node1Pnt, du, dv );
- gp_Vec geomNormal = du ^ dv;
- if ( geomNormal.SquareMagnitude() < DBL_MIN )
- return findShapeID( mesh, node2, node3, node1, toMeshHoles );
- if ( geomFace.Orientation() == TopAbs_REVERSED )
- geomNormal.Reverse();
+ if ( !geomNormalOK)
+ return invalidID;
// compare normals
bool isReverse = ( meshNormal * geomNormal ) < 0;
else
return meshDS->ShapeToIndex( solids(2) );
}
-
+
//=======================================================================
//function : readResultFile
//purpose :
//=======================================================================
static bool readResultFile(const int fileOpen,
+#ifdef WNT
+ const char* fileName,
+#endif
SMESH_Mesh& theMesh,
TopoDS_Shape tabShape[],
double** tabBox,
const int nbShape,
map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
- bool toMeshHoles)
+ bool toMeshHoles,
+ int nbEnforcedVertices)
{
+ MESSAGE("GHS3DPlugin_GHS3D::readResultFile()");
+ Kernel_Utils::Localizer loc;
struct stat status;
size_t length;
length = status.st_size;
// Mapping the result file into memory
+#ifdef WNT
+ HANDLE fd = CreateFile(fileName, GENERIC_READ, FILE_SHARE_READ,
+ NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
+ HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
+ 0, (DWORD)length, NULL);
+ ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
+#else
ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
+#endif
mapPtr = ptr;
ptr = readMapIntLine(ptr, tab);
for (int i=0; i < 4*nbElems; i++)
nodeId = strtol(ptr, &ptr, 10);
+ MESSAGE("nbInputNodes: "<<nbInputNodes);
+ MESSAGE("nbEnforcedVertices: "<<nbEnforcedVertices);
// Reading the nodeCoor and update the nodeMap
for (int iNode=1; iNode <= nbNodes; iNode++) {
for (int iCoor=0; iCoor < 3; iCoor++)
coord[ iCoor ] = strtod(ptr, &ptr);
nodeAssigne[ iNode ] = 1;
- if ( iNode > nbInputNodes ) {
+ if ( iNode > (nbInputNodes-nbEnforcedVertices) ) {
+ // Creating SMESH nodes
+ // - for enforced vertices
+ // - for vertices of forced edges
+ // - for ghs3d nodes
nodeAssigne[ iNode ] = 0;
aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
theGhs3dIdToNodeMap.insert(theGhs3dIdToNodeMap.end(), make_pair( iNode, aNewNode ));
}
// END -- 0020330: Pb with ghs3d as a submesh
#ifdef _DEBUG_
- cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << endl;
+ std::cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << std::endl;
+#endif
+ }
+ catch ( Standard_Failure & ex)
+ {
+#ifdef _DEBUG_
+ std::cout << i+1 << " subdomain: Exception caugt: " << ex.GetMessageString() << std::endl;
+#endif
+ }
+ catch (...) {
+#ifdef _DEBUG_
+ std::cout << i+1 << " subdomain: unknown exception caught " << std::endl;
#endif
- } catch ( Standard_Failure ) {
- } catch (...) {}
+ }
}
}
node[ iNode ] = itOnNode->second;
nodeID[ iNode ] = ID;
}
- // We always run GHS3D with "to mesh holes'==TRUE but we must not create
+ // We always run GHS3D with "to mesh holes"==TRUE but we must not create
// tetras within holes depending on hypo option,
// so we first check if aTet is inside a hole and then create it
//aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
if ( nbElems )
cout << nbElems << " tetrahedrons have been associated to " << nbShape << " shapes" << endl;
+#ifdef WNT
+ UnmapViewOfFile(mapPtr);
+ CloseHandle(hMapObject);
+ CloseHandle(fd);
+#else
munmap(mapPtr, length);
+#endif
close(fileOpen);
delete [] tab;
delete [] nodeAssigne;
#ifdef _DEBUG_
+ shapeIDs.erase(-1);
if ( shapeIDs.size() != nbShape ) {
- cout << "Only " << shapeIDs.size() << " solids of " << nbShape << " found" << endl;
+ std::cout << "Only " << shapeIDs.size() << " solids of " << nbShape << " found" << std::endl;
for (int i=0; i<nbShape; i++) {
shapeID = theMeshDS->ShapeToIndex( tabShape[i] );
if ( shapeIDs.find( shapeID ) == shapeIDs.end() )
- cout << " Solid #" << shapeID << " not found" << endl;
+ std::cout << " Solid #" << shapeID << " not found" << std::endl;
}
}
#endif
//=======================================================================
static bool readResultFile(const int fileOpen,
- SMESHDS_Mesh* theMeshDS,
+#ifdef WNT
+ const char* fileName,
+#endif
+ SMESH_Mesh& theMesh,
TopoDS_Shape aSolid,
- vector <const SMDS_MeshNode*>& theNodeByGhs3dId) {
+ vector <const SMDS_MeshNode*>& theNodeByGhs3dId,
+ int nbEnforcedVertices)
+{
+ SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
+ Kernel_Utils::Localizer loc;
struct stat status;
size_t length;
length = status.st_size;
// Mapping the result file into memory
+#ifdef WNT
+ HANDLE fd = CreateFile(fileName, GENERIC_READ, FILE_SHARE_READ,
+ NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL);
+ HANDLE hMapObject = CreateFileMapping(fd, NULL, PAGE_READONLY,
+ 0, (DWORD)length, NULL);
+ ptr = ( char* ) MapViewOfFile(hMapObject, FILE_MAP_READ, 0, 0, 0 );
+#else
ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
+#endif
mapPtr = ptr;
ptr = readMapIntLine(ptr, tab);
for (int i=0; i < 4*nbElems; i++)
nodeId = strtol(ptr, &ptr, 10);
+ // Issue 0020682. Avoid creating nodes and tetras at place where
+ // volumic elements already exist
+ SMESH_ElementSearcher* elemSearcher = 0;
+ vector< const SMDS_MeshElement* > foundVolumes;
+ if ( theMesh.NbVolumes() > 0 )
+ elemSearcher = SMESH_MeshEditor( &theMesh ).GetElementSearcher();
+
// Reading the nodeCoor and update the nodeMap
shapeID = theMeshDS->ShapeToIndex( aSolid );
for (int iNode=0; iNode < nbNodes; iNode++) {
for (int iCoor=0; iCoor < 3; iCoor++)
coord[ iCoor ] = strtod(ptr, &ptr);
- if ((iNode+1) > nbInputNodes) {
- aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
- theMeshDS->SetNodeInVolume( aNewNode, shapeID );
- theNodeByGhs3dId[ iNode ] = aNewNode;
+ if ((iNode+1) > (nbInputNodes-nbEnforcedVertices)) {
+ // Issue 0020682. Avoid creating nodes and tetras at place where
+ // volumic elements already exist
+ if ( elemSearcher &&
+ elemSearcher->FindElementsByPoint( gp_Pnt(coord[0],coord[1],coord[2]),
+ SMDSAbs_Volume, foundVolumes ))
+ {
+ theNodeByGhs3dId[ iNode ] = 0;
+ }
+ else
+ {
+ aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
+ theMeshDS->SetNodeInVolume( aNewNode, shapeID );
+ theNodeByGhs3dId[ iNode ] = aNewNode;
+ }
}
}
ID = strtol(tetraPtr, &tetraPtr, 10);
node[ iNode ] = theNodeByGhs3dId[ ID-1 ];
}
+ if ( elemSearcher )
+ {
+ // Issue 0020682. Avoid creating nodes and tetras at place where
+ // volumic elements already exist
+ if ( !node[1] || !node[0] || !node[2] || !node[3] )
+ continue;
+ if ( elemSearcher->FindElementsByPoint(( SMESH_MeshEditor::TNodeXYZ(node[0]) +
+ SMESH_MeshEditor::TNodeXYZ(node[1]) +
+ SMESH_MeshEditor::TNodeXYZ(node[2]) +
+ SMESH_MeshEditor::TNodeXYZ(node[3]) ) / 4.,
+ SMDSAbs_Volume, foundVolumes ))
+ continue;
+ }
aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
shapeID = theMeshDS->ShapeToIndex( aSolid );
theMeshDS->SetMeshElementOnShape( aTet, shapeID );
}
if ( nbElems )
cout << nbElems << " tetrahedrons have been associated to " << nbTriangle << " shapes" << endl;
+#ifdef WNT
+ UnmapViewOfFile(mapPtr);
+ CloseHandle(hMapObject);
+ CloseHandle(fd);
+#else
munmap(mapPtr, length);
+#endif
close(fileOpen);
delete [] tab;
tabBox[i] = new double[6];
Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
- // TopExp_Explorer expBox (meshDS->ShapeToMesh(), TopAbs_SOLID); -- 0020330:...ghs3d as a submesh
for (expBox.ReInit(); expBox.More(); expBox.Next()) {
tabShape[iShape] = expBox.Current();
Bnd_Box BoundingBox;
ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
Ok =
-#ifdef WIN32
- aFacesFile->is_open() && aPointsFile->is_open();
-#else
aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
-#endif
if (!Ok) {
INFOS( "Can't write into " << aFacesFileName);
return error(SMESH_Comment("Can't write into ") << aFacesFileName);
}
map <int,int> aSmdsToGhs3dIdMap;
map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
+ map<vector<double>,double> enforcedVertices;
+ int nbEnforcedVertices = 0;
+ try {
+ enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
+ nbEnforcedVertices = enforcedVertices.size();
+ }
+ catch(...) {
+ }
+
+ SMESH_MesherHelper helper( theMesh );
+ helper.SetSubShape( theShape );
+
+ // make prisms on quadrangles
+ vector<StdMeshers_QuadToTriaAdaptor> aQuad2Trias;
+ if ( theMesh.NbQuadrangles() > 0 )
+ {
+ aQuad2Trias.resize( _nbShape );
+ for (_iShape = 0, expBox.ReInit(); expBox.More(); expBox.Next())
+ aQuad2Trias[ _iShape++ ].Compute( theMesh, expBox.Current() );
+ }
+
+ Ok = (writePoints( aPointsFile, helper, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, enforcedVertices) &&
+ writeFaces ( aFacesFile, meshDS, theShape, aQuad2Trias, aSmdsToGhs3dIdMap ));
- Ok = writePoints( aPointsFile, meshDS, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap ) &&
- writeFaces ( aFacesFile, meshDS, aSmdsToGhs3dIdMap );
+ // Write aSmdsToGhs3dIdMap to temp file
+ TCollection_AsciiString aSmdsToGhs3dIdMapFileName;
+ aSmdsToGhs3dIdMapFileName = aGenericName + ".ids"; // ids relation
+ ofstream aIdsFile ( aSmdsToGhs3dIdMapFileName.ToCString() , ios::out);
+ Ok =
+ aIdsFile.rdbuf()->is_open();
+ if (!Ok) {
+ INFOS( "Can't write into " << aSmdsToGhs3dIdMapFileName);
+ return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName);
+ }
+ aIdsFile << "Smds Ghs3d" << std::endl;
+ map <int,int>::const_iterator myit;
+ for (myit=aSmdsToGhs3dIdMap.begin() ; myit != aSmdsToGhs3dIdMap.end() ; ++myit) {
+ aIdsFile << myit->first << " " << myit->second << std::endl;
+ }
aFacesFile.close();
aPointsFile.close();
-
+ aIdsFile.close();
+
if ( ! Ok ) {
if ( !_keepFiles ) {
- OSD_File( aFacesFileName ).Remove();
- OSD_File( aPointsFileName ).Remove();
+ removeFile( aFacesFileName );
+ removeFile( aPointsFileName );
+ removeFile( aSmdsToGhs3dIdMapFileName );
}
return error(COMPERR_BAD_INPUT_MESH);
}
- OSD_File( aResultFileName ).Remove(); // needed for boundary recovery module usage
+ removeFile( aResultFileName ); // needed for boundary recovery module usage
// -----------------
// run ghs3d mesher
cmd += TCollection_AsciiString(" -f ") + aGenericName; // file to read
cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file
- cout << endl;
- cout << "Ghs3d execution..." << endl;
- cout << cmd << endl;
+ std::cout << std::endl;
+ std::cout << "Ghs3d execution..." << std::endl;
+ std::cout << cmd << std::endl;
system( cmd.ToCString() ); // run
- cout << endl;
- cout << "End of Ghs3d execution !" << endl;
+ std::cout << std::endl;
+ std::cout << "End of Ghs3d execution !" << std::endl;
// --------------
// read a result
int fileOpen;
fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
if ( fileOpen < 0 ) {
- cout << endl;
- cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << endl;
- cout << "Log: " << aLogFileName << endl;
+ std::cout << std::endl;
+ std::cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << std::endl;
+ std::cout << "Log: " << aLogFileName << std::endl;
Ok = false;
}
else {
bool toMeshHoles =
_hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
- Ok = readResultFile( fileOpen, theMesh, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
- toMeshHoles );
+ Ok = readResultFile( fileOpen,
+#ifdef WNT
+ aResultFileName.ToCString(),
+#endif
+ theMesh, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
+ toMeshHoles, nbEnforcedVertices );
}
// ---------------------
if ( Ok )
{
if ( !_keepFiles )
- OSD_File( aLogFileName ).Remove();
+ removeFile( aLogFileName );
}
else if ( OSD_File( aLogFileName ).Size() > 0 )
{
else
{
// the log file is empty
- OSD_File( aLogFileName ).Remove();
+ removeFile( aLogFileName );
INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
}
if ( !_keepFiles ) {
- OSD_File( aFacesFileName ).Remove();
- OSD_File( aPointsFileName ).Remove();
- OSD_File( aResultFileName ).Remove();
- OSD_File( aBadResFileName ).Remove();
- OSD_File( aBbResFileName ).Remove();
+ removeFile( aFacesFileName );
+ removeFile( aPointsFileName );
+ removeFile( aResultFileName );
+ removeFile( aBadResFileName );
+ removeFile( aBbResFileName );
+ removeFile( aSmdsToGhs3dIdMapFileName );
}
- cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
+ std::cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
if ( !Ok )
- cout << "not ";
- cout << "treated !" << endl;
- cout << endl;
+ std::cout << "not ";
+ std::cout << "treated !" << std::endl;
+ std::cout << std::endl;
_nbShape = 0; // re-initializing _nbShape for the next Compute() method call
delete [] tabShape;
ofstream aFacesFile ( aFacesFileName.ToCString() , ios::out);
ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
bool Ok =
-#ifdef WIN32
- aFacesFile->is_open() && aPointsFile->is_open();
-#else
aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
-#endif
if (!Ok)
return error( SMESH_Comment("Can't write into ") << aPointsFileName);
+
+ GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices;
+ int nbEnforcedVertices = 0;
+ try {
+ enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp);
+ nbEnforcedVertices = enforcedVertices.size();
+ }
+ catch(...) {
+ }
vector <const SMDS_MeshNode*> aNodeByGhs3dId;
- Ok = (writeFaces ( aFacesFile, meshDS, aNodeByGhs3dId ) &&
- writePoints( aPointsFile, meshDS, aNodeByGhs3dId));
+ StdMeshers_QuadToTriaAdaptor aQuad2Trias;
+ if ( theMesh.NbQuadrangles() > 0 )
+ aQuad2Trias.Compute( theMesh );
+
+ Ok = (writeFaces ( aFacesFile, meshDS, aQuad2Trias, aNodeByGhs3dId ) &&
+ writePoints( aPointsFile, &theMesh, aNodeByGhs3dId,enforcedVertices));
aFacesFile.close();
aPointsFile.close();
if ( ! Ok ) {
if ( !_keepFiles ) {
- OSD_File( aFacesFileName ).Remove();
- OSD_File( aPointsFileName ).Remove();
+ removeFile( aFacesFileName );
+ removeFile( aPointsFileName );
}
return error(COMPERR_BAD_INPUT_MESH);
}
- OSD_File( aResultFileName ).Remove(); // needed for boundary recovery module usage
+ removeFile( aResultFileName ); // needed for boundary recovery module usage
// -----------------
// run ghs3d mesher
int fileOpen;
fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
if ( fileOpen < 0 ) {
- cout << endl;
- cout << "Error when opening the " << aResultFileName.ToCString() << " file" << endl;
- cout << "Log: " << aLogFileName << endl;
- cout << endl;
+ std::cout << std::endl;
+ std::cout << "Error when opening the " << aResultFileName.ToCString() << " file" << std::endl;
+ std::cout << "Log: " << aLogFileName << std::endl;
+ std::cout << std::endl;
Ok = false;
}
else {
- Ok = readResultFile( fileOpen, meshDS, theShape ,aNodeByGhs3dId );
+ Ok = readResultFile( fileOpen,
+#ifdef WNT
+ aResultFileName.ToCString(),
+#endif
+ theMesh, theShape ,aNodeByGhs3dId, nbEnforcedVertices );
}
// ---------------------
if ( Ok )
{
if ( !_keepFiles )
- OSD_File( aLogFileName ).Remove();
+ removeFile( aLogFileName );
}
else if ( OSD_File( aLogFileName ).Size() > 0 )
{
}
else {
// the log file is empty
- OSD_File( aLogFileName ).Remove();
+ removeFile( aLogFileName );
INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
}
if ( !_keepFiles )
{
- OSD_File( aFacesFileName ).Remove();
- OSD_File( aPointsFileName ).Remove();
- OSD_File( aResultFileName ).Remove();
- OSD_File( aBadResFileName ).Remove();
- OSD_File( aBbResFileName ).Remove();
+ removeFile( aFacesFileName );
+ removeFile( aPointsFileName );
+ removeFile( aResultFileName );
+ removeFile( aBadResFileName );
+ removeFile( aBbResFileName );
}
return Ok;
return 0;
}
+
+//=============================================================================
+/*!
+ *
+ */
+//=============================================================================
+bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh,
+ const TopoDS_Shape& aShape,
+ MapShapeNbElems& aResMap)
+{
+ int nbtri = 0, nbqua = 0;
+ double fullArea = 0.0;
+ for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
+ TopoDS_Face F = TopoDS::Face( exp.Current() );
+ SMESH_subMesh *sm = aMesh.GetSubMesh(F);
+ MapShapeNbElemsItr anIt = aResMap.find(sm);
+ if( anIt==aResMap.end() ) {
+ SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+ smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
+ "Submesh can not be evaluated",this));
+ return false;
+ }
+ std::vector<int> aVec = (*anIt).second;
+ nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
+ nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
+ GProp_GProps G;
+ BRepGProp::SurfaceProperties(F,G);
+ double anArea = G.Mass();
+ fullArea += anArea;
+ }
+
+ // collect info from edges
+ int nb0d_e = 0, nb1d_e = 0;
+ bool IsQuadratic = false;
+ bool IsFirst = true;
+ TopTools_MapOfShape tmpMap;
+ for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
+ TopoDS_Edge E = TopoDS::Edge(exp.Current());
+ if( tmpMap.Contains(E) )
+ continue;
+ tmpMap.Add(E);
+ SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+ MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+ std::vector<int> aVec = (*anIt).second;
+ nb0d_e += aVec[SMDSEntity_Node];
+ nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+ if(IsFirst) {
+ IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
+ IsFirst = false;
+ }
+ }
+ tmpMap.Clear();
+
+ double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
+
+ GProp_GProps G;
+ BRepGProp::VolumeProperties(aShape,G);
+ double aVolume = G.Mass();
+ double tetrVol = 0.1179*ELen*ELen*ELen;
+ double CoeffQuality = 0.9;
+ int nbVols = int(aVolume/tetrVol/CoeffQuality);
+ int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
+ int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
+ std::vector<int> aVec(SMDSEntity_Last);
+ for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
+ if( IsQuadratic ) {
+ aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
+ aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
+ aVec[SMDSEntity_Quad_Pyramid] = nbqua;
+ }
+ else {
+ aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
+ aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
+ aVec[SMDSEntity_Pyramid] = nbqua;
+ }
+ SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
+ aResMap.insert(std::make_pair(sm,aVec));
+
+ return true;
+}
+
