//purpose :
//=======================================================================
-static bool writeFaces (ofstream & theFile,
- SMESHDS_Mesh* theMesh,
- const TopoDS_Shape& theShape,
- vector<StdMeshers_QuadToTriaAdaptor>& theQuad2Trias,
- StdMeshers_ProxyMesh::Ptr& theCdfMesh,
- const map <int,int> & theSmdsToGhs3dIdMap)
+static bool writeFaces (ofstream & theFile,
+ const StdMeshers_ProxyMesh& theMesh,
+ const TopoDS_Shape& theShape,
+ const map <int,int> & theSmdsToGhs3dIdMap)
{
// record structure:
//
TopTools_IndexedMapOfShape facesMap;
TopExp::MapShapes( theShape, TopAbs_FACE, facesMap );
- // 2 adaptors for each face in facesMap, as a face can belong to 2 solids
- typedef vector< StdMeshers_QuadToTriaAdaptor* > TwoAdaptors;
- vector< TwoAdaptors > qttaByFace;
- if ( theCdfMesh )
- {
- for ( int i = 1; i <= facesMap.Extent(); ++i )
- if (( theSubMesh = theCdfMesh->GetSubMesh( facesMap(i))))
- nbTriangles += theSubMesh->NbElements();
- }
- else 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 = 1; i <= facesMap.Extent(); ++i )
+ if (( theSubMesh = theMesh.GetSubMesh( facesMap(i))))
+ nbTriangles += theSubMesh->NbElements();
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 = theCdfMesh ? theCdfMesh->GetSubMesh( aShape ) : theMesh->MeshElements( aShape );
+ theSubMesh = theMesh.GetSubMesh(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;
-
- itOnSubFace = aFace->nodesIterator();
- while ( itOnSubFace->more() ) {
- // 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;
- }
+ nbNodes = aFace->NbNodes();
- // (NB_NODES + 1) times: DUMMY_INT
- for ( int j=0; j<=nbNodes; j++)
- theFile << space << dummyint;
+ theFile << space << nbNodes;
- theFile << std::endl;
+ itOnSubFace = aFace->nodesIterator();
+ while ( itOnSubFace->more() ) {
+ // 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;
}
+
+ // (NB_NODES + 1) times: DUMMY_INT
+ for ( int j=0; j<=nbNodes; j++)
+ theFile << space << dummyint;
+
+ theFile << std::endl;
}
}
-
+
return true;
}
//=======================================================================
static bool writeFaces (ofstream & theFile,
- SMESHDS_Mesh * theMesh,
- StdMeshers_QuadToTriaAdaptor& theQuad2Trias,
+ const StdMeshers_ProxyMesh& theMesh,
vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
{
// record structure:
const int dummyint = 0;
// count faces
-
- nbTriangles = theQuad2Trias.TotalNbOfTriangles();
- if ( nbTriangles == 0 )
- // theQuad2Trias not computed as there are no quadrangles in mesh
- nbTriangles = theMesh->NbFaces();
+ nbTriangles = theMesh.NbFaces();
if ( nbTriangles == 0 )
return false;
map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
- 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();
+ SMDS_ElemIteratorPtr eIt = theMesh.GetFaces();
while ( eIt->more() )
{
elem = eIt->next();
- // get triangles
- if ( elem->NbCornerNodes() == 4 )
- {
- 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
+ // NB_NODES PER FACE
+ nbNodes = elem->NbNodes();
+ theFile << space << nbNodes;
+
+ // NODE_NB_1 NODE_NB_2 ...
+ nodeIt = elem->nodesIterator();
+ while ( nodeIt->more() )
{
- nbTriangles = 1;
- trias[0] = elem;
+ // 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;
}
- // 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 << std::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
return true;
}
+//================================================================================
+/*!
+ * \brief returns true if a triangle defined by the nodes is a temporary face on a
+ * side facet of pyramid and defines sub-domian inside the pyramid
+ */
+//================================================================================
+
+static bool isTmpFace(const SMDS_MeshNode* node1,
+ const SMDS_MeshNode* node2,
+ const SMDS_MeshNode* node3)
+{
+ // find a pyramid sharing the 3 nodes
+ //const SMDS_MeshElement* pyram = 0;
+ SMDS_ElemIteratorPtr vIt1 = node1->GetInverseElementIterator(SMDSAbs_Volume);
+ while ( vIt1->more() )
+ {
+ const SMDS_MeshElement* pyram = vIt1->next();
+ if ( pyram->NbCornerNodes() != 5 ) continue;
+ int i2, i3;
+ if ( (i2 = pyram->GetNodeIndex( node2 )) >= 0 &&
+ (i3 = pyram->GetNodeIndex( node3 )) >= 0 )
+ {
+ // Triangle defines sub-domian inside the pyramid if it's
+ // normal points out of the pyram
+
+ // make i2 and i3 hold indices of base nodes of the pyram while
+ // keeping the nodes order in the triangle
+ const int iApex = 4;
+ if ( i2 == iApex )
+ i2 = i3, i3 = pyram->GetNodeIndex( node1 );
+ else if ( i3 == iApex )
+ i3 = i2, i2 = pyram->GetNodeIndex( node1 );
+
+ int i3base = (i2+1) % 4; // next index after i2 within the pyramid base
+ return ( i3base != i3 );
+ }
+ }
+ return false;
+}
+
//=======================================================================
//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"
+// the technique proposed in GHS3D manual (available within
+// ghs3d installation) in chapter "B.4 Subdomain (sub-region) assignment".
+// In brief: normal of the triangle defined by the given nodes
+// points out of the domain it is associated to
//=======================================================================
static int findShapeID(SMESH_Mesh& mesh,
// face the nodes belong to
const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3);
if ( !face )
- return invalidID;
+ return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
#ifdef _DEBUG_
std::cout << "bnd face " << face->GetID() << " - ";
#endif
SMESH_MeshEditor editor(&mesh);
int geomFaceID = editor.FindShape( face );
if ( !geomFaceID )
- return invalidID;
+ return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID;
TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID );
if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE )
return invalidID;
if ( theMesh.NbVolumes() > 0 )
elemSearcher = SMESH_MeshEditor( &theMesh ).GetElementSearcher();
- // Reading the nodeCoor and update the nodeMap
+ // Reading the nodeCoord and update the nodeMap
shapeID = theMeshDS->ShapeToIndex( aSolid );
for (int iNode=0; iNode < nbNodes; iNode++) {
for (int iCoor=0; iCoor < 3; iCoor++)
const TopoDS_Shape& theShape)
{
bool Ok(false);
- SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
+ //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
// we count the number of shapes
// _nbShape = countShape( meshDS, TopAbs_SOLID ); -- 0020330: Pb with ghs3d as a submesh
SMESH_MesherHelper helper( theMesh );
helper.SetSubShape( theShape );
+ StdMeshers_ProxyMesh::Ptr proxyMesh( new StdMeshers_ProxyMesh( theMesh ));
+
// 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() );
+ vector<StdMeshers_ProxyMesh::Ptr> components;
+ for (expBox.ReInit(); expBox.More(); expBox.Next())
+ {
+ if ( _viscousLayersHyp )
+ {
+ proxyMesh = _viscousLayersHyp->Compute( theMesh, expBox.Current() );
+ if ( !proxyMesh )
+ return false;
+ }
+ StdMeshers_QuadToTriaAdaptor* q2t = new StdMeshers_QuadToTriaAdaptor;
+ q2t->Compute( theMesh, expBox.Current(), proxyMesh.get() );
+ components.push_back( StdMeshers_ProxyMesh::Ptr( q2t ));
+ }
+ proxyMesh.reset( new StdMeshers_ProxyMesh( components ));
}
- // viscous layers
- StdMeshers_ProxyMesh::Ptr cdfMesh;
- if ( _viscousLayersHyp && aQuad2Trias.empty() )
+ // build viscous layers
+ else if ( _viscousLayersHyp )
{
- cdfMesh = _viscousLayersHyp->Compute( theMesh, theShape );
- if ( !cdfMesh )
+ proxyMesh = _viscousLayersHyp->Compute( theMesh, theShape );
+ if ( !proxyMesh )
return false;
}
Ok = (writePoints( aPointsFile, helper, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, enforcedVertices) &&
- writeFaces ( aFacesFile, meshDS, theShape, aQuad2Trias, cdfMesh, aSmdsToGhs3dIdMap ));
+ writeFaces ( aFacesFile, *proxyMesh, theShape, aSmdsToGhs3dIdMap ));
// Write aSmdsToGhs3dIdMap to temp file
TCollection_AsciiString aSmdsToGhs3dIdMapFileName;
{
MESSAGE("GHS3DPlugin_GHS3D::Compute()");
- SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
+ //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
TopoDS_Shape theShape = aHelper->GetSubShape();
// a unique working file name
vector <const SMDS_MeshNode*> aNodeByGhs3dId;
- StdMeshers_QuadToTriaAdaptor aQuad2Trias;
+ StdMeshers_ProxyMesh::Ptr proxyMesh( new StdMeshers_ProxyMesh( theMesh ));
if ( theMesh.NbQuadrangles() > 0 )
- aQuad2Trias.Compute( theMesh );
+ {
+ StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor;
+ aQuad2Trias->Compute( theMesh );
+ proxyMesh.reset( aQuad2Trias );
+ }
- Ok = (writeFaces ( aFacesFile, meshDS, aQuad2Trias, aNodeByGhs3dId ) &&
+ Ok = (writeFaces ( aFacesFile, *proxyMesh, aNodeByGhs3dId ) &&
writePoints( aPointsFile, &theMesh, aNodeByGhs3dId,enforcedVertices));
aFacesFile.close();
