#include "SMESH_Mesh.hxx"
#include "SMESH_Comment.hxx"
#include "SMESH_MesherHelper.hxx"
+#include "SMESH_MeshEditor.hxx"
#include "SMDS_MeshElement.hxx"
#include "SMDS_MeshNode.hxx"
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepClass3d_SolidClassifier.hxx>
+#include <BRep_Tool.hxx>
+#include <Bnd_Box.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+#include <OSD_File.hxx>
+#include <Precision.hxx>
+#include <Quantity_Parameter.hxx>
+#include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
-#include <OSD_File.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
+#include <TopoDS.hxx>
+//#include <BRepClass_FaceClassifier.hxx>
+//#include <BRepGProp.hxx>
+//#include <GProp_GProps.hxx>
#include "utilities.h"
//#include <Standard_Stream.hxx>
-#include <BRepGProp.hxx>
-#include <BRepBndLib.hxx>
-#include <BRepClass_FaceClassifier.hxx>
-#include <BRepClass3d_SolidClassifier.hxx>
-#include <TopAbs.hxx>
-#include <Bnd_Box.hxx>
-#include <GProp_GProps.hxx>
-#include <Precision.hxx>
#define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
TopoDS_Shape aShape,
const TopoDS_Shape shape[],
double** box,
- const int nShape) {
- double *pntCoor;
- int iShape, nbNode = 4;
-
- pntCoor = new double[3];
- for ( int i=0; i<3; i++ ) {
- pntCoor[i] = 0;
- for ( int j=0; j<nbNode; j++ ) {
- if ( i == 0) pntCoor[i] += aNode[j]->X();
- if ( i == 1) pntCoor[i] += aNode[j]->Y();
- if ( i == 2) pntCoor[i] += aNode[j]->Z();
- }
- pntCoor[i] /= nbNode;
- }
+ const int nShape,
+ TopAbs_State * state = 0)
+{
+ gp_XYZ aPnt(0,0,0);
+ int j, iShape, nbNode = 4;
+
+ for ( j=0; j<nbNode; j++ )
+ aPnt += gp_XYZ( aNode[j]->X(), aNode[j]->Y(), aNode[j]->Z() );
+ aPnt /= nbNode;
- gp_Pnt aPnt(pntCoor[0], pntCoor[1], pntCoor[2]);
BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
+ if (state) *state = SC.State();
if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) {
for (iShape = 0; iShape < nShape; iShape++) {
aShape = shape[iShape];
- if ( !( pntCoor[0] < box[iShape][0] || box[iShape][1] < pntCoor[0] ||
- pntCoor[1] < box[iShape][2] || box[iShape][3] < pntCoor[1] ||
- pntCoor[2] < box[iShape][4] || box[iShape][5] < pntCoor[2]) ) {
+ if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() ||
+ aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() ||
+ aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) {
BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
+ if (state) *state = SC.State();
if (SC.State() == TopAbs_IN)
break;
}
return true;
}
+//=======================================================================
+//function : findShapeID
+//purpose : find the solid corresponding to GHS3D sub-domain following
+// the technique proposed in GHS3D manual in chapter
+// "B.4 Subdomain (sub-region) assignment"
+//=======================================================================
+
+static int findShapeID(SMESH_Mesh& mesh,
+ const SMDS_MeshNode* node1,
+ const SMDS_MeshNode* node2,
+ const SMDS_MeshNode* node3)
+{
+ const int invalidID = 0;
+ SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
+
+ // face th enodes belong to
+ const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
+ if ( !face )
+ return invalidID;
+
+ // geom face the face assigned to
+ SMESH_MeshEditor editor(&mesh);
+ int geomFaceID = editor.FindShape( face );
+ if ( !geomFaceID )
+ return invalidID;
+ TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
+ if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
+ return invalidID;
+ TopoDS_Face geomFace = TopoDS::Face( shape );
+
+ // solids bounded by geom face
+ TopTools_IndexedMapOfShape solids;
+ TopTools_ListIteratorOfListOfShape ansIt = mesh.GetAncestors(geomFace);
+ for ( ; ansIt.More(); ansIt.Next() ) {
+ if ( ansIt.Value().ShapeType() == TopAbs_SOLID )
+ solids.Add( ansIt.Value() );
+ }
+ // analyse found solids
+ if ( solids.Extent() == 0 )
+ return invalidID;
+ if ( solids.Extent() == 1 )
+ return meshDS->ShapeToIndex( solids(1) );
+
+ // find orientation of geom face within the first solid
+ TopoDS_Shape solid1 = solids(1);
+ TopExp_Explorer fExp( solid1, TopAbs_FACE );
+ for ( ; fExp.More(); fExp.Next() )
+ if ( geomFace.IsSame( fExp.Current() )) {
+ geomFace = TopoDS::Face( fExp.Current() );
+ break;
+ }
+ if ( !fExp.More() )
+ return invalidID; // face not found
+
+ // find UV of node1 on geomFace
+ SMESH_MesherHelper helper( mesh );
+ gp_XY uv = helper.GetNodeUV( geomFace, node1 );
+
+ // check that uv is correct
+ gp_Pnt node1Pnt ( node1->X(), node1->Y(), node1->Z() );
+ double tol = BRep_Tool::Tolerance( geomFace );
+ BRepAdaptor_Surface surface( geomFace );
+ if ( node1Pnt.Distance( surface.Value( uv.X(), uv.Y() )) > 2 * tol ) {
+ // project node1 onto geomFace to get right UV
+ GeomAPI_ProjectPointOnSurf projector( node1Pnt, surface.Surface().Surface() );
+ if ( !projector.IsDone() || projector.NbPoints() < 1 )
+ return invalidID;
+ Quantity_Parameter U,V;
+ projector.LowerDistanceParameters(U,V);
+ uv = gp_XY( U,V );
+ }
+ // normale to face at node1
+ gp_Pnt node2Pnt ( node2->X(), node2->Y(), node2->Z() );
+ gp_Pnt node3Pnt ( node3->X(), node3->Y(), node3->Z() );
+ gp_Vec vec12( node1Pnt, node2Pnt );
+ gp_Vec vec13( node1Pnt, node3Pnt );
+ gp_Vec meshNormal = vec12 ^ vec13;
+ if ( meshNormal.SquareMagnitude() < DBL_MIN )
+ return invalidID;
+
+ // 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 );
+ if ( geomFace.Orientation() == TopAbs_REVERSED )
+ geomNormal.Reverse();
+
+ // compare normals
+ bool isReverse = ( meshNormal * geomNormal ) < 0;
+ if ( isReverse )
+ return meshDS->ShapeToIndex( solids(2) );
+ else
+ return meshDS->ShapeToIndex( solid1 );
+}
+
//=======================================================================
//function : readResultFile
//purpose :
//=======================================================================
static bool readResultFile(const int fileOpen,
- SMESHDS_Mesh* theMeshDS,
+ SMESH_Mesh& theMesh,
TopoDS_Shape tabShape[],
double** tabBox,
const int nbShape,
- map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap) {
-
+ map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap,
+ bool toMeshHoles)
+{
struct stat status;
size_t length;
char *tetraPtr;
char *shapePtr;
+ SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
+
int fileStat;
int nbElems, nbNodes, nbInputNodes;
- int nodeId, triangleId;
+ int nodeId/*, triangleId*/;
int nbTriangle;
int ID, shapeID, ghs3dShapeID;
int IdShapeRef = 1;
const SMDS_MeshNode **node;
tab = new int[3];
- tabID = new int[nbShape];
+ //tabID = new int[nbShape];
nodeID = new int[4];
coord = new double[3];
node = new const SMDS_MeshNode*[4];
SMDS_MeshNode * aNewNode;
map <int,const SMDS_MeshNode*>::iterator itOnNode;
SMDS_MeshElement* aTet;
-
- for (int i=0; i<nbShape; i++)
- tabID[i] = 0;
+#ifdef _DEBUG_
+ set<int> shapeIDs;
+#endif
// Read the file state
fileStat = fstat(fileOpen, &status);
nodeAssigne = new int[ nbNodes+1 ];
+ if (nbShape > 0)
+ aSolid = tabShape[0];
+
// Reading the nodeId
for (int i=0; i < 4*nbElems; i++)
nodeId = strtol(ptr, &ptr, 10);
// Reading the nodeCoor and update the nodeMap
- for (int iNode=0; iNode < nbNodes; iNode++) {
+ for (int iNode=1; iNode <= nbNodes; iNode++) {
for (int iCoor=0; iCoor < 3; iCoor++)
coord[ iCoor ] = strtod(ptr, &ptr);
- nodeAssigne[ iNode+1 ] = 1;
- if ( (iNode+1) > nbInputNodes ) {
- nodeAssigne[ iNode+1 ] = 0;
+ nodeAssigne[ iNode ] = 1;
+ if ( iNode > nbInputNodes ) {
+ nodeAssigne[ iNode ] = 0;
aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
- theGhs3dIdToNodeMap.insert(make_pair( (iNode+1), aNewNode ));
+ theGhs3dIdToNodeMap.insert(theGhs3dIdToNodeMap.end(), make_pair( iNode, aNewNode ));
}
}
- // Reading the number of triangles which corresponds to the number of shapes
+ // Reading the number of triangles which corresponds to the number of sub-domains
nbTriangle = strtol(ptr, &ptr, 10);
- for (int i=0; i < 3*nbTriangle; i++)
- triangleId = strtol(ptr, &ptr, 10);
+ tabID = new int[nbTriangle];
+ for (int i=0; i < nbTriangle; i++) {
+ tabID[i] = 0;
+ // find the solid corresponding to GHS3D sub-domain following
+ // the technique proposed in GHS3D manual in chapter
+ // "B.4 Subdomain (sub-region) assignment"
+ int nodeId1 = strtol(ptr, &ptr, 10);
+ int nodeId2 = strtol(ptr, &ptr, 10);
+ int nodeId3 = strtol(ptr, &ptr, 10);
+ const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ];
+ const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ];
+ const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
+ try {
+ OCC_CATCH_SIGNALS;
+ tabID[i] = findShapeID( theMesh, n1, n2, n3 );
+#ifdef _DEBUG_
+ cout << i+1 << " subdomain: findShapeID() returns " << tabID[i] << endl;
+#endif
+ } catch ( Standard_Failure ) {
+ } catch (...) {}
+ }
shapePtr = ptr;
+ if ( nbTriangle <= nbShape ) // no holes
+ toMeshHoles = true; // not avoid creating tetras in holes
+
// Associating the tetrahedrons to the shapes
shapeID = compoundID;
for (int iElem = 0; iElem < nbElems; iElem++) {
node[ iNode ] = itOnNode->second;
nodeID[ iNode ] = ID;
}
- aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
- if ( nbShape > 1 ) {
- if ( nbTriangle > 1 ) {
- ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef;
- if ( tabID[ ghs3dShapeID ] == 0 ) {
- if (iElem == 0)
- aSolid = tabShape[0];
- aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
+ // 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 ( nbTriangle > 1 ) {
+ shapeID = -1; // negative shapeID means not to create tetras if !toMeshHoles
+ ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef;
+ if ( tabID[ ghs3dShapeID ] == 0 ) {
+ TopAbs_State state;
+ aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state);
+ if ( toMeshHoles || state == TopAbs_IN )
shapeID = theMeshDS->ShapeToIndex( aSolid );
- tabID[ ghs3dShapeID ] = shapeID;
- }
- else
- shapeID = tabID[ ghs3dShapeID ];
+ tabID[ ghs3dShapeID ] = shapeID;
}
- else {
- // Case where nbTriangle == 1 while nbShape == 2 encountered
- // with compound of 2 boxes and "To mesh holes"==False,
- // so there are no subdomains specified for each tetrahedron.
- // Try to guess a solid by a node already bound to shape
- shapeID = 0;
- for ( int i=0; i<4 && shapeID==0; i++ ) {
- if ( nodeAssigne[ nodeID[i] ] == 1 &&
- node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
- node[i]->GetPosition()->GetShapeId() > 1 )
- {
- shapeID = node[i]->GetPosition()->GetShapeId();
- }
- }
- if ( shapeID==0 ) {
- aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
- shapeID = theMeshDS->ShapeToIndex( aSolid );
+ else
+ shapeID = tabID[ ghs3dShapeID ];
+ }
+ else if ( nbShape > 1 ) {
+ // Case where nbTriangle == 1 while nbShape == 2 encountered
+ // with compound of 2 boxes and "To mesh holes"==False,
+ // so there are no subdomains specified for each tetrahedron.
+ // Try to guess a solid by a node already bound to shape
+ shapeID = 0;
+ for ( int i=0; i<4 && shapeID==0; i++ ) {
+ if ( nodeAssigne[ nodeID[i] ] == 1 &&
+ node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE &&
+ node[i]->GetPosition()->GetShapeId() > 1 )
+ {
+ shapeID = node[i]->GetPosition()->GetShapeId();
}
}
+ if ( shapeID==0 ) {
+ aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
+ shapeID = theMeshDS->ShapeToIndex( aSolid );
+ }
}
- // set new nodes and tetrahedron on to the shape
+ // set new nodes and tetrahedron onto the shape
for ( int i=0; i<4; i++ ) {
if ( nodeAssigne[ nodeID[i] ] == 0 ) {
- theMeshDS->SetNodeInVolume( node[i], shapeID );
- nodeAssigne[ nodeID[i] ] = 1;
+ if ( shapeID > 0 )
+ theMeshDS->SetNodeInVolume( node[i], shapeID );
+ nodeAssigne[ nodeID[i] ] = shapeID;
}
}
- theMeshDS->SetMeshElementOnShape( aTet, shapeID );
+ if ( toMeshHoles || shapeID > 0 ) {
+ aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
+ theMeshDS->SetMeshElementOnShape( aTet, shapeID );
+ }
+#ifdef _DEBUG_
+ shapeIDs.insert( shapeID );
+#endif
+ }
+ // Remove nodes of tetras inside holes if !toMeshHoles
+ if ( !toMeshHoles ) {
+ itOnNode = theGhs3dIdToNodeMap.find( nbInputNodes );
+ for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) {
+ ID = itOnNode->first;
+ if ( nodeAssigne[ ID ] < 0 )
+ theMeshDS->RemoveFreeNode( itOnNode->second, 0 );
+ }
}
+
if ( nbElems )
cout << nbElems << " tetrahedrons have been associated to " << nbShape << " shapes" << endl;
munmap(mapPtr, length);
delete [] node;
delete [] nodeAssigne;
+#ifdef _DEBUG_
+ if ( shapeIDs.size() != nbShape ) {
+ cout << "Only " << shapeIDs.size() << " solids of " << nbShape << " found" << 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;
+ }
+ }
+#endif
+
return true;
}
}
return error(COMPERR_BAD_INPUT_MESH);
}
+ OSD_File( aResultFileName ).Remove(); // needed for boundary recovery module usage
// -----------------
// run ghs3d mesher
cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << endl;
Ok = false;
}
- else
- Ok = readResultFile( fileOpen, meshDS, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap );
+ else {
+ bool toMeshHoles =
+ _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
+ Ok = readResultFile( fileOpen, theMesh, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
+ toMeshHoles );
+ }
// ---------------------
// remove working files
}
return error(COMPERR_BAD_INPUT_MESH);
}
+ OSD_File( aResultFileName ).Remove(); // needed for boundary recovery module usage
// -----------------
// run ghs3d mesher
// -----------------
- TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() );
+ TCollection_AsciiString cmd =
+ (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str();
cmd += TCollection_AsciiString(" -f ") + aGenericName; // file to read
cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file
