From: gdd Date: Mon, 2 May 2011 16:52:27 +0000 (+0000) Subject: Merge from BR_ENFORCED_MESH X-Git-Tag: V6_3_0~28 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=2a26ed174871bbe62f3077de4c7faa5d41e86d34;p=plugins%2Fhybridplugin.git Merge from BR_ENFORCED_MESH --- diff --git a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx index 48caf14..c726643 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx +++ b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx @@ -27,14 +27,15 @@ #include "GHS3DPlugin_GHS3D.hxx" #include "GHS3DPlugin_Hypothesis.hxx" - #include -#include "SMESH_Gen.hxx" +//#include "SMESH_Gen.hxx" +#include #include "SMESH_Mesh.hxx" #include "SMESH_Comment.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_MeshEditor.hxx" +#include "SMESH_OctreeNode.hxx" #include "SMDS_MeshElement.hxx" #include "SMDS_MeshNode.hxx" @@ -46,7 +47,9 @@ #include #include +#include #include +#include #include #include #include @@ -74,6 +77,7 @@ #else #include #endif +#include using namespace std; @@ -101,6 +105,10 @@ extern "C" #define HOLE_ID -1 +#ifndef GHS3D_VERSION +#define GHS3D_VERSION 41 +#endif + typedef const list TTriaList; static void removeFile( const TCollection_AsciiString& fileName ) @@ -179,6 +187,7 @@ bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh& aMesh, return true; } + //======================================================================= //function : findShape //purpose : @@ -238,394 +247,6 @@ static char* readMapIntLine(char* ptr, int tab[]) { return ptr; } -//======================================================================= -//function : countShape -//purpose : -//======================================================================= - -// template < class Mesh, class Shape > -// static int countShape( Mesh* mesh, Shape shape ) { -// TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape ); -// int nbShape = 0; -// for ( ; expShape.More(); expShape.Next() ) { -// nbShape++; -// } -// return nbShape; -// } - -//======================================================================= -//function : writeFaces -//purpose : -//======================================================================= - -static bool writeFaces (ofstream & theFile, - const SMESH_ProxyMesh& theMesh, - const TopoDS_Shape& theShape, - const map & theSmdsToGhs3dIdMap) -{ - // record structure: - // - // NB_ELEMS DUMMY_INT - // Loop from 1 to NB_ELEMS - // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT - - TopoDS_Shape aShape; - const SMESHDS_SubMesh* theSubMesh; - const SMDS_MeshElement* aFace; - const char* space = " "; - const int dummyint = 0; - map::const_iterator itOnMap; - SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace; - int nbNodes, aSmdsID; - - // count triangles bound to geometry - int nbTriangles = 0; - - TopTools_IndexedMapOfShape facesMap; - TopExp::MapShapes( theShape, TopAbs_FACE, facesMap ); - - 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; - - for ( int i = 1; i <= facesMap.Extent(); i++ ) - { - aShape = facesMap(i); - theSubMesh = theMesh.GetSubMesh(aShape); - if ( !theSubMesh ) continue; - itOnSubMesh = theSubMesh->GetElements(); - while ( itOnSubMesh->more() ) - { - aFace = itOnSubMesh->next(); - nbNodes = aFace->NbCornerNodes(); - - theFile << space << nbNodes; - - itOnSubFace = aFace->nodesIterator(); - while ( itOnSubFace->more() && 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; - } - - // (NB_NODES + 1) times: DUMMY_INT - nbNodes = aFace->NbCornerNodes(); - for ( int j=0; j<=nbNodes; j++) - theFile << space << dummyint; - - theFile << std::endl; - } - } - - return true; -} - -//======================================================================= -//function : writeFaces -//purpose : Write Faces in case if generate 3D mesh w/o geometry -//======================================================================= - -static bool writeFaces (ofstream & theFile, - const SMESH_ProxyMesh& theMesh, - vector & theNodeByGhs3dId) -{ - // record structure: - // - // NB_ELEMS DUMMY_INT - // Loop from 1 to NB_ELEMS - // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT - - int nbNodes, nbTriangles = 0; - map< const SMDS_MeshNode*,int >::iterator it; - SMDS_ElemIteratorPtr nodeIt; - const SMDS_MeshElement* elem; - - const char* space = " "; - const int dummyint = 0; - - // count faces - nbTriangles = theMesh.NbFaces(); - - if ( nbTriangles == 0 ) - return false; - - std::cout << "The initial 2D mesh contains " << nbTriangles << " faces and "; - - // NB_ELEMS DUMMY_INT - theFile << space << nbTriangles << space << dummyint << std::endl; - - - map aNodeToGhs3dIdMap; - - // Loop from 1 to NB_ELEMS - - SMDS_ElemIteratorPtr eIt = theMesh.GetFaces(); - while ( eIt->more() ) - { - elem = eIt->next(); - // NB_NODES PER FACE - nbNodes = elem->NbCornerNodes(); - theFile << space << nbNodes; - - // NODE_NB_1 NODE_NB_2 ... - nodeIt = elem->nodesIterator(); - while ( nodeIt->more() && nbNodes-- ) - { - // 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 - nbNodes = elem->NbCornerNodes(); - for ( int i=0; i<=nbNodes; i++) - theFile << space << dummyint; - theFile << std::endl; - } - - // put nodes to theNodeByGhs3dId vector - theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() ); - map::const_iterator n2id = aNodeToGhs3dIdMap.begin(); - for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id) - { - theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1 - } - - return true; -} - -//======================================================================= -//function : writePoints -//purpose : -//======================================================================= - -static bool writePoints (ofstream & theFile, - SMESH_MesherHelper& theHelper, - map & theSmdsToGhs3dIdMap, - map & theGhs3dIdToNodeMap, - map,double> & theEnforcedVertices) -{ - // record structure: - // - // NB_NODES - // 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(); - - int aGhs3dID = 1; - SMDS_NodeIteratorPtr it = theMesh->nodesIterator(); - const SMDS_MeshNode* node; - - // 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 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 bool isQuadMesh = - theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) || - theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) || - theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC ); - if ( isQuadMesh ) - { - // descrease nbNodes by nb of medium nodes - while ( it->more() ) - { - node = it->next(); - if ( !theHelper.IsDegenShape( node->getshapeId() )) - nbNodes -= int( theHelper.IsMedium( node )); - } - it = theMesh->nodesIterator(); - } - - const char* space = " "; - const int dummyint = 0; - - // NB_NODES - 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 (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238 - theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674 - 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; - - theFile << std::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; -} - -//======================================================================= -//function : writePoints -//purpose : -//======================================================================= - -static bool writePoints (ofstream & theFile, - SMESH_Mesh * theMesh, - const vector & theNodeByGhs3dId, - map,double> & theEnforcedVertices) -{ - // record structure: - // - // NB_NODES - // Loop from 1 to NB_NODES - // X Y Z DUMMY_INT - - 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 - 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 - - vector::const_iterator nodeIt = theNodeByGhs3dId.begin(); - vector::const_iterator after = theNodeByGhs3dId.end(); - for ( ; nodeIt != after; ++nodeIt ) - { - node = *nodeIt; - - // X Y Z DUMMY_INT - theFile - << space << node->X() - << space << node->Y() - << space << node->Z() - << space << dummyint; - - theFile << std::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; -} - //================================================================================ /*! * \brief returns true if a triangle defined by the nodes is a temporary face on a @@ -789,34 +410,2029 @@ static int findShapeID(SMESH_Mesh& mesh, 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 ); + 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 ); + } + } + } + if ( !geomNormalOK) + return invalidID; + + // compare normals + bool isReverse = ( meshNormal * geomNormal ) < 0; + if ( !isReverse ) + return meshDS->ShapeToIndex( solid1 ); + + if ( solids.Extent() == 1 ) + return HOLE_ID; // we are inside a hole + else + return meshDS->ShapeToIndex( solids(2) ); +} + +//======================================================================= +//function : countShape +//purpose : +//======================================================================= + +// template < class Mesh, class Shape > +// static int countShape( Mesh* mesh, Shape shape ) { +// TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape ); +// TopTools_MapOfShape mapShape; +// int nbShape = 0; +// for ( ; expShape.More(); expShape.Next() ) { +// if (mapShape.Add(expShape.Current())) { +// nbShape++; +// } +// } +// return nbShape; +// } + +//======================================================================= +//function : getShape +//purpose : +//======================================================================= + +// template < class Mesh, class Shape, class Tab > +// void getShape(Mesh* mesh, Shape shape, Tab *t_Shape) { +// TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape ); +// TopTools_MapOfShape mapShape; +// for ( int i=0; expShape.More(); expShape.Next() ) { +// if (mapShape.Add(expShape.Current())) { +// t_Shape[i] = expShape.Current(); +// i++; +// } +// } +// return; +// } + +// //======================================================================= +// //function : findEdgeID +// //purpose : +// //======================================================================= +// +// static int findEdgeID(const SMDS_MeshNode* aNode, +// const SMESHDS_Mesh* theMesh, +// const int nEdge, +// const TopoDS_Shape* t_Edge) { +// +// TopoDS_Shape aPntShape, foundEdge; +// TopoDS_Vertex aVertex; +// gp_Pnt aPnt( aNode->X(), aNode->Y(), aNode->Z() ); +// +// int foundInd, ind; +// double nearest = RealLast(), *t_Dist; +// double epsilon = Precision::Confusion(); +// +// t_Dist = new double[ nEdge ]; +// aPntShape = BRepBuilderAPI_MakeVertex( aPnt ).Shape(); +// aVertex = TopoDS::Vertex( aPntShape ); +// +// for ( ind=0; ind < nEdge; ind++ ) { +// BRepExtrema_DistShapeShape aDistance ( aVertex, t_Edge[ind] ); +// t_Dist[ind] = aDistance.Value(); +// if ( t_Dist[ind] < nearest ) { +// nearest = t_Dist[ind]; +// foundEdge = t_Edge[ind]; +// foundInd = ind; +// if ( nearest < epsilon ) +// ind = nEdge; +// } +// } +// +// delete [] t_Dist; +// return theMesh->ShapeToIndex( foundEdge ); +// } + + +//======================================================================= +//function : readGMFFile +//purpose : read GMF file with geometry associated to mesh +// TODO +//======================================================================= + +// static bool readGMFFile( +// const int fileOpen, +// const char* theFileName, +// SMESH_Mesh& theMesh, +// const int nbShape, +// const TopoDS_Shape* tabShape, +// double** tabBox, +// map & theGhs3dIdToNodeMap, +// bool toMeshHoles, +// int nbEnforcedVertices, +// int nbEnforcedNodes, +// TIDSortedNodeSet & theEnforcedNodes, +// TIDSortedElemSet & theEnforcedTriangles, +// TIDSortedElemSet & theEnforcedQuadrangles) +// { +// TopoDS_Shape aShape; +// TopoDS_Vertex aVertex; +// SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS(); +// int nbElem = 0, nbRef = 0, IdShapeRef = 1; +// int *tabID; +// int aGMFNodeID = 0; +// int compoundID = +// nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() ); +// int tetraShapeID = compoundID; +// double epsilon = Precision::Confusion(); +// int *nodeAssigne, *GMFNodeAssigne; +// SMDS_MeshNode** GMFNode; +// TopoDS_Shape *tabCorner, *tabEdge; +// std::map tabRef; +// +// +// int ver, dim; +// MESSAGE("Read " << theFileName << " file"); +// int InpMsh = GmfOpenMesh(theFileName, GmfRead, &ver, &dim); +// if (!InpMsh) +// return false; +// +// // =========================== +// // Fill the tabID array: BEGIN +// // =========================== +// +// /* +// The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) +// */ +// Kernel_Utils::Localizer loc; +// struct stat status; +// size_t length; +// +// char *ptr, *mapPtr; +// char *tetraPtr; +// int *tab = new int[3]; +// +// // Read the file state +// fstat(fileOpen, &status); +// length = status.st_size; +// +// // Mapping the result file into memory +// #ifdef WNT +// HANDLE fd = CreateFile(theFileName, 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); +// tetraPtr = ptr; +// +// nbElem = tab[0]; +// int nbNodes = tab[1]; +// +// for (int i=0; i < 4*nbElem; i++) +// strtol(ptr, &ptr, 10); +// +// for (int iNode=1; iNode <= nbNodes; iNode++) +// for (int iCoor=0; iCoor < 3; iCoor++) +// strtod(ptr, &ptr); +// +// +// // Reading the number of triangles which corresponds to the number of sub-domains +// int nbTriangle = strtol(ptr, &ptr, 10); +// +// +// // The keyword does not exist yet => to update when it is created +// // int nbTriangle = GmfStatKwd(InpMsh, GmfSubdomain); +// // int id_tri[3]; +// +// +// tabID = new int[nbTriangle]; +// for (int i=0; i < nbTriangle; i++) { +// tabID[i] = 0; +// int nodeId1, nodeId2, nodeId3; +// // find the solid corresponding to GHS3D sub-domain following +// // the technique proposed in GHS3D manual in chapter +// // "B.4 Subdomain (sub-region) assignment" +// +// nodeId1 = strtol(ptr, &ptr, 10); +// nodeId2 = strtol(ptr, &ptr, 10); +// nodeId3 = strtol(ptr, &ptr, 10); +// +// // // The keyword does not exist yet => to update when it is created +// // GmfGetLin(InpMsh, GmfSubdomain, &id_tri[0], &id_tri[1], &id_tri[2]); +// // nodeId1 = id_tri[0]; +// // nodeId2 = id_tri[1]; +// // nodeId3 = id_tri[2]; +// +// if ( nbTriangle > 1 ) { +// // get the nodes indices +// 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, toMeshHoles ); +// // -- 0020330: Pb with ghs3d as a submesh +// // check that found shape is to be meshed +// if ( tabID[i] > 0 ) { +// const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( tabID[i] ); +// bool isToBeMeshed = false; +// for ( int iS = 0; !isToBeMeshed && iS < nbShape; ++iS ) +// isToBeMeshed = foundShape.IsSame( tabShape[ iS ]); +// if ( !isToBeMeshed ) +// tabID[i] = HOLE_ID; +// } +// // END -- 0020330: Pb with ghs3d as a submesh +// #ifdef _DEBUG_ +// 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 +// } +// } +// } +// +// // =========================== +// // Fill the tabID array: END +// // =========================== +// +// +// tabRef[GmfVertices] = 3; +// tabRef[GmfCorners] = 1; +// tabRef[GmfEdges] = 2; +// tabRef[GmfRidges] = 1; +// tabRef[GmfTriangles] = 3; +// // tabRef[GmfQuadrilaterals] = 4; +// tabRef[GmfTetrahedra] = 4; +// // tabRef[GmfHexahedra] = 8; +// +// SMDS_NodeIteratorPtr itOnGMFInputNode = theMeshDS->nodesIterator(); +// while ( itOnGMFInputNode->more() ) +// theMeshDS->RemoveNode( itOnGMFInputNode->next() ); +// +// +// int nbVertices = GmfStatKwd(InpMsh, GmfVertices); +// int nbCorners = max(countShape( theMeshDS, TopAbs_VERTEX ) , GmfStatKwd(InpMsh, GmfCorners)); +// int nbShapeEdge = countShape( theMeshDS, TopAbs_EDGE ); +// +// tabCorner = new TopoDS_Shape[ nbCorners ]; +// tabEdge = new TopoDS_Shape[ nbShapeEdge ]; +// nodeAssigne = new int[ nbVertices + 1 ]; +// GMFNodeAssigne = new int[ nbVertices + 1 ]; +// GMFNode = new SMDS_MeshNode*[ nbVertices + 1 ]; +// +// getShape(theMeshDS, TopAbs_VERTEX, tabCorner); +// getShape(theMeshDS, TopAbs_EDGE, tabEdge); +// +// std::map ::const_iterator it = tabRef.begin(); +// for ( ; it != tabRef.end() ; ++it) +// { +// // int dummy; +// GmfKwdCod token = it->first; +// nbRef = it->second; +// +// nbElem = GmfStatKwd(InpMsh, token); +// if (nbElem > 0) { +// GmfGotoKwd(InpMsh, token); +// std::cout << "Read " << nbElem; +// } +// else +// continue; +// +// int id[nbElem*tabRef[token]]; +// int ghs3dShapeID[nbElem]; +// +// if (token == GmfVertices) { +// std::cout << " vertices" << std::endl; +// int aGMFID; +// +// float VerTab_f[nbElem][3]; +// double VerTab_d[nbElem][3]; +// SMDS_MeshNode * aGMFNode; +// +// for ( int iElem = 0; iElem < nbElem; iElem++ ) { +// aGMFID = iElem + 1; +// if (ver == GmfFloat) { +// GmfGetLin(InpMsh, token, &VerTab_f[nbElem][0], &VerTab_f[nbElem][1], &VerTab_f[nbElem][2], &ghs3dShapeID[iElem]); +// aGMFNode = theMeshDS->AddNode(VerTab_f[nbElem][0], VerTab_f[nbElem][1], VerTab_f[nbElem][2]); +// } +// else { +// GmfGetLin(InpMsh, token, &VerTab_d[nbElem][0], &VerTab_d[nbElem][1], &VerTab_d[nbElem][2], &ghs3dShapeID[iElem]); +// aGMFNode = theMeshDS->AddNode(VerTab_d[nbElem][0], VerTab_d[nbElem][1], VerTab_d[nbElem][2]); +// } +// GMFNode[ aGMFID ] = aGMFNode; +// nodeAssigne[ aGMFID ] = 0; +// GMFNodeAssigne[ aGMFID ] = 0; +// } +// } +// else if (token == GmfCorners && nbElem > 0) { +// std::cout << " corners" << std::endl; +// for ( int iElem = 0; iElem < nbElem; iElem++ ) +// GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]); +// } +// else if (token == GmfRidges && nbElem > 0) { +// std::cout << " ridges" << std::endl; +// for ( int iElem = 0; iElem < nbElem; iElem++ ) +// GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]); +// } +// else if (token == GmfEdges && nbElem > 0) { +// std::cout << " edges" << std::endl; +// for ( int iElem = 0; iElem < nbElem; iElem++ ) +// GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &ghs3dShapeID[iElem]); +// } +// else if (token == GmfTriangles && nbElem > 0) { +// std::cout << " triangles" << std::endl; +// for ( int iElem = 0; iElem < nbElem; iElem++ ) +// GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &ghs3dShapeID[iElem]); +// } +// // else if (token == GmfQuadrilaterals && nbElem > 0) { +// // std::cout << " Quadrilaterals" << std::endl; +// // for ( int iElem = 0; iElem < nbElem; iElem++ ) +// // GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &ghs3dShapeID[iElem]); +// // } +// else if (token == GmfTetrahedra && nbElem > 0) { +// std::cout << " Tetrahedra" << std::endl; +// for ( int iElem = 0; iElem < nbElem; iElem++ ) +// GmfGetLin(InpMsh, token, +// &id[iElem*tabRef[token]], +// &id[iElem*tabRef[token]+1], +// &id[iElem*tabRef[token]+2], +// &id[iElem*tabRef[token]+3], +// &ghs3dShapeID[iElem]); +// } +// // else if (token == GmfHexahedra && nbElem > 0) { +// // std::cout << " Hexahedra" << std::endl; +// // for ( int iElem = 0; iElem < nbElem; iElem++ ) +// // GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], +// // &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &ghs3dShapeID[iElem]); +// // } +// +// switch (token) { +// case GmfCorners: +// case GmfRidges: +// case GmfEdges: +// case GmfTriangles: +// // case GmfQuadrilaterals: +// case GmfTetrahedra: +// // case GmfHexahedra: +// { +// int nodeDim, shapeID, *nodeID; +// SMDS_MeshNode** node; +// // std::vector< SMDS_MeshNode* > enfNode( nbRef ); +// SMDS_MeshElement * aGMFElement; +// +// node = new SMDS_MeshNode*[nbRef]; +// nodeID = new int[ nbRef ]; +// +// for ( int iElem = 0; iElem < nbElem; iElem++ ) +// { +// for ( int iRef = 0; iRef < nbRef; iRef++ ) +// { +// aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID +// node [ iRef ] = GMFNode[ aGMFNodeID ]; +// nodeID[ iRef ] = aGMFNodeID; +// } +// +// switch (token) +// { +// case GmfCorners: { +// nodeDim = 1; +// gp_Pnt GMFPnt ( node[0]->X(), node[0]->Y(), node[0]->Z() ); +// for ( int i=0; iAddEdge( node[0], node[1] ); +// int iNode = 1; +// if ( GMFNodeAssigne[ nodeID[0] ] == 0 || GMFNodeAssigne[ nodeID[0] ] == 2 ) +// iNode = 0; +// shapeID = findEdgeID( node[iNode], theMeshDS, nbShapeEdge, tabEdge ); +// break; +// } +// case GmfRidges: +// break; +// case GmfTriangles: { +// nodeDim = 3; +// aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2]); +// shapeID = -1; +// break; +// } +// // case GmfQuadrilaterals: { +// // nodeDim = 4; +// // aGMFElement = theMeshDS->AddFace( node[0], node[1], node[2], node[3] ); +// // shapeID = -1; +// // break; +// // } +// case GmfTetrahedra: { +// +// // IN WORK +// TopoDS_Shape aSolid; +// // 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 +// if ( nbTriangle > 1 ) { +// tetraShapeID = HOLE_ID; // negative tetraShapeID means not to create tetras if !toMeshHoles +// int aGhs3dShapeID = ghs3dShapeID[iElem] - IdShapeRef; +// if ( tabID[ aGhs3dShapeID ] == 0 ) { +// TopAbs_State state; +// aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state); +// if ( toMeshHoles || state == TopAbs_IN ) +// tetraShapeID = theMeshDS->ShapeToIndex( aSolid ); +// tabID[ aGhs3dShapeID ] = tetraShapeID; +// } +// else +// tetraShapeID = tabID[ aGhs3dShapeID ]; +// } +// 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 +// tetraShapeID = 0; +// for ( int i=0; i<4 && tetraShapeID==0; i++ ) { +// if ( nodeAssigne[ nodeID[i] ] == 1 && +// node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE && +// node[i]->getshapeId() > 1 ) +// { +// tetraShapeID = node[i]->getshapeId(); +// } +// } +// if ( tetraShapeID==0 ) { +// aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape); +// tetraShapeID = theMeshDS->ShapeToIndex( aSolid ); +// } +// } +// // set new nodes and tetrahedron onto the shape +// for ( int i=0; i<4; i++ ) { +// if ( nodeAssigne[ nodeID[i] ] == 0 ) { +// if ( tetraShapeID != HOLE_ID ) +// theMeshDS->SetNodeInVolume( node[i], tetraShapeID ); +// nodeAssigne[ nodeID[i] ] = tetraShapeID; +// } +// } +// if ( toMeshHoles || tetraShapeID != HOLE_ID ) { +// aGMFElement = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] ); +// theMeshDS->SetMeshElementOnShape( aGMFElement, tetraShapeID ); +// } +// +// // IN WORK +// +// nodeDim = 5; +// break; +// } +// // case GmfHexahedra: { +// // nodeDim = 6; +// // aGMFElement = theMeshDS->AddVolume( node[0], node[3], node[2], node[1], +// // node[4], node[7], node[6], node[5] ); +// // break; +// // } +// default: continue; +// } +// if (token != GmfRidges) +// { +// for ( int i=0; iSetNodeOnVertex( node[0], aVertex ); +// else if ( token == GmfEdges ) theMeshDS->SetNodeOnEdge( node[i], shapeID ); +// else if ( token == GmfTriangles ) theMeshDS->SetNodeOnFace( node[i], shapeID ); +// GMFNodeAssigne[ nodeID[i] ] = nodeDim; +// } +// } +// if ( token != "Corners" ) +// theMeshDS->SetMeshElementOnShape( aGMFElement, shapeID ); +// } +// } // for +// +// if ( !toMeshHoles ) { +// map ::iterator itOnNode = theGhs3dIdToNodeMap.find( nbVertices-(nbEnforcedVertices+nbEnforcedNodes) ); +// for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) { +// if ( nodeAssigne[ itOnNode->first ] == HOLE_ID ) +// theMeshDS->RemoveFreeNode( itOnNode->second, 0 ); +// } +// } +// +// delete [] node; +// delete [] nodeID; +// break; +// } // case GmfTetrahedra +// } // switch(token) +// } // for +// cout << std::endl; +// +// #ifdef WNT +// UnmapViewOfFile(mapPtr); +// CloseHandle(hMapObject); +// CloseHandle(fd); +// #else +// munmap(mapPtr, length); +// #endif +// close(fileOpen); +// +// delete [] tabID; +// delete [] tabCorner; +// delete [] tabEdge; +// delete [] nodeAssigne; +// delete [] GMFNodeAssigne; +// delete [] GMFNode; +// +// return true; +// } + + +//======================================================================= +//function : readGMFFile +//purpose : read GMF file w/o geometry associated to mesh +//======================================================================= + + +static bool readGMFFile(const char* theFile, +#ifdef WITH_SMESH_CANCEL_COMPUTE + GHS3DPlugin_GHS3D* theAlgo, +#endif + SMESH_MesherHelper* theHelper, + TIDSortedNodeSet & theEnforcedNodes, + TIDSortedElemSet & theEnforcedTriangles, + TIDSortedElemSet & theEnforcedQuadrangles) +{ + SMESHDS_Mesh* theMesh = theHelper->GetMeshDS(); + + // --------------------------------- + // Read generated elements and nodes + // --------------------------------- + + int nbElem = 0, nbRef = 0; + int aGMFNodeID = 0, shapeID; + int *nodeAssigne; + SMDS_MeshNode** GMFNode; + std::map tabRef; + + tabRef[GmfVertices] = 3; + tabRef[GmfCorners] = 1; + tabRef[GmfEdges] = 2; + tabRef[GmfRidges] = 1; + tabRef[GmfTriangles] = 3; + tabRef[GmfQuadrilaterals] = 4; + tabRef[GmfTetrahedra] = 4; + tabRef[GmfHexahedra] = 8; + + theHelper->GetMesh()->Clear(); + + int ver, dim; + MESSAGE("Read " << theFile << " file"); + int InpMsh = GmfOpenMesh(theFile, GmfRead, &ver, &dim); + if (!InpMsh) + return false; + + int nbVertices = GmfStatKwd(InpMsh, GmfVertices); + GMFNode = new SMDS_MeshNode*[ nbVertices + 1 ]; + nodeAssigne = new int[ nbVertices + 1 ]; + + std::map ::const_iterator it = tabRef.begin(); + for ( ; it != tabRef.end() ; ++it) + { +#ifdef WITH_SMESH_CANCEL_COMPUTE + if(theAlgo->computeCanceled()) { + GmfCloseMesh(InpMsh); + delete [] GMFNode; + delete [] nodeAssigne; + return false; + } +#endif + int dummy; + GmfKwdCod token = it->first; + nbRef = it->second; + + nbElem = GmfStatKwd(InpMsh, token); + if (nbElem > 0) { + GmfGotoKwd(InpMsh, token); + std::cout << "Read " << nbElem; + } + else + continue; + + int id[nbElem*tabRef[token]]; + + if (token == GmfVertices) { + std::cout << " vertices" << std::endl; + int aGMFID; + + float VerTab_f[nbElem][3]; + double VerTab_d[nbElem][3]; + SMDS_MeshNode * aGMFNode; + + for ( int iElem = 0; iElem < nbElem; iElem++ ) { +#ifdef WITH_SMESH_CANCEL_COMPUTE + if(theAlgo->computeCanceled()) { + GmfCloseMesh(InpMsh); + delete [] GMFNode; + delete [] nodeAssigne; + return false; + } +#endif + aGMFID = iElem + 1; + if (ver == GmfFloat) { + GmfGetLin(InpMsh, token, &VerTab_f[nbElem][0], &VerTab_f[nbElem][1], &VerTab_f[nbElem][2], &dummy); + aGMFNode = theMesh->AddNode(VerTab_f[nbElem][0], VerTab_f[nbElem][1], VerTab_f[nbElem][2]); + } + else { + GmfGetLin(InpMsh, token, &VerTab_d[nbElem][0], &VerTab_d[nbElem][1], &VerTab_d[nbElem][2], &dummy); + aGMFNode = theMesh->AddNode(VerTab_d[nbElem][0], VerTab_d[nbElem][1], VerTab_d[nbElem][2]); + } + GMFNode[ aGMFID ] = aGMFNode; + nodeAssigne[ aGMFID ] = 0; + } + } + else if (token == GmfCorners && nbElem > 0) { + std::cout << " corners" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]); + } + else if (token == GmfRidges && nbElem > 0) { + std::cout << " ridges" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]); + } + else if (token == GmfEdges && nbElem > 0) { + std::cout << " edges" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &dummy); + } + else if (token == GmfTriangles && nbElem > 0) { + std::cout << " triangles" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &dummy); + } + else if (token == GmfQuadrilaterals && nbElem > 0) { + std::cout << " Quadrilaterals" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &dummy); + } + else if (token == GmfTetrahedra && nbElem > 0) { + std::cout << " Tetrahedra" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &dummy); + } + else if (token == GmfHexahedra && nbElem > 0) { + std::cout << " Hexahedra" << std::endl; + for ( int iElem = 0; iElem < nbElem; iElem++ ) + GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], + &id[iElem*tabRef[token]+4], &id[iElem*tabRef[token]+5], &id[iElem*tabRef[token]+6], &id[iElem*tabRef[token]+7], &dummy); + } + std::cout << std::endl; + + switch (token) { + case GmfCorners: + case GmfRidges: + case GmfEdges: + case GmfTriangles: + case GmfQuadrilaterals: + case GmfTetrahedra: + case GmfHexahedra: { + std::vector< SMDS_MeshNode* > node( nbRef ); + std::vector< int > nodeID( nbRef ); + std::vector< SMDS_MeshNode* > enfNode( nbRef ); + + for ( int iElem = 0; iElem < nbElem; iElem++ ) + { +#ifdef WITH_SMESH_CANCEL_COMPUTE + if(theAlgo->computeCanceled()) { + GmfCloseMesh(InpMsh); + delete [] GMFNode; + delete [] nodeAssigne; + return false; + } +#endif + for ( int iRef = 0; iRef < nbRef; iRef++ ) + { + aGMFNodeID = id[iElem*tabRef[token]+iRef]; // read nbRef aGMFNodeID + node [ iRef ] = GMFNode[ aGMFNodeID ]; + nodeID[ iRef ] = aGMFNodeID; + } + + switch (token) + { + case GmfEdges: + theHelper->AddEdge( node[0], node[1] ); break; + case GmfTriangles: { + theMesh->AddFace( node[0], node[1], node[2]); + break; + } + case GmfQuadrilaterals: { + theMesh->AddFace( node[0], node[1], node[2], node[3] ); + break; + } + case GmfTetrahedra: + theHelper->AddVolume( node[0], node[1], node[2], node[3] ); break; + case GmfHexahedra: + theHelper->AddVolume( node[0], node[3], node[2], node[1], + node[4], node[7], node[6], node[5] ); break; + default: continue; + } + if ( token == GmfTriangles || token == GmfQuadrilaterals ) // "Quadrilaterals" and "Triangles" + for ( int iRef = 0; iRef < nbRef; iRef++ ) + nodeAssigne[ nodeID[ iRef ]] = 1; + } + break; + } + } + } + + shapeID = theHelper->GetSubShapeID(); + for ( int i = 0; i < nbVertices; ++i ) { +#ifdef WITH_SMESH_CANCEL_COMPUTE + if(theAlgo->computeCanceled()) { + GmfCloseMesh(InpMsh); + delete [] GMFNode; + delete [] nodeAssigne; + return false; + } +#endif + if ( !nodeAssigne[ i+1 ]) + theMesh->SetNodeInVolume( GMFNode[ i+1 ], shapeID ); + } + + GmfCloseMesh(InpMsh); + delete [] GMFNode; + delete [] nodeAssigne; + return true; +} + +static bool writeGMFFile(const char* theMeshFileName, + const char* theRequiredFileName, + const char* theSolFileName, + const SMESH_ProxyMesh& theProxyMesh, + SMESH_Mesh * theMesh, + vector & theNodeByGhs3dId, + vector & theEnforcedNodeByGhs3dId, + TIDSortedNodeSet & theEnforcedNodes, + TIDSortedElemSet & theEnforcedEdges, + TIDSortedElemSet & theEnforcedTriangles, + TIDSortedElemSet & theEnforcedQuadrangles, + GHS3DPlugin_Hypothesis::TEnforcedVertexValues & theEnforcedVertices) +{ + MESSAGE("writeGMFFile w/o geometry"); + int idx, idxRequired, idxSol; + const int dummyint = 0; + GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt; + std::vector enfVertexSizes; + const SMDS_MeshElement* elem; + TIDSortedElemSet anElemSet, anEnforcedEdgeSet, anEnforcedTriangleSet, anEnforcedQuadrangleSet; + SMDS_ElemIteratorPtr nodeIt; + map aNodeToGhs3dIdMap, anEnforcedNodeToGhs3dIdMap; + TIDSortedElemSet::iterator elemIt; + bool isOK; + auto_ptr< SMESH_ElementSearcher > pntCls ( SMESH_MeshEditor( theMesh ).GetElementSearcher()); + + int nbEnforcedVertices = theEnforcedVertices.size(); +// int nbEnforcedNodes = theEnforcedNodes.size(); +// int nbEnforcedEdges = theEnforcedEdges.size(); +// int nbEnforcedTriangles = theEnforcedTriangles.size(); +// int nbEnforcedQuadrangles = theEnforcedQuadrangles.size(); + + // count faces + int nbFaces = theProxyMesh.NbFaces(); + int nbNodes; + + if ( nbFaces == 0 ) + return false; + + idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + if (!idx) + return false; + + /* ========================== FACES ========================== */ + /* TRIANGLES ========================== */ + SMDS_ElemIteratorPtr eIt = theProxyMesh.GetFaces(); + while ( eIt->more() ) + { + elem = eIt->next(); + anElemSet.insert(elem); + // NODE_NB_1 NODE_NB_2 ... + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) + { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1 + aNodeToGhs3dIdMap.insert( make_pair( node, newId )); + } + } + + /* EDGES ========================== */ + + // Iterate over the enforced edges + for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) { + elem = (*elemIt); + isOK = true; + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + // Test if point is inside shape to mesh + gp_Pnt myPoint(node->X(),node->Y(),node->Z()); + TopAbs_State result = pntCls->GetPointState( myPoint ); + if ( result != TopAbs_IN ) { + isOK = false; + break; + } + } + if (isOK) { + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + int newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1 + anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId )); + } + anEnforcedEdgeSet.insert(elem); + } + } + + /* ENFORCED TRIANGLES ========================== */ + + // Iterate over the enforced triangles + for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) { + elem = (*elemIt); + isOK = true; + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + // Test if point is inside shape to mesh + gp_Pnt myPoint(node->X(),node->Y(),node->Z()); + TopAbs_State result = pntCls->GetPointState( myPoint ); + if ( result != TopAbs_IN ) { + isOK = false; + break; + } + } + if (isOK) { + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + int newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1 + anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId )); + } + anEnforcedTriangleSet.insert(elem); + } + } + + /* ENFORCED QUADRANGLES ========================== */ + + // Iterate over the enforced quadrangles + for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) { + elem = (*elemIt); + isOK = true; + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + // Test if point is inside shape to mesh + gp_Pnt myPoint(node->X(),node->Y(),node->Z()); + TopAbs_State result = pntCls->GetPointState( myPoint ); + if ( result != TopAbs_IN ) { + isOK = false; + break; + } + } + if (isOK) { + nodeIt = elem->nodesIterator(); + nbNodes = elem->NbCornerNodes(); + while ( nodeIt->more() && nbNodes-- ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + int newId = aNodeToGhs3dIdMap.size() + anEnforcedNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1 + anEnforcedNodeToGhs3dIdMap.insert( make_pair( node, newId )); + } + anEnforcedQuadrangleSet.insert(elem); + } + } + + + // put nodes to theNodeByGhs3dId vector + std::cout << "aNodeToGhs3dIdMap.size(): "<::const_iterator n2id = aNodeToGhs3dIdMap.begin(); + for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id) + { +// std::cout << "n2id->first: "<first<second - 1 ] = n2id->first; // ghs3d ids count from 1 + } + + // put nodes to anEnforcedNodeToGhs3dIdMap vector + std::cout << "anEnforcedNodeToGhs3dIdMap.size(): "<first: "<first<second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // ghs3d ids count from 1 + } + + + /* ========================== NODES ========================== */ + vector theOrderedNodes, theRequiredNodes; + std::set< std::vector > nodesCoords; + vector::const_iterator ghs3dNodeIt = theNodeByGhs3dId.begin(); + vector::const_iterator after = theNodeByGhs3dId.end(); + + std::cout << theNodeByGhs3dId.size() << " nodes from mesh ..." << std::endl; + for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt ) + { + const SMDS_MeshNode* node = *ghs3dNodeIt; + std::vector coords; + coords.push_back(node->X()); + coords.push_back(node->Y()); + coords.push_back(node->Z()); + nodesCoords.insert(coords); + theOrderedNodes.push_back(node); + } + + // Iterate over the enforced nodes + TIDSortedNodeSet::const_iterator enfNodeIt; + std::cout << theEnforcedNodes.size() << " nodes from enforced nodes ..." << std::endl; + for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt) + { + const SMDS_MeshNode* node = *enfNodeIt; + std::vector coords; + coords.push_back(node->X()); + coords.push_back(node->Y()); + coords.push_back(node->Z()); + + if (nodesCoords.find(coords) != nodesCoords.end()) { + std::cout << "Node at " << node->X()<<", " <Y()<<", " <Z() << " found" << std::endl; + continue; + } + + if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) + continue; + + // Test if point is inside shape to mesh + gp_Pnt myPoint(node->X(),node->Y(),node->Z()); + TopAbs_State result = pntCls->GetPointState( myPoint ); + if ( result != TopAbs_IN ) + continue; + + nodesCoords.insert(coords); + theOrderedNodes.push_back(node); + theRequiredNodes.push_back(node); + } + // Iterate over the enforced nodes given by enforced elements + ghs3dNodeIt = theEnforcedNodeByGhs3dId.begin(); + after = theEnforcedNodeByGhs3dId.end(); + std::cout << theEnforcedNodeByGhs3dId.size() << " nodes from enforced elements ..." << std::endl; + for ( ; ghs3dNodeIt != after; ++ghs3dNodeIt ) + { + const SMDS_MeshNode* node = *ghs3dNodeIt; + std::vector coords; + coords.push_back(node->X()); + coords.push_back(node->Y()); + coords.push_back(node->Z()); + + if (nodesCoords.find(coords) != nodesCoords.end()) { + std::cout << "Node at " << node->X()<<", " <Y()<<", " <Z() << " found" << std::endl; + continue; + } + + if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) + continue; + + nodesCoords.insert(coords); + theOrderedNodes.push_back(node); + theRequiredNodes.push_back(node); + } + + int requiredNodes = theRequiredNodes.size(); + int solSize = 0; + std::vector > ReqVerTab; + if (nbEnforcedVertices) { +// ReqVerTab.clear(); + std::cout << theEnforcedVertices.size() << " nodes from enforced vertices ..." << std::endl; + // Iterate over the enforced vertices + 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 ) + continue; + std::vector coords; + coords.push_back(x); + coords.push_back(y); + coords.push_back(z); + ReqVerTab.push_back(coords); + enfVertexSizes.push_back(vertexIt->second); + solSize++; + } + } + + // GmfVertices + std::cout << "Begin writting required nodes in GmfVertices" << std::endl; + GmfSetKwd(idx, GmfVertices, theOrderedNodes.size()+solSize); + for (ghs3dNodeIt = theOrderedNodes.begin();ghs3dNodeIt != theOrderedNodes.end();++ghs3dNodeIt) + GmfSetLin(idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint); + for (int i=0;inodesIterator(); + int index=0; + while ( nodeIt->more() ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node); + if (it == anEnforcedNodeToGhs3dIdMap.end()) + throw "Node not found"; + nedge[index] = it->second; + index++; + } + GmfSetLin(idx, GmfEdges, nedge[0], nedge[1], dummyint); +// GmfSetLin(idxRequired, GmfEdges, nedge[0], nedge[1], dummyint); + usedEnforcedEdges++; + } +// GmfCloseMesh(idxRequired); + } + + if (usedEnforcedEdges) { + GmfSetKwd(idx, GmfRequiredEdges, usedEnforcedEdges); + for (int enfID=1;enfID<=usedEnforcedEdges;enfID++) { + GmfSetLin(idx, GmfRequiredEdges, enfID); + } + } + + if (anElemSet.size()+anEnforcedTriangleSet.size()) { + GmfSetKwd(idx, GmfTriangles, anElemSet.size()+anEnforcedTriangleSet.size()); + for(elemIt = anElemSet.begin() ; elemIt != anElemSet.end() ; ++elemIt) { + elem = (*elemIt); + nodeIt = elem->nodesIterator(); + int index=0; + while ( nodeIt->more() ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + map< const SMDS_MeshNode*,int >::iterator it = aNodeToGhs3dIdMap.find(node); + if (it == aNodeToGhs3dIdMap.end()) + throw "Node not found"; + ntri[index] = it->second; + index++; + } + GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint); + } + if (anEnforcedTriangleSet.size()) { + int usedEnforcedTriangles = 0; + for(elemIt = anEnforcedTriangleSet.begin() ; elemIt != anEnforcedTriangleSet.end() ; ++elemIt) { + elem = (*elemIt); + nodeIt = elem->nodesIterator(); + int index=0; + while ( nodeIt->more() ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node); + if (it == anEnforcedNodeToGhs3dIdMap.end()) + throw "Node not found"; + ntri[index] = it->second; + index++; + } + GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint); + usedEnforcedTriangles++; + } + if (usedEnforcedTriangles) { + GmfSetKwd(idx, GmfRequiredTriangles, usedEnforcedTriangles); + for (int enfID=1;enfID<=usedEnforcedTriangles;enfID++) + GmfSetLin(idx, GmfRequiredTriangles, anElemSet.size()+enfID); + } + } + } + + if (anEnforcedQuadrangleSet.size()) { + int usedEnforcedQuadrilaterals = 0; + GmfSetKwd(idx, GmfQuadrilaterals, anEnforcedQuadrangleSet.size()); + for(elemIt = anEnforcedQuadrangleSet.begin() ; elemIt != anEnforcedQuadrangleSet.end() ; ++elemIt) { + elem = (*elemIt); + nodeIt = elem->nodesIterator(); + int index=0; + while ( nodeIt->more() ) { + // find GHS3D ID + const SMDS_MeshNode* node = castToNode( nodeIt->next() ); + map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node); + if (it == anEnforcedNodeToGhs3dIdMap.end()) + throw "Node not found"; + nquad[index] = it->second; + index++; + } + GmfSetLin(idx, GmfQuadrilaterals, nquad[0], nquad[1], nquad[2], nquad[3], dummyint); + usedEnforcedQuadrilaterals++; + } + if (usedEnforcedQuadrilaterals) { + GmfSetKwd(idx, GmfRequiredQuadrilaterals, usedEnforcedQuadrilaterals); + for (int enfID=1;enfID<=usedEnforcedQuadrilaterals;enfID++) + GmfSetLin(idx, GmfRequiredQuadrilaterals, enfID); + } + } + + GmfCloseMesh(idx); + if (idxRequired) + GmfCloseMesh(idxRequired); + if (idxSol) + GmfCloseMesh(idxSol); + + return true; + +} + +static bool writeGMFFile(const char* theMeshFileName, + const char* theRequiredFileName, + const char* theSolFileName, + SMESH_MesherHelper& theHelper, + const SMESH_ProxyMesh& theProxyMesh, + map & theSmdsToGhs3dIdMap, + map & theGhs3dIdToNodeMap, + TIDSortedNodeSet & theEnforcedNodes, + TIDSortedElemSet & theEnforcedEdges, + TIDSortedElemSet & theEnforcedTriangles, + TIDSortedElemSet & theEnforcedQuadrangles, + GHS3DPlugin_Hypothesis::TEnforcedVertexValues & theEnforcedVertices) +{ + MESSAGE("writeGMFFile with geometry"); + int idx, nbv, nbev, nben, aGhs3dID = 0; + const int dummyint = 0; + GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt; + std::vector enfVertexSizes; + TIDSortedNodeSet::const_iterator enfNodeIt; + const SMDS_MeshNode* node; + SMDS_NodeIteratorPtr nodeIt; + std::map theNodeId2NodeIndexMap; + + idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + if (!idx) + return false; + + SMESHDS_Mesh * theMeshDS = theHelper.GetMeshDS(); + + /* ========================== NODES ========================== */ + // NB_NODES + nbv = theMeshDS->NbNodes(); + if ( nbv == 0 ) + return false; + nbev = theEnforcedVertices.size(); + nben = theEnforcedNodes.size(); + + nodeIt = theMeshDS->nodesIterator(); + + // 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 checkedSM; + for (TopExp_Explorer e(theMeshDS->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() ) + nbv -= sm->GetSubMeshDS()->NbNodes(); + } + } + } + + const bool isQuadMesh = + theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) || + theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) || + theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC ); + if ( isQuadMesh ) + { + // descrease nbNodes by nb of medium nodes + while ( nodeIt->more() ) + { + node = nodeIt->next(); + if ( !theHelper.IsDegenShape( node->getshapeId() )) + nbv -= int( theHelper.IsMedium( node )); + } + nodeIt = theMeshDS->nodesIterator(); + } + + std::vector > VerTab; + VerTab.clear(); + std::vector aVerTab; + // Loop from 1 to NB_NODES + + while ( nodeIt->more() ) + { + node = nodeIt->next(); + if (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238 + theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674 + continue; + + aVerTab.clear(); + aVerTab.push_back(node->X()); + aVerTab.push_back(node->Y()); + aVerTab.push_back(node->Z()); + VerTab.push_back(aVerTab); + aGhs3dID++; + theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID )); + theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node )); + } + + /* ENFORCED NODES ========================== */ + if (nben) { + for(enfNodeIt = theEnforcedNodes.begin() ; enfNodeIt != theEnforcedNodes.end() ; ++enfNodeIt) { + double x = (*enfNodeIt)->X(); + double y = (*enfNodeIt)->Y(); + double z = (*enfNodeIt)->Z(); + // Test if point is inside shape to mesh + gp_Pnt myPoint(x,y,z); + BRepClass3d_SolidClassifier scl(theMeshDS->ShapeToMesh()); + scl.Perform(myPoint, 1e-7); + TopAbs_State result = scl.State(); + if ( result != TopAbs_IN ) + continue; + std::vector coords; + coords.push_back(x); + coords.push_back(y); + coords.push_back(z); + if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) + continue; + aVerTab.clear(); + aVerTab.push_back(x); + aVerTab.push_back(y); + aVerTab.push_back(z); + VerTab.push_back(aVerTab); + aGhs3dID++; + theNodeId2NodeIndexMap.insert( make_pair( (*enfNodeIt)->GetID(), aGhs3dID )); + } + } + + /* Write vertices number */ + MESSAGE("Number of vertices: "< > ReqVerTab; + ReqVerTab.clear(); + std::vector aReqVerTab; + int solSize = 0; + 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(theMeshDS->ShapeToMesh()); + scl.Perform(myPoint, 1e-7); + TopAbs_State result = scl.State(); + if ( result != TopAbs_IN ) + continue; + enfVertexSizes.push_back(vertexIt->second); + aReqVerTab.clear(); + aReqVerTab.push_back(x); + aReqVerTab.push_back(y); + aReqVerTab.push_back(z); + ReqVerTab.push_back(aReqVerTab); + solSize++; + } + + if (solSize) { + int idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + if (!idxRequired) + return false; + int idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + if (!idxSol) + return false; + + int TypTab[] = {GmfSca}; + GmfSetKwd(idxRequired, GmfVertices, solSize); + GmfSetKwd(idxSol, GmfSolAtVertices, solSize, 1, TypTab); + + for (int i=0;i::const_iterator itOnMap; + std::vector > tt, qt,et; + tt.clear(); + qt.clear(); + et.clear(); + std::vector att, aqt, aet; + + TopExp::MapShapes( theMeshDS->ShapeToMesh(), TopAbs_FACE, facesMap ); + + for ( int i = 1; i <= facesMap.Extent(); ++i ) + if (( theSubMesh = theProxyMesh.GetSubMesh( facesMap(i)))) + { + SMDS_ElemIteratorPtr it = theSubMesh->GetElements(); + while (it->more()) + { + const SMDS_MeshElement *elem = it->next(); + if (elem->NbNodes() == 3) + { + trianglesMap.Add(facesMap(i)); + nbTriangles ++; + } + else if (elem->NbNodes() == 4) + { + quadranglesMap.Add(facesMap(i)); + nbQuadrangles ++; + } + } + } + + /* TRIANGLES ========================== */ + if (nbTriangles) { + for ( int i = 1; i <= trianglesMap.Extent(); i++ ) + { + aShape = trianglesMap(i); + theSubMesh = theProxyMesh.GetSubMesh(aShape); + if ( !theSubMesh ) continue; + itOnSubMesh = theSubMesh->GetElements(); + while ( itOnSubMesh->more() ) + { + aFace = itOnSubMesh->next(); + itOnSubFace = aFace->nodesIterator(); + att.clear(); + while ( itOnSubFace->more() ) { + // find GHS3D ID + aSmdsID = itOnSubFace->next()->GetID(); + itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID ); + ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() ); + att.push_back((*itOnMap).second); + } + tt.push_back(att); + } + } + } + + if (theEnforcedTriangles.size()) { + // Iterate over the enforced triangles + for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) { + aFace = (*elemIt); + bool isOK = true; + itOnSubFace = aFace->nodesIterator(); + att.clear(); + while ( itOnSubFace->more() ) { + int aNodeID = itOnSubFace->next()->GetID(); + itOnMap = theNodeId2NodeIndexMap.find(aNodeID); + if (itOnMap != theNodeId2NodeIndexMap.end()) + att.push_back((*itOnMap).second); + else { + isOK = false; + theEnforcedTriangles.erase(elemIt); + break; + } + } + if (isOK) + tt.push_back(att); + } + } + + if (tt.size()) { + GmfSetKwd(idx, GmfTriangles, tt.size()); + for (int i=0;iGetElements(); + while ( itOnSubMesh->more() ) + { + aFace = itOnSubMesh->next(); + itOnSubFace = aFace->nodesIterator(); + aqt.clear(); + while ( itOnSubFace->more() ) { + // find GHS3D ID + aSmdsID = itOnSubFace->next()->GetID(); + itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID ); + ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() ); + aqt.push_back((*itOnMap).second); + } + qt.push_back(aqt); + } + } + } + + if (theEnforcedQuadrangles.size()) { + // Iterate over the enforced triangles + for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) { + aFace = (*elemIt); + bool isOK = true; + itOnSubFace = aFace->nodesIterator(); + aqt.clear(); + while ( itOnSubFace->more() ) { + int aNodeID = itOnSubFace->next()->GetID(); + itOnMap = theNodeId2NodeIndexMap.find(aNodeID); + if (itOnMap != theNodeId2NodeIndexMap.end()) + aqt.push_back((*itOnMap).second); + else { + isOK = false; + theEnforcedQuadrangles.erase(elemIt); + break; + } + } + if (isOK) + qt.push_back(aqt); + } + } + + if (qt.size()) { + GmfSetKwd(idx, GmfQuadrilaterals, qt.size()); + for (int i=0;inodesIterator(); + aet.clear(); + while ( itOnSubFace->more() ) { + int aNodeID = itOnSubFace->next()->GetID(); + itOnMap = theNodeId2NodeIndexMap.find(aNodeID); + if (itOnMap != theNodeId2NodeIndexMap.end()) + aet.push_back((*itOnMap).second); + else { + isOK = false; + theEnforcedEdges.erase(elemIt); + break; + } } + if (isOK) + et.push_back(aet); } } - if ( !geomNormalOK) - return invalidID; - // compare normals - bool isReverse = ( meshNormal * geomNormal ) < 0; - if ( !isReverse ) - return meshDS->ShapeToIndex( solid1 ); + if (et.size()) { + GmfSetKwd(idx, GmfEdges, et.size()); + for (int i=0;iShapeToIndex( solids(2) ); + GmfCloseMesh(idx); + return true; } +//======================================================================= +//function : writeFaces +//purpose : Write Faces in case if generate 3D mesh with geometry +//======================================================================= + +// static bool writeFaces (ofstream & theFile, +// const SMESH_ProxyMesh& theMesh, +// const TopoDS_Shape& theShape, +// const map & theSmdsToGhs3dIdMap, +// const map & theEnforcedNodeIdToGhs3dIdMap, +// TIDSortedElemSet & theEnforcedEdges, +// TIDSortedElemSet & theEnforcedTriangles, +// TIDSortedElemSet & theEnforcedQuadrangles) +// { +// // record structure: +// // +// // NB_ELEMS DUMMY_INT +// // Loop from 1 to NB_ELEMS +// // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT +// +// TopoDS_Shape aShape; +// const SMESHDS_SubMesh* theSubMesh; +// const SMDS_MeshElement* aFace; +// const char* space = " "; +// const int dummyint = 0; +// map::const_iterator itOnMap; +// SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace; +// int nbNodes, aSmdsID; +// +// TIDSortedElemSet::const_iterator elemIt; +// int nbEnforcedEdges = theEnforcedEdges.size(); +// int nbEnforcedTriangles = theEnforcedTriangles.size(); +// int nbEnforcedQuadrangles = theEnforcedQuadrangles.size(); +// // count triangles bound to geometry +// int nbTriangles = 0; +// +// TopTools_IndexedMapOfShape facesMap, trianglesMap, quadranglesMap; +// TopExp::MapShapes( theShape, TopAbs_FACE, facesMap ); +// +// 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 and" << std::endl; +// if (nbEnforcedEdges+nbEnforcedTriangles+nbEnforcedQuadrangles) +// std::cout << " " << nbEnforcedEdges+nbEnforcedTriangles+nbEnforcedQuadrangles +// << " enforced shapes:" << std::endl; +// if (nbEnforcedEdges) +// std::cout << " " << nbEnforcedEdges << " enforced edges" << std::endl; +// if (nbEnforcedTriangles) +// std::cout << " " << nbEnforcedTriangles << " enforced triangles" << std::endl; +// if (nbEnforcedQuadrangles) +// std::cout << " " << nbEnforcedQuadrangles << " enforced quadrangles" << std::endl; +// std::cout << std::endl; +// +// // theFile << space << nbTriangles << space << dummyint << std::endl; +// std::ostringstream globalStream, localStream, aStream; +// +// // +// // FACES : BEGIN +// // +// +// for ( int i = 1; i <= facesMap.Extent(); i++ ) +// { +// aShape = facesMap(i); +// theSubMesh = theMesh.GetSubMesh(aShape); +// if ( !theSubMesh ) continue; +// itOnSubMesh = theSubMesh->GetElements(); +// while ( itOnSubMesh->more() ) +// { +// aFace = itOnSubMesh->next(); +// nbNodes = aFace->NbNodes(); +// +// localStream << nbNodes << space; +// +// 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() ); +// +// localStream << (*itOnMap).second << space ; +// } +// +// // (NB_NODES + 1) times: DUMMY_INT +// for ( int j=0; j<=nbNodes; j++) +// localStream << dummyint << space ; +// +// localStream << std::endl; +// } +// } +// +// globalStream << localStream.str(); +// localStream.str(""); +// +// // +// // FACES : END +// // +// +// // +// // ENFORCED EDGES : BEGIN +// // +// +// // Iterate over the enforced edges +// int usedEnforcedEdges = 0; +// bool isOK; +// for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) { +// aFace = (*elemIt); +// isOK = true; +// itOnSubFace = aFace->nodesIterator(); +// aStream.str(""); +// aStream << "2" << space ; +// while ( itOnSubFace->more() ) { +// aSmdsID = itOnSubFace->next()->GetID(); +// itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID); +// if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end()) +// aStream << (*itOnMap).second << space; +// else { +// isOK = false; +// break; +// } +// } +// if (isOK) { +// for ( int j=0; j<=2; j++) +// aStream << dummyint << space ; +// // aStream << dummyint << space << dummyint; +// localStream << aStream.str() << std::endl; +// usedEnforcedEdges++; +// } +// } +// +// if (usedEnforcedEdges) { +// globalStream << localStream.str(); +// localStream.str(""); +// } +// +// // +// // ENFORCED EDGES : END +// // +// // +// +// // +// // ENFORCED TRIANGLES : BEGIN +// // +// // Iterate over the enforced triangles +// int usedEnforcedTriangles = 0; +// for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) { +// aFace = (*elemIt); +// isOK = true; +// itOnSubFace = aFace->nodesIterator(); +// aStream.str(""); +// aStream << "3" << space ; +// while ( itOnSubFace->more() ) { +// aSmdsID = itOnSubFace->next()->GetID(); +// itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID); +// if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end()) +// aStream << (*itOnMap).second << space; +// else { +// isOK = false; +// break; +// } +// } +// if (isOK) { +// for ( int j=0; j<=3; j++) +// aStream << dummyint << space ; +// localStream << aStream.str() << std::endl; +// usedEnforcedTriangles++; +// } +// } +// +// if (usedEnforcedTriangles) { +// globalStream << localStream.str(); +// localStream.str(""); +// } +// +// // +// // ENFORCED TRIANGLES : END +// // +// +// // +// // ENFORCED QUADRANGLES : BEGIN +// // +// // Iterate over the enforced quadrangles +// int usedEnforcedQuadrangles = 0; +// for(elemIt = theEnforcedQuadrangles.begin() ; elemIt != theEnforcedQuadrangles.end() ; ++elemIt) { +// aFace = (*elemIt); +// isOK = true; +// itOnSubFace = aFace->nodesIterator(); +// aStream.str(""); +// aStream << "4" << space ; +// while ( itOnSubFace->more() ) { +// aSmdsID = itOnSubFace->next()->GetID(); +// itOnMap = theEnforcedNodeIdToGhs3dIdMap.find(aSmdsID); +// if (itOnMap != theEnforcedNodeIdToGhs3dIdMap.end()) +// aStream << (*itOnMap).second << space; +// else { +// isOK = false; +// break; +// } +// } +// if (isOK) { +// for ( int j=0; j<=4; j++) +// aStream << dummyint << space ; +// localStream << aStream.str() << std::endl; +// usedEnforcedQuadrangles++; +// } +// } +// +// if (usedEnforcedQuadrangles) { +// globalStream << localStream.str(); +// localStream.str(""); +// } +// // +// // ENFORCED QUADRANGLES : END +// // +// +// theFile +// << nbTriangles+usedEnforcedQuadrangles+usedEnforcedTriangles+usedEnforcedEdges +// << " 0" << std::endl +// << globalStream.str(); +// +// return true; +// } + +//======================================================================= +//function : writePoints +//purpose : +//======================================================================= + +// static bool writePoints (ofstream & theFile, +// SMESH_MesherHelper& theHelper, +// map & theSmdsToGhs3dIdMap, +// map & theEnforcedNodeIdToGhs3dIdMap, +// map & theGhs3dIdToNodeMap, +// GHS3DPlugin_Hypothesis::TID2SizeMap & theNodeIDToSizeMap, +// GHS3DPlugin_Hypothesis::TEnforcedVertexValues & theEnforcedVertices, +// TIDSortedNodeSet & theEnforcedNodes) +// { +// // record structure: +// // +// // NB_NODES +// // 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(); +// int nbEnforcedNodes = theEnforcedNodes.size(); +// +// int aGhs3dID = 1; +// SMDS_NodeIteratorPtr it = theMesh->nodesIterator(); +// const SMDS_MeshNode* node; +// +// // 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 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 bool isQuadMesh = +// theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) || +// theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) || +// theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC ); +// if ( isQuadMesh ) +// { +// // descrease nbNodes by nb of medium nodes +// while ( it->more() ) +// { +// node = it->next(); +// if ( !theHelper.IsDegenShape( node->getshapeId() )) +// nbNodes -= int( theHelper.IsMedium( node )); +// } +// it = theMesh->nodesIterator(); +// } +// +// const char* space = " "; +// const int dummyint = 0; +// +// // NB_NODES +// 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; +// if (nbEnforcedNodes > 0) +// std::cout << " " << nbEnforcedNodes << " enforced nodes" << std::endl; +// +// // std::cout << std::endl; +// // std::cout << "Start writing in 'points' file ..." << std::endl; +// theFile << nbNodes << space << std::endl; +// +// // Loop from 1 to NB_NODES +// +// while ( it->more() ) +// { +// node = it->next(); +// if (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238 +// theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674 +// continue; +// +// theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID )); +// theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node )); +// aGhs3dID++; +// +// // X Y Z DUMMY_INT +// theFile +// << node->X() << space +// << node->Y() << space +// << node->Z() << space +// << dummyint << space ; +// theFile << std::endl; +// +// } +// +// // Iterate over the enforced nodes +// TIDSortedNodeSet::const_iterator nodeIt; +// std::map enfVertexIndexSizeMap; +// if (nbEnforcedNodes) { +// for(nodeIt = theEnforcedNodes.begin() ; nodeIt != theEnforcedNodes.end() ; ++nodeIt) { +// double x = (*nodeIt)->X(); +// double y = (*nodeIt)->Y(); +// double z = (*nodeIt)->Z(); +// // Test if point is inside shape to mesh +// gp_Pnt myPoint(x,y,z); +// BRepClass3d_SolidClassifier scl(theMesh->ShapeToMesh()); +// scl.Perform(myPoint, 1e-7); +// TopAbs_State result = scl.State(); +// if ( result != TopAbs_IN ) +// continue; +// std::vector coords; +// coords.push_back(x); +// coords.push_back(y); +// coords.push_back(z); +// if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) +// continue; +// +// double size = theNodeIDToSizeMap.find((*nodeIt)->GetID())->second; +// // theGhs3dIdToNodeMap.insert( make_pair( nbNodes + i, (*nodeIt) )); +// // MESSAGE("Adding enforced node (" << x << "," << y <<"," << z << ")"); +// // X Y Z PHY_SIZE DUMMY_INT +// theFile +// << x << space +// << y << space +// << z << space +// << size << space +// << dummyint << space; +// theFile << std::endl; +// theEnforcedNodeIdToGhs3dIdMap.insert( make_pair( (*nodeIt)->GetID(), aGhs3dID )); +// enfVertexIndexSizeMap[aGhs3dID] = -1; +// aGhs3dID++; +// // else +// // MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added "); +// } +// } +// +// if (nbEnforcedVertices) { +// // Iterate over the enforced vertices +// GHS3DPlugin_Hypothesis::TEnforcedVertexValues::const_iterator vertexIt; +// // int i = 1; +// 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(theMesh->ShapeToMesh()); +// scl.Perform(myPoint, 1e-7); +// TopAbs_State result = scl.State(); +// if ( result != TopAbs_IN ) +// continue; +// // MESSAGE("Adding enforced vertex (" << x << "," << y <<"," << z << ") = " << vertexIt->second); +// // X Y Z PHY_SIZE DUMMY_INT +// theFile +// << x << space +// << y << space +// << z << space +// << vertexIt->second << space +// << dummyint << space; +// theFile << std::endl; +// enfVertexIndexSizeMap[aGhs3dID] = vertexIt->second; +// aGhs3dID++; +// +// // else +// // MESSAGE("Enforced vertex (" << x << "," << y <<"," << z << ") is not inside the geometry: it was not added "); +// } +// } +// theFile << std::endl; +// +// +// // std::cout << std::endl; +// // std::cout << "End writing in 'points' file." << std::endl; +// +// return true; +// } + + //======================================================================= //function : readResultFile //purpose : @@ -830,12 +2446,17 @@ static bool readResultFile(const int fileOpen, GHS3DPlugin_GHS3D* theAlgo, #endif SMESH_MesherHelper& theHelper, +// SMESH_Mesh& theMesh, TopoDS_Shape tabShape[], double** tabBox, const int nbShape, map & theGhs3dIdToNodeMap, bool toMeshHoles, - int nbEnforcedVertices) + int nbEnforcedVertices, + int nbEnforcedNodes, + TIDSortedElemSet & theEnforcedEdges, + TIDSortedElemSet & theEnforcedTriangles, + TIDSortedElemSet & theEnforcedQuadrangles) { MESSAGE("GHS3DPlugin_GHS3D::readResultFile()"); Kernel_Utils::Localizer loc; @@ -917,7 +2538,7 @@ static bool readResultFile(const int fileOpen, for (int iCoor=0; iCoor < 3; iCoor++) coord[ iCoor ] = strtod(ptr, &ptr); nodeAssigne[ iNode ] = 1; - if ( iNode > (nbInputNodes-nbEnforcedVertices) ) { + if ( iNode > (nbInputNodes-(nbEnforcedVertices+nbEnforcedNodes)) ) { // Creating SMESH nodes // - for enforced vertices // - for vertices of forced edges @@ -950,6 +2571,7 @@ static bool readResultFile(const int fileOpen, const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ]; try { OCC_CATCH_SIGNALS; +// tabID[i] = findShapeID( theHelper, n1, n2, n3, toMeshHoles ); tabID[i] = findShapeID( *theHelper.GetMesh(), n1, n2, n3, toMeshHoles ); // -- 0020330: Pb with ghs3d as a submesh // check that found shape is to be meshed @@ -1094,174 +2716,10 @@ static bool readResultFile(const int fileOpen, return true; } -//======================================================================= -//function : readResultFile -//purpose : -//======================================================================= - -static bool readResultFile(const int fileOpen, -#ifdef WNT - const char* fileName, -#endif -#ifdef WITH_SMESH_CANCEL_COMPUTE - GHS3DPlugin_GHS3D* theAlgo, -#endif - SMESH_MesherHelper& theHelper, - TopoDS_Shape aSolid, - vector & theNodeByGhs3dId, - int nbEnforcedVertices) -{ - SMESHDS_Mesh* theMeshDS = theHelper.GetMeshDS(); - - Kernel_Utils::Localizer loc; - struct stat status; - size_t length; - - char *ptr, *mapPtr; - char *tetraPtr; - char *shapePtr; - - int fileStat; - int nbElems, nbNodes, nbInputNodes; - int nodeId, triangleId; - int nbTriangle; - int ID, shapeID; - - int *tab; - double *coord; - const SMDS_MeshNode **node; - - tab = new int[3]; - coord = new double[3]; - node = new const SMDS_MeshNode*[4]; - - SMDS_MeshNode * aNewNode; - map ::iterator IdNode; - SMDS_MeshElement* aTet; - - // Read the file state - fileStat = fstat(fileOpen, &status); - 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); - tetraPtr = ptr; - - nbElems = tab[0]; - nbNodes = tab[1]; - nbInputNodes = tab[2]; - - theNodeByGhs3dId.resize( nbNodes ); - - // Reading the nodeId - 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 ( theHelper.GetMesh()->NbVolumes() > 0 ) - elemSearcher = SMESH_MeshEditor( theHelper.GetMesh() ).GetElementSearcher(); - - // Reading the nodeCoord and update the nodeMap - shapeID = theMeshDS->ShapeToIndex( aSolid ); - for (int iNode=0; iNode < nbNodes; iNode++) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if(theAlgo->computeCanceled()) - return false; -#endif - for (int iCoor=0; iCoor < 3; iCoor++) - coord[ iCoor ] = strtod(ptr, &ptr); - 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; - } - } - } - - // Reading the triangles - nbTriangle = strtol(ptr, &ptr, 10); - - for (int i=0; i < 3*nbTriangle; i++) - triangleId = strtol(ptr, &ptr, 10); - - shapePtr = ptr; - - if ( theHelper.IsQuadraticMesh() != SMESH_MesherHelper::LINEAR ) - theHelper.IsQuadraticSubMesh( aSolid ); - - // Associating the tetrahedrons to the shapes - - for (int iElem = 0; iElem < nbElems; iElem++) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if(theAlgo->computeCanceled()) - return false; -#endif - for (int iNode = 0; iNode < 4; iNode++) { - 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_TNodeXYZ(node[0]) + - SMESH_TNodeXYZ(node[1]) + - SMESH_TNodeXYZ(node[2]) + - SMESH_TNodeXYZ(node[3]) ) / 4., - SMDSAbs_Volume, foundVolumes )) - continue; - } - aTet = theHelper.AddVolume( node[1], node[0], node[2], node[3], - /*id=*/0, /*force3d=*/true); - 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; - delete [] coord; - delete [] node; - - return true; -} //============================================================================= /*! - *Here we are going to use the GHS3D mesher + *Here we are going to use the GHS3D mesher with geometry */ //============================================================================= @@ -1305,39 +2763,66 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, TCollection_AsciiString aGenericName = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str(); - TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName; - TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName; - aFacesFileName = aGenericName + ".faces"; // in faces - aPointsFileName = aGenericName + ".points"; // in points - aResultFileName = aGenericName + ".noboite";// out points and volumes - aBadResFileName = aGenericName + ".boite"; // out bad result - aBbResFileName = aGenericName + ".bb"; // out vertex stepsize - aLogFileName = aGenericName + ".log"; // log - - // ----------------- - // make input files - // ----------------- - - ofstream aFacesFile ( aFacesFileName.ToCString() , ios::out); - ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out); + TCollection_AsciiString aResultFileName; + TCollection_AsciiString aLogFileName = aGenericName + ".log"; // log +// #if GHS3D_VERSION < 42 +// TCollection_AsciiString aFacesFileName, aPointsFileName; +// TCollection_AsciiString aBadResFileName, aBbResFileName; +// aFacesFileName = aGenericName + ".faces"; // in faces +// aPointsFileName = aGenericName + ".points"; // in points +// aResultFileName = aGenericName + ".noboite";// out points and volumes +// aBadResFileName = aGenericName + ".boite"; // out bad result +// aBbResFileName = aGenericName + ".bb"; // out vertex stepsize +// +// // ----------------- +// // make input files +// // ----------------- +// +// ofstream aFacesFile ( aFacesFileName.ToCString() , ios::out); +// ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out); +// +// Ok = +// aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open(); +// if (!Ok) { +// INFOS( "Can't write into " << aFacesFileName); +// return error(SMESH_Comment("Can't write into ") << aFacesFileName); +// } +// #else + TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName; + TCollection_AsciiString aResultGMFFileName; - Ok = - aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open(); - if (!Ok) { - INFOS( "Can't write into " << aFacesFileName); - return error(SMESH_Comment("Can't write into ") << aFacesFileName); - } - map aSmdsToGhs3dIdMap; +#ifdef _DEBUG_ + aGMFFileName = aGenericName + ".mesh"; // GMF mesh file + // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) + aResultGMFFileName = aGenericName + "Vol.mesh"; // GMF mesh file + aResultFileName = aGenericName + ".noboite";// out points and volumes +// aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file + aRequiredVerticesFileName = aGenericName + "_required.mesh"; // GMF required vertices mesh file + aSolFileName = aGenericName + "_required.sol"; // GMF solution file +#else + aGMFFileName = aGenericName + ".mesh"; // GMF mesh file +// aGMFFileName = aGenericName + ".meshb"; // GMF mesh file + // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) + aResultGMFFileName = aGenericName + "Vol.meshb"; // GMF mesh file + aResultFileName = aGenericName + ".noboite";// out points and volumes +// aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file + aRequiredVerticesFileName = aGenericName + "_required.meshb"; // GMF required vertices mesh file + aSolFileName = aGenericName + "_required.solb"; // GMF solution file +#endif + map aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap; map aGhs3dIdToNodeMap; - map,double> enforcedVertices; - int nbEnforcedVertices = 0; - try { - enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp); - nbEnforcedVertices = enforcedVertices.size(); - } - catch(...) { - } - + std::map nodeID2nodeIndexMap; + GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp); + TIDSortedNodeSet enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp); + TIDSortedElemSet enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp); + TIDSortedElemSet enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp); + TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp); + GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp); +// GHS3DPlugin_Hypothesis::TID2SizeMap elemIDToSizeMap = GHS3DPlugin_Hypothesis::GetElementIDToSizeMap(_hyp); + + int nbEnforcedVertices = enforcedVertices.size(); + int nbEnforcedNodes = enforcedNodes.size(); + SMESH_MesherHelper helper( theMesh ); helper.SetSubShape( theShape ); @@ -1369,18 +2854,29 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, if ( !proxyMesh ) return false; } - - Ok = (writePoints( aPointsFile, helper, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, enforcedVertices) - && - writeFaces ( aFacesFile, *proxyMesh, theShape, aSmdsToGhs3dIdMap )); +// #if GHS3D_VERSION < 42 +// Ok = (writePoints( aPointsFile, helper, +// aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, aGhs3dIdToNodeMap, +// nodeIDToSizeMap, +// enforcedVertices, enforcedNodes) +// && +// writeFaces ( aFacesFile, *proxyMesh, theShape, +// aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, +// enforcedEdges, enforcedTriangles, enforcedQuadrangles)); +// #else + Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(), + helper, *proxyMesh, + aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, + enforcedNodes, enforcedEdges, enforcedTriangles, enforcedQuadrangles, + enforcedVertices); +// #endif } // Write aSmdsToGhs3dIdMap to temp file TCollection_AsciiString aSmdsToGhs3dIdMapFileName; aSmdsToGhs3dIdMapFileName = aGenericName + ".ids"; // ids relation ofstream aIdsFile ( aSmdsToGhs3dIdMapFileName.ToCString() , ios::out); - Ok = - aIdsFile.rdbuf()->is_open(); + Ok = aIdsFile.rdbuf()->is_open(); if (!Ok) { INFOS( "Can't write into " << aSmdsToGhs3dIdMapFileName); return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName); @@ -1391,14 +2887,22 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, aIdsFile << myit->first << " " << myit->second << std::endl; } - aFacesFile.close(); - aPointsFile.close(); aIdsFile.close(); +// #if GHS3D_VERSION < 42 +// aFacesFile.close(); +// aPointsFile.close(); +// #endif if ( ! Ok ) { if ( !_keepFiles ) { - removeFile( aFacesFileName ); - removeFile( aPointsFileName ); +// #if GHS3D_VERSION < 42 +// removeFile( aFacesFileName ); +// removeFile( aPointsFileName ); +// #else + removeFile( aGMFFileName ); + removeFile( aRequiredVerticesFileName ); + removeFile( aSolFileName ); +// #endif removeFile( aSmdsToGhs3dIdMapFileName ); } return error(COMPERR_BAD_INPUT_MESH); @@ -1409,9 +2913,14 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, // run ghs3d mesher // ----------------- - TCollection_AsciiString cmd( (char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() ); + TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() ); + // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) cmd += TCollection_AsciiString(" -f ") + aGenericName; // file to read - cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file + cmd += TCollection_AsciiString(" -IM "); +// cmd += TCollection_AsciiString(" --in ") + aGenericName; +// cmd += TCollection_AsciiString(" --required_vertices ") + aRequiredVerticesFileName; +// cmd += TCollection_AsciiString(" --out ") + aGenericName; + cmd += TCollection_AsciiString(" -Om 1>" ) + aLogFileName; // dump into file std::cout << std::endl; std::cout << "Ghs3d execution..." << std::endl; @@ -1430,6 +2939,7 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, // read a result // -------------- +// #if GHS3D_VERSION < 42 // Mapping the result file int fileOpen; @@ -1454,9 +2964,37 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, #ifdef WITH_SMESH_CANCEL_COMPUTE this, #endif - helper, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap, - toMeshHoles, nbEnforcedVertices ); + /*theMesh, */helper, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap, + toMeshHoles, + nbEnforcedVertices, nbEnforcedNodes, + enforcedEdges, enforcedTriangles, enforcedQuadrangles ); } +// /*/*#else +// #ifndef GMF_HAS_SUBDOMAIN_INFO +// // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) +// +// int fileOpen = open( aResultFileName.ToCString(), O_RDONLY); +// if ( fileOpen < 0 ) { +// 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 { +// #endif +// Ok = readGMFFile( +// #ifndef GMF_HAS_SUBDOMAIN_INFO +// fileOpen, +// #endif +// aGenericName.ToCString(), theMesh, +// _nbShape, tabShape, tabBox, +// aGhs3dIdToNodeMap, toMeshHoles, +// nbEnforcedVertices, nbEnforcedNodes, +// enforcedNodes, enforcedTriangles, enforcedQuadrangles); +// #ifndef GMF_HAS_SUBDOMAIN_INFO +// } +// #endif +// #endif*/*/ // --------------------- // remove working files @@ -1482,17 +3020,21 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, } if ( !_keepFiles ) { +// #if GHS3D_VERSION < 42 +// removeFile( aFacesFileName ); +// removeFile( aPointsFileName ); +// removeFile( aResultFileName ); +// removeFile( aBadResFileName ); +// removeFile( aBbResFileName ); +// #endif + removeFile( aSmdsToGhs3dIdMapFileName ); + // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) + #ifdef WITH_SMESH_CANCEL_COMPUTE if (! Ok) if(_compute_canceled) removeFile( aLogFileName ); #endif - removeFile( aFacesFileName ); - removeFile( aPointsFileName ); - removeFile( aResultFileName ); - removeFile( aBadResFileName ); - removeFile( aBbResFileName ); - removeFile( aSmdsToGhs3dIdMapFileName ); } std::cout << "<" << aResultFileName.ToCString() << "> GHS3D output file "; if ( !Ok ) @@ -1525,37 +3067,53 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, TCollection_AsciiString aGenericName = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str(); - TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName; - TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName; - aFacesFileName = aGenericName + ".faces"; // in faces - aPointsFileName = aGenericName + ".points"; // in points - aResultFileName = aGenericName + ".noboite";// out points and volumes - aBadResFileName = aGenericName + ".boite"; // out bad result - aBbResFileName = aGenericName + ".bb"; // out vertex stepsize - aLogFileName = aGenericName + ".log"; // log - - // ----------------- - // make input files - // ----------------- - - ofstream aFacesFile ( aFacesFileName.ToCString() , ios::out); - ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out); - bool Ok = - aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open(); - - if (!Ok) - return error( SMESH_Comment("Can't write into ") << aPointsFileName); + TCollection_AsciiString aLogFileName = aGenericName + ".log"; // log + TCollection_AsciiString aResultFileName; + bool Ok; +// #if GHS3D_VERSION < 42 +// TCollection_AsciiString aFacesFileName, aPointsFileName; +// TCollection_AsciiString aBadResFileName, aBbResFileName; +// aFacesFileName = aGenericName + ".faces"; // in faces +// aPointsFileName = aGenericName + ".points"; // in points +// aResultFileName = aGenericName + ".noboite";// out points and volumes +// aBadResFileName = aGenericName + ".boite"; // out bad result +// aBbResFileName = aGenericName + ".bb"; // out vertex stepsize +// +// // ----------------- +// // make input files +// // ----------------- +// +// ofstream aFacesFile ( aFacesFileName.ToCString() , ios::out); +// ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out); +// Ok = aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open(); +// if (!Ok) { +// INFOS( "Can't write into " << aFacesFileName); +// return error( SMESH_Comment("Can't write into ") << aFacesFileName); +// } +// #else + TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName; +#ifdef _DEBUG_ + aGMFFileName = aGenericName + ".mesh"; // GMF mesh file + aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file + aRequiredVerticesFileName = aGenericName + "_required.mesh"; // GMF required vertices mesh file + aSolFileName = aGenericName + "_required.sol"; // GMF solution file +#else + aGMFFileName = aGenericName + ".meshb"; // GMF mesh file + aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file + aRequiredVerticesFileName = aGenericName + "_required.meshb"; // GMF required vertices mesh file + aSolFileName = aGenericName + ".solb"; // GMF solution file +#endif +// #endif - GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices; - int nbEnforcedVertices = 0; - try { - enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp); - nbEnforcedVertices = enforcedVertices.size(); - } - catch(...) { - } - - vector aNodeByGhs3dId; + std::map nodeID2nodeIndexMap; + GHS3DPlugin_Hypothesis::TEnforcedVertexValues enforcedVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp); + TIDSortedNodeSet enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp); + TIDSortedElemSet enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp); + TIDSortedElemSet enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp); + TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp); + GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp); + + vector aNodeByGhs3dId, anEnforcedNodeByGhs3dId; { SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh )); if ( theMesh.NbQuadrangles() > 0 ) @@ -1564,30 +3122,57 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, aQuad2Trias->Compute( theMesh ); proxyMesh.reset( aQuad2Trias ); } - - Ok = (writeFaces ( aFacesFile, *proxyMesh, aNodeByGhs3dId ) && - writePoints( aPointsFile, &theMesh, aNodeByGhs3dId,enforcedVertices)); - } - aFacesFile.close(); - aPointsFile.close(); - - if ( ! Ok ) { - if ( !_keepFiles ) { - removeFile( aFacesFileName ); - removeFile( aPointsFileName ); - } - return error(COMPERR_BAD_INPUT_MESH); +// #if GHS3D_VERSION < 42 +// Ok = (writeFaces ( aFacesFile, *proxyMesh, &theMesh, aNodeByGhs3dId, anEnforcedNodeByGhs3dId, +// enforcedEdges, enforcedTriangles, enforcedQuadrangles ) && +// writePoints( aPointsFile, &theMesh, aNodeByGhs3dId, anEnforcedNodeByGhs3dId, +// nodeIDToSizeMap, enforcedVertices, enforcedNodes)); +// int nbEnforcedVertices = enforcedVertices.size(); +// int nbEnforcedNodes = enforcedNodes.size(); +// #else + Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(), + *proxyMesh, &theMesh, + aNodeByGhs3dId, anEnforcedNodeByGhs3dId, + enforcedNodes, enforcedEdges, enforcedTriangles, enforcedQuadrangles, + enforcedVertices); +// #endif } - removeFile( aResultFileName ); // needed for boundary recovery module usage + + TIDSortedNodeSet enforcedNodesFromEnforcedElem; + for (int i=0;i" ) + aLogFileName; // dump into file + + std::cout << std::endl; + std::cout << "Ghs3d execution..." << std::endl; + std::cout << cmd << std::endl; #ifdef WITH_SMESH_CANCEL_COMPUTE _compute_canceled = false; @@ -1595,28 +3180,41 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, system( cmd.ToCString() ); // run + 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 ) { - 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, -#ifdef WNT - aResultFileName.ToCString(), -#endif +// #if GHS3D_VERSION < 42 +// int fileOpen = open( aResultFileName.ToCString(), O_RDONLY); +// if ( fileOpen < 0 ) { +// 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, +// #ifdef WNT +// aResultFileName.ToCString(), +// #endif +// #ifdef WITH_SMESH_CANCEL_COMPUTE +// this, +// #endif +// theMesh, theShape ,aNodeByGhs3dId, anEnforcedNodeByGhs3dId, +// nbEnforcedVertices, nbEnforcedNodes, +// enforcedEdges, enforcedTriangles, enforcedQuadrangles ); +// } +// #else + Ok = readGMFFile(aResultFileName.ToCString(), #ifdef WITH_SMESH_CANCEL_COMPUTE - this, + this, #endif - *theHelper, theShape ,aNodeByGhs3dId, nbEnforcedVertices ); - } + theHelper, + enforcedNodesFromEnforcedElem, enforcedTriangles, enforcedQuadrangles); +// #endif // --------------------- // remove working files @@ -1639,7 +3237,7 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" ); error(COMPERR_ALGO_FAILED, "ghs3d: command not found" ); } - +// #if GHS3D_VERSION < 42 if ( !_keepFiles ) { #ifdef WITH_SMESH_CANCEL_COMPUTE @@ -1647,13 +3245,13 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, if(_compute_canceled) removeFile( aLogFileName ); #endif - removeFile( aFacesFileName ); - removeFile( aPointsFileName ); - removeFile( aResultFileName ); - removeFile( aBadResFileName ); - removeFile( aBbResFileName ); +// removeFile( aFacesFileName ); +// removeFile( aPointsFileName ); +// removeFile( aResultFileName ); +// removeFile( aBadResFileName ); +// removeFile( aBbResFileName ); } - +// #endif return Ok; } @@ -2253,3 +3851,14 @@ bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh, return true; } +bool GHS3DPlugin_GHS3D::importGMFMesh(const char* theGMFFileName, SMESH_Mesh& theMesh) +{ + SMESH_MesherHelper* helper = new SMESH_MesherHelper(theMesh ); + TIDSortedElemSet dummyElemSet; + TIDSortedNodeSet dummyNodeSet; + return readGMFFile(theGMFFileName, +#ifdef WITH_SMESH_CANCEL_COMPUTE + this, +#endif + helper, dummyNodeSet , dummyElemSet, dummyElemSet); +} diff --git a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.hxx b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.hxx index 1c2c320..3888942 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.hxx +++ b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.hxx @@ -27,10 +27,21 @@ #define _GHS3DPlugin_GHS3D_HXX_ #include "SMESH_3D_Algo.hxx" +#include "SMESH_Gen.hxx" #include #include +extern "C" +{ + #include "libmesh5.h" +} + +#ifndef GMFVERSION +#define GMFVERSION GmfDouble +#endif +#define GMFDIMENSION 3 + class GHS3DPlugin_Hypothesis; class SMDS_MeshNode; class SMESH_Mesh; @@ -62,6 +73,8 @@ public: virtual bool Compute(SMESH_Mesh& theMesh, SMESH_MesherHelper* aHelper); + bool importGMFMesh(const char* aGMFFileName, SMESH_Mesh& aMesh); + private: bool storeErrorDescription(const TCollection_AsciiString& logFile, diff --git a/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.cxx b/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.cxx index c63525a..cee9677 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.cxx +++ b/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.cxx @@ -25,9 +25,13 @@ // #include "GHS3DPlugin_GHS3D_i.hxx" #include "SMESH_Gen.hxx" +#include "SMESH_Mesh_i.hxx" +#include "SMESH_Gen_i.hxx" #include "GHS3DPlugin_GHS3D.hxx" +#include "SMESH_PythonDump.hxx" #include "utilities.h" +#include using namespace std; @@ -80,3 +84,36 @@ GHS3DPlugin_GHS3D_i::~GHS3DPlugin_GHS3D_i() return ( ::GHS3DPlugin_GHS3D* )myBaseImpl; } +//============================================================================= +/*! + * GHS3DPlugin_GHS3D_i::~GHS3DPlugin_GHS3D_i + * + * Destructor + */ +//============================================================================= + +bool GHS3DPlugin_GHS3D_i::importGMFMesh(const char* theGMFFileName) +{ + MESSAGE( "GHS3DPlugin_GHS3D_i::importGMFMesh" ); + + SMESH::SMESH_Mesh_ptr theMesh = SMESH_Gen_i::GetSMESHGen()->CreateEmptyMesh(); + SMESH_Gen_i::GetSMESHGen()->RemoveLastFromPythonScript(SMESH_Gen_i::GetSMESHGen()->GetCurrentStudy()->StudyId()); + SALOMEDS::SObject_ptr theSMesh = SMESH_Gen_i::GetSMESHGen()->ObjectToSObject(SMESH_Gen_i::GetSMESHGen()->GetCurrentStudy(), theMesh); +#ifdef WINNT +#define SEP '\\' +#else +#define SEP '/' +#endif + string strFileName (theGMFFileName); + strFileName = strFileName.substr(strFileName.rfind(SEP)+1); + strFileName.erase(strFileName.rfind('.')); + SMESH_Gen_i::GetSMESHGen()->SetName(theSMesh, strFileName.c_str()); + SMESH_Mesh_i* meshServant = dynamic_cast( SMESH_Gen_i::GetSMESHGen()->GetServant( theMesh ).in() ); + ASSERT( meshServant ); + if ( meshServant ) { + bool res = GetImpl()->importGMFMesh(theGMFFileName, meshServant->GetImpl()); + SMESH::TPythonDump() << "isDone = " << _this() << ".importGMFMesh( \"" << theGMFFileName << "\")"; + return res; + } + return false; +} diff --git a/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.hxx b/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.hxx index e61fb45..9de897f 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.hxx +++ b/src/GHS3DPlugin/GHS3DPlugin_GHS3D_i.hxx @@ -49,6 +49,8 @@ public: // Get implementation ::GHS3DPlugin_GHS3D* GetImpl(); + + virtual bool importGMFMesh(const char* theGMFFileName); }; #endif diff --git a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.cxx b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.cxx index bdfa2b9..7498235 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.cxx +++ b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.cxx @@ -24,6 +24,8 @@ //============================================================================= // #include "GHS3DPlugin_Hypothesis.hxx" +#include +#include #include @@ -49,10 +51,19 @@ GHS3DPlugin_Hypothesis::GHS3DPlugin_Hypothesis(int hypId, int studyId, SMESH_Gen myToUseBoundaryRecoveryVersion(DefaultToUseBoundaryRecoveryVersion()), myToUseFemCorrection(DefaultToUseFEMCorrection()), myToRemoveCentralPoint(DefaultToRemoveCentralPoint()), - myEnforcedVertices(DefaultEnforcedVertices()) + myEnforcedVertices(DefaultEnforcedVertices()), + _enfNodes(DefaultIDSortedNodeSet()), + _enfEdges(DefaultIDSortedElemSet()), + _enfTriangles(DefaultIDSortedElemSet()), + _enfQuadrangles(DefaultIDSortedElemSet()) { _name = "GHS3D_Parameters"; _param_algo_dim = 3; + _edgeID2nodeIDMap.clear(); + _triID2nodeIDMap.clear(); + _quadID2nodeIDMap.clear(); + _nodeIDToSizeMap.clear(); + _elementIDToSizeMap.clear(); } //======================================================================= @@ -328,6 +339,181 @@ void GHS3DPlugin_Hypothesis::SetEnforcedVertex(double x, double y, double z, dou NotifySubMeshesHypothesisModification(); } + +//======================================================================= +//function : SetEnforcedMesh +//======================================================================= +void GHS3DPlugin_Hypothesis::SetEnforcedMesh(SMESH_Mesh& theMesh, SMESH::ElementType elementType, double size) +{ + TIDSortedElemSet theElemSet; + SMDS_ElemIteratorPtr eIt; +/* + if ((elementType == SMESH::FACE) && (theMesh.NbQuadrangles() > 0)) { + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh )); + + StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor; + aQuad2Trias->Compute( theMesh ); + proxyMesh.reset(aQuad2Trias ); + +// std::cout << "proxyMesh->NbFaces(): " << proxyMesh->NbFaces() << std::endl; +// eIt = proxyMesh->GetFaces(); +// if (eIt) +// while ( eIt->more() ) +// theElemSet.insert( eIt->next() ); +// else { +// std::cout << "********************** eIt == 0 *****************" << std::endl; + eIt = theMesh.GetMeshDS()->elementsIterator(SMDSAbs_ElementType(elementType)); + while ( eIt->more() ) { + const SMDS_MeshElement* elem = eIt->next(); + theElemSet.insert( elem ); + } + } + + else + { + */ + eIt = theMesh.GetMeshDS()->elementsIterator(SMDSAbs_ElementType(elementType)); + while ( eIt->more() ) + theElemSet.insert( eIt->next() ); +/* + } +*/ + MESSAGE("Add "<GetID()); + const SMDS_MeshNode* node = dynamic_cast(elem); + switch (elementType) { + case SMESH::NODE: + if (node) { +// MESSAGE("This is a node"); + _enfNodes.insert(node); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + else { +// MESSAGE("This is an element"); + _enfNodes.insert(elem->begin_nodes(),elem->end_nodes()); + SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); + for (;nodeIt->more();) + _nodeIDToSizeMap.insert(make_pair(nodeIt->next()->GetID(), size)); + } + NotifySubMeshesHypothesisModification(); + break; + case SMESH::EDGE: + if (node) { +// MESSAGE("This is a node"); + } + else { +// MESSAGE("This is an element"); + if (elem->GetType() == SMDSAbs_Edge) { + _enfEdges.insert(elem); +// _enfNodes.insert(elem->begin_nodes(),elem->end_nodes()); + _elementIDToSizeMap.insert(make_pair(elem->GetID(), size)); + SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); + for (;nodeIt->more();) { + node = dynamic_cast(nodeIt->next()); + _edgeID2nodeIDMap[elem->GetID()].push_back(node->GetID()); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + } + else if (elem->GetType() > SMDSAbs_Edge) { + SMDS_ElemIteratorPtr it = elem->edgesIterator(); + for (;it->more();) { + const SMDS_MeshElement* anEdge = it->next(); + _enfEdges.insert(anEdge); +// _enfNodes.insert(anEdge->begin_nodes(),anEdge->end_nodes()); + _elementIDToSizeMap.insert(make_pair(anEdge->GetID(), size)); + SMDS_ElemIteratorPtr nodeIt = anEdge->nodesIterator(); + for (;nodeIt->more();) { + node = dynamic_cast(nodeIt->next()); + _edgeID2nodeIDMap[anEdge->GetID()].push_back(node->GetID()); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + } + } + } + NotifySubMeshesHypothesisModification(); + break; + case SMESH::FACE: + if (node) { +// MESSAGE("This is a node"); + } + else { +// MESSAGE("This is an element"); + if (elem->GetType() == SMDSAbs_Face) + { + if (elem->NbNodes() == 3) { + _enfTriangles.insert(elem); +// _enfNodes.insert(elem->begin_nodes(),elem->end_nodes()); + _elementIDToSizeMap.insert(make_pair(elem->GetID(), size)); + SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); + for (;nodeIt->more();) { + node = dynamic_cast(nodeIt->next()); + _triID2nodeIDMap[elem->GetID()].push_back(node->GetID()); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + } + else if (elem->NbNodes() == 4) { + _enfQuadrangles.insert(elem); +// _enfNodes.insert(elem->begin_nodes(),elem->end_nodes()); + _elementIDToSizeMap.insert(make_pair(elem->GetID(), size)); + SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); + for (;nodeIt->more();) { + node = dynamic_cast(nodeIt->next()); + _quadID2nodeIDMap[elem->GetID()].push_back(node->GetID()); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + } + } + else if (elem->GetType() > SMDSAbs_Face) { + SMDS_ElemIteratorPtr it = elem->facesIterator(); + for (;it->more();) { + const SMDS_MeshElement* aFace = it->next(); + if (aFace->NbNodes() == 3) { + _enfTriangles.insert(aFace); +// _enfNodes.insert(aFace->begin_nodes(),aFace->end_nodes()); + _elementIDToSizeMap.insert(make_pair(aFace->GetID(), size)); + SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator(); + for (;nodeIt->more();) { + node = dynamic_cast(nodeIt->next()); + _triID2nodeIDMap[aFace->GetID()].push_back(node->GetID()); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + } + else if (aFace->NbNodes() == 4) { + _enfQuadrangles.insert(aFace); +// _enfNodes.insert(aFace->begin_nodes(),aFace->end_nodes()); + _elementIDToSizeMap.insert(make_pair(aFace->GetID(), size)); + SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator(); + for (;nodeIt->more();) { + node = dynamic_cast(nodeIt->next()); + _quadID2nodeIDMap[aFace->GetID()].push_back(node->GetID()); + _nodeIDToSizeMap.insert(make_pair(node->GetID(), size)); + } + } + } + } + } + NotifySubMeshesHypothesisModification(); + break; + }; + } +} + //======================================================================= //function : GetEnforcedVertex //======================================================================= @@ -377,6 +563,24 @@ void GHS3DPlugin_Hypothesis::ClearEnforcedVertices() NotifySubMeshesHypothesisModification(); } +//======================================================================= +//function : ClearEnforcedMeshes +//======================================================================= +void GHS3DPlugin_Hypothesis::ClearEnforcedMeshes() +{ + _enfNodes.clear(); + _enfEdges.clear(); + _enfTriangles.clear(); + _enfQuadrangles.clear(); +// _edgeID2nodeIDMap.clear(); +// _triID2nodeIDMap.clear(); +// _quadID2nodeIDMap.clear(); +// _nodeIDToSizeMap.clear(); +// _elementIDToSizeMap.clear(); + NotifySubMeshesHypothesisModification(); +} + + //======================================================================= //function : DefaultMeshHoles //======================================================================= @@ -524,6 +728,15 @@ GHS3DPlugin_Hypothesis::TEnforcedVertexValues GHS3DPlugin_Hypothesis::DefaultEnf return GHS3DPlugin_Hypothesis::TEnforcedVertexValues(); } +TIDSortedNodeSet GHS3DPlugin_Hypothesis::DefaultIDSortedNodeSet() +{ + return TIDSortedNodeSet(); +} + +TIDSortedElemSet GHS3DPlugin_Hypothesis::DefaultIDSortedElemSet() +{ + return TIDSortedElemSet(); +} //======================================================================= //function : SaveTo @@ -744,7 +957,7 @@ bool GHS3DPlugin_Hypothesis::SetParametersByDefaults(const TDefaults& /*dflts*/ std::string GHS3DPlugin_Hypothesis::CommandToRun(const GHS3DPlugin_Hypothesis* hyp, const bool hasShapeToMesh) { - TCollection_AsciiString cmd( "ghs3d" ); + TCollection_AsciiString cmd( "ghs3d" ); // ghs3dV4.1_32bits (to permit the .mesh output) // check if any option is overridden by hyp->myTextOption bool m = hyp ? ( hyp->myTextOption.find("-m") == std::string::npos ) : true; bool M = hyp ? ( hyp->myTextOption.find("-M") == std::string::npos ) : true; @@ -881,3 +1094,47 @@ GHS3DPlugin_Hypothesis::TEnforcedVertexValues GHS3DPlugin_Hypothesis::GetEnforce return hyp ? hyp->_GetEnforcedVertices():DefaultEnforcedVertices(); } +TIDSortedNodeSet GHS3DPlugin_Hypothesis::GetEnforcedNodes(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetEnforcedNodes():DefaultIDSortedNodeSet(); +} + +TIDSortedElemSet GHS3DPlugin_Hypothesis::GetEnforcedEdges(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetEnforcedEdges():DefaultIDSortedElemSet(); +} + +TIDSortedElemSet GHS3DPlugin_Hypothesis::GetEnforcedTriangles(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetEnforcedTriangles():DefaultIDSortedElemSet(); +} + +TIDSortedElemSet GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetEnforcedQuadrangles():DefaultIDSortedElemSet(); +} + +GHS3DPlugin_Hypothesis::TElemID2NodeIDMap GHS3DPlugin_Hypothesis::GetEdgeID2NodeIDMap(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetEdgeID2NodeIDMap(): GHS3DPlugin_Hypothesis::TElemID2NodeIDMap(); +} + +GHS3DPlugin_Hypothesis::TElemID2NodeIDMap GHS3DPlugin_Hypothesis::GetTri2NodeMap(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetTri2NodeMap(): GHS3DPlugin_Hypothesis::TElemID2NodeIDMap(); +} + +GHS3DPlugin_Hypothesis::TElemID2NodeIDMap GHS3DPlugin_Hypothesis::GetQuad2NodeMap(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetQuad2NodeMap(): GHS3DPlugin_Hypothesis::TElemID2NodeIDMap(); +} + +GHS3DPlugin_Hypothesis::TID2SizeMap GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetNodeIDToSizeMap(): GHS3DPlugin_Hypothesis::TID2SizeMap(); +} + +GHS3DPlugin_Hypothesis::TID2SizeMap GHS3DPlugin_Hypothesis::GetElementIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp) +{ + return hyp ? hyp->_GetElementIDToSizeMap(): GHS3DPlugin_Hypothesis::TID2SizeMap(); +} diff --git a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.hxx b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.hxx index 28fe2bb..ac54676 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.hxx +++ b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis.hxx @@ -27,8 +27,11 @@ #include "GHS3DPlugin_Defs.hxx" -#include -#include +#include "SMESH_Hypothesis.hxx" +#include "SMESH_Mesh_i.hxx" +#include "SMESH_Gen_i.hxx" +#include "SMESH_TypeDefs.hxx" +#include "utilities.h" #include #include @@ -117,6 +120,39 @@ public: /*! * To set an enforced vertex */ +// struct TEnforcedNode { +// std::vector coords; +// double size; +// std::string geomEntry; +// std::string groupName; +// }; +// +// struct CompareEnfNodes { +// bool operator () (const TEnforcedNode* e1, const TEnforcedNode* e2) const { +// if (e1 && e2) { +// if (e1->coords.size() && e2->coords.size()) +// return (e1->coords < e2->coords); +// else +// return (e1->geomEntry < e2->geomEntry); +// } +// return false; +// } +// }; +// typedef std::set< TEnforcedNode*, CompareEnfNodes > TEnforcedNodeList; +// // Map Coords / Enforced node +// typedef std::map< std::vector, TEnforcedNode* > TCoordsEnfNodeMap; +// // Map geom entry / Enforced ndoe +// typedef std::map< std::string, TEnforcedNode* > TGeomEntryEnfNodeMap; +// +// +// struct TEnforcedEdge { +// long ID; +// long node1; +// long node2; +// std::string groupName; +// }; + + typedef std::map,double> TEnforcedVertexValues; void SetEnforcedVertex(double x, double y, double z, double size); double GetEnforcedVertex(double x, double y, double z) throw (std::invalid_argument); @@ -124,6 +160,24 @@ public: const TEnforcedVertexValues _GetEnforcedVertices() const { return myEnforcedVertices; } void ClearEnforcedVertices(); + /*! + * To set enforced elements + */ + void SetEnforcedMesh(SMESH_Mesh& theMesh, SMESH::ElementType elementType, double size); + void SetEnforcedElements(TIDSortedElemSet theElemSet, SMESH::ElementType elementType, double size); + void ClearEnforcedMeshes(); + const TIDSortedNodeSet _GetEnforcedNodes() const { return _enfNodes; } + const TIDSortedElemSet _GetEnforcedEdges() const { return _enfEdges; } + const TIDSortedElemSet _GetEnforcedTriangles() const { return _enfTriangles; } + const TIDSortedElemSet _GetEnforcedQuadrangles() const { return _enfQuadrangles; } + typedef std::map > TElemID2NodeIDMap; + const TElemID2NodeIDMap _GetEdgeID2NodeIDMap() const {return _edgeID2nodeIDMap; } + const TElemID2NodeIDMap _GetTri2NodeMap() const {return _triID2nodeIDMap; } + const TElemID2NodeIDMap _GetQuad2NodeMap() const {return _quadID2nodeIDMap; } + typedef std::map TID2SizeMap; + const TID2SizeMap _GetNodeIDToSizeMap() const {return _nodeIDToSizeMap; } + const TID2SizeMap _GetElementIDToSizeMap() const {return _elementIDToSizeMap; } + static bool DefaultMeshHoles(); static short DefaultMaximumMemory(); static short DefaultInitialMemory(); @@ -136,6 +190,8 @@ public: static bool DefaultToUseFEMCorrection(); static bool DefaultToRemoveCentralPoint(); static TEnforcedVertexValues DefaultEnforcedVertices(); + static TIDSortedNodeSet DefaultIDSortedNodeSet(); + static TIDSortedElemSet DefaultIDSortedElemSet(); /*! * \brief Return command to run ghs3d mesher excluding file prefix (-f) @@ -150,7 +206,16 @@ public: * \brief Return the enforced vertices */ static TEnforcedVertexValues GetEnforcedVertices(const GHS3DPlugin_Hypothesis* hyp); - + static TIDSortedNodeSet GetEnforcedNodes(const GHS3DPlugin_Hypothesis* hyp); + static TIDSortedElemSet GetEnforcedEdges(const GHS3DPlugin_Hypothesis* hyp); + static TIDSortedElemSet GetEnforcedTriangles(const GHS3DPlugin_Hypothesis* hyp); + static TIDSortedElemSet GetEnforcedQuadrangles(const GHS3DPlugin_Hypothesis* hyp); + static TElemID2NodeIDMap GetEdgeID2NodeIDMap(const GHS3DPlugin_Hypothesis* hyp); + static TElemID2NodeIDMap GetTri2NodeMap(const GHS3DPlugin_Hypothesis* hyp); + static TElemID2NodeIDMap GetQuad2NodeMap(const GHS3DPlugin_Hypothesis* hyp); + static TID2SizeMap GetNodeIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp); + static TID2SizeMap GetElementIDToSizeMap(const GHS3DPlugin_Hypothesis* hyp); + // Persistence virtual std::ostream & SaveTo(std::ostream & save); virtual std::istream & LoadFrom(std::istream & load); @@ -182,7 +247,15 @@ private: bool myToRemoveCentralPoint; std::string myTextOption; TEnforcedVertexValues myEnforcedVertices; - + TIDSortedNodeSet _enfNodes; + TIDSortedElemSet _enfEdges; + TIDSortedElemSet _enfTriangles; + TIDSortedElemSet _enfQuadrangles; + TElemID2NodeIDMap _edgeID2nodeIDMap; + TElemID2NodeIDMap _triID2nodeIDMap; + TElemID2NodeIDMap _quadID2nodeIDMap; + TID2SizeMap _nodeIDToSizeMap; + TID2SizeMap _elementIDToSizeMap; }; diff --git a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.cxx b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.cxx index ff3953e..63952ba 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.cxx +++ b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.cxx @@ -23,13 +23,22 @@ // #include "GHS3DPlugin_Hypothesis_i.hxx" -#include -#include - -#include -#include -#include - +#include "SMESH_Gen.hxx" +#include "SMESH_PythonDump.hxx" +//#include "SMESH_Mesh.hxx" +//#include "SMESH_ProxyMesh.hxx" +//#include + +#include "Utils_CorbaException.hxx" +#include "utilities.h" +#include "SMESH_Mesh_i.hxx" +#include "SMESH_Group_i.hxx" +#include "SMESH_Gen_i.hxx" +#include "SMESH_TypeDefs.hxx" + +#ifndef GHS3D_VERSION +#define GHS3D_VERSION 41 +#endif //======================================================================= //function : GHS3DPlugin_Hypothesis_i //======================================================================= @@ -418,7 +427,7 @@ void GHS3DPlugin_Hypothesis_i::RemoveEnforcedVertex(CORBA::Double x, CORBA::Doub ExDescription.text = ex.what(); ExDescription.type = SALOME::BAD_PARAM; ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::RemoveEnforcedVertex(x,y,z)"; - ExDescription.lineNumber = 0; + ExDescription.lineNumber = 408; throw SALOME::SALOME_Exception(ExDescription); } catch (SALOME_Exception& ex) { @@ -434,8 +443,195 @@ void GHS3DPlugin_Hypothesis_i::ClearEnforcedVertices() { ASSERT(myBaseImpl); this->GetImpl()->ClearEnforcedVertices(); + SMESH::TPythonDump () << _this() << ".ClearEnforcedVertices() "; +} + +//======================================================================= +//function : ClearEnforcedMeshes +//======================================================================= + +void GHS3DPlugin_Hypothesis_i::ClearEnforcedMeshes() +{ + ASSERT(myBaseImpl); + this->GetImpl()->ClearEnforcedMeshes(); + SMESH::TPythonDump () << _this() << ".ClearEnforcedMeshes() "; +} + +/*! + * \brief Adds enforced elements of type elementType using another mesh/sub-mesh/mesh group theSource. + */ +void GHS3DPlugin_Hypothesis_i::SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType theType) + throw (SALOME::SALOME_Exception) +{ +#if GHS3D_VERSION >= 42 + _SetEnforcedMesh(theSource, theType, -1.0); + SMESH_Mesh_i* theMesh_i = SMESH::DownCast( theSource); + SMESH_Group_i* theGroup_i = SMESH::DownCast( theSource); + if (theGroup_i) + { + SMESH::TPythonDump () << _this() << ".SetEnforcedMesh( " + << theSource << ", " << theType << " )"; + } + else if (theMesh_i) + { + SMESH::TPythonDump () << _this() << ".SetEnforcedMesh( " + << theSource << ".GetMesh(), " << theType << " )"; + } +#else + SALOME::ExceptionStruct ExDescription; + ExDescription.text = "Bad version of GHS3D. It must >= 4.2."; + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::SetEnforcedMesh(theSource, theType)"; + ExDescription.lineNumber = 463; + throw SALOME::SALOME_Exception(ExDescription); +#endif +} + +/*! + * \brief Adds enforced elements of type elementType using another mesh/sub-mesh/mesh group theSource and a size. + */ +void GHS3DPlugin_Hypothesis_i::SetEnforcedMeshSize(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType theType, double theSize) + throw (SALOME::SALOME_Exception) +{ + if (theSize <= 0) { + SALOME::ExceptionStruct ExDescription; + ExDescription.text = "Size cannot be negative"; + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::SetEnforcedMeshSize(theSource, theType)"; + ExDescription.lineNumber = 475; + throw SALOME::SALOME_Exception(ExDescription); + } + + _SetEnforcedMesh(theSource, theType, theSize); + SMESH_Mesh_i* theMesh_i = SMESH::DownCast( theSource); + SMESH_Group_i* theGroup_i = SMESH::DownCast( theSource); + if (theGroup_i) + { + SMESH::TPythonDump () << _this() << ".SetEnforcedMeshSize( " + << theSource << ", " << theType << ", " << theSize << " )"; + } + else if (theMesh_i) + { + SMESH::TPythonDump () << _this() << ".SetEnforcedMeshSize( " + << theSource << ".GetMesh(), " << theType << ", " << theSize << " )"; + } } +void GHS3DPlugin_Hypothesis_i::_SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType theType, double theSize) + throw (SALOME::SALOME_Exception) +{ + ASSERT(myBaseImpl); + + if (CORBA::is_nil( theSource )) + { + SALOME::ExceptionStruct ExDescription; + ExDescription.text = "The source mesh CORBA object is NULL"; + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)"; + ExDescription.lineNumber = 502; + throw SALOME::SALOME_Exception(ExDescription); + } + + if ((theType != SMESH::NODE) && (theType != SMESH::EDGE) && (theType != SMESH::FACE)) + { + SALOME::ExceptionStruct ExDescription; + ExDescription.text = "Bad elementType"; + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)"; + ExDescription.lineNumber = 502; + throw SALOME::SALOME_Exception(ExDescription); + } + + SMESH::array_of_ElementType_var types = theSource->GetTypes(); +// MESSAGE("Required type is "<length();i++){MESSAGE(types[i]);} + if ( types->length() >= 1 && types[types->length()-1] < theType) + { + SALOME::ExceptionStruct ExDescription; + ExDescription.text = "The source mesh has bad type"; + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)"; + ExDescription.lineNumber = 502; + throw SALOME::SALOME_Exception(ExDescription); + } + + SMESHDS_Mesh* theMeshDS; + SMESH_Mesh_i* anImplPtr = SMESH::DownCast(theSource->GetMesh()); + if (anImplPtr) + theMeshDS = anImplPtr->GetImpl().GetMeshDS(); + else + return; + + SMESH_Mesh_i* theMesh_i = SMESH::DownCast( theSource); + SMESH_Group_i* theGroup_i = SMESH::DownCast( theSource); + TIDSortedElemSet theElemSet; + + if (theMesh_i) + { + try { + this->GetImpl()->SetEnforcedMesh(anImplPtr->GetImpl(), theType, theSize); + } + catch (const std::invalid_argument& ex) { + SALOME::ExceptionStruct ExDescription; + ExDescription.text = ex.what(); + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)"; + ExDescription.lineNumber = 502; + throw SALOME::SALOME_Exception(ExDescription); + } + catch (SALOME_Exception& ex) { + THROW_SALOME_CORBA_EXCEPTION( ex.what() ,SALOME::BAD_PARAM ); + } +// +// SMESH::long_array_var anIDs = theMesh_i->GetElementsByType(theType); +// if ( anIDs->length() == 0 ){MESSAGE("The source mesh is empty");} +// for (int i=0; ilength(); i++) { +// CORBA::Long ind = anIDs[i]; +// const SMDS_MeshElement * elem = theMeshDS->FindElement(ind); +// if (elem) +// theElemSet.insert( elem ); +// } +//// } +// MESSAGE("Add "<length() == 1 && types[0] == theType) + { + SMESH::long_array_var anIDs = theGroup_i->GetListOfID(); + if ( anIDs->length() == 0 ){MESSAGE("The source group is empty");} + for (int i=0; ilength(); i++) { + CORBA::Long ind = anIDs[i]; + if (theType == SMESH::NODE) + { + const SMDS_MeshNode * node = theMeshDS->FindNode(ind); + if (node) + theElemSet.insert( node ); + } + else + { + const SMDS_MeshElement * elem = theMeshDS->FindElement(ind); + if (elem) + theElemSet.insert( elem ); + } + } + MESSAGE("Add "<GetName()); + + try { + this->GetImpl()->SetEnforcedElements(theElemSet, theType, theSize); + } + catch (const std::invalid_argument& ex) { + SALOME::ExceptionStruct ExDescription; + ExDescription.text = ex.what(); + ExDescription.type = SALOME::BAD_PARAM; + ExDescription.sourceFile = "GHS3DPlugin_Hypothesis::_SetEnforcedMesh(theSource, theType, theSize)"; + ExDescription.lineNumber = 502; + throw SALOME::SALOME_Exception(ExDescription); + } + catch (SALOME_Exception& ex) { + THROW_SALOME_CORBA_EXCEPTION( ex.what() ,SALOME::BAD_PARAM ); + } + } +} //============================================================================= /*! * Get implementation diff --git a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.hxx b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.hxx index 0cf9a26..d9b5941 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.hxx +++ b/src/GHS3DPlugin/GHS3DPlugin_Hypothesis_i.hxx @@ -31,6 +31,7 @@ #include CORBA_SERVER_HEADER(GHS3DPlugin_Algorithm) #include "SMESH_Hypothesis_i.hxx" +#include "SMESH_Mesh_i.hxx" #include "GHS3DPlugin_Hypothesis.hxx" class SMESH_Gen; @@ -129,6 +130,13 @@ class GHS3DPLUGIN_EXPORT GHS3DPlugin_Hypothesis_i: void RemoveEnforcedVertex(CORBA::Double x, CORBA::Double y, CORBA::Double z) throw (SALOME::SALOME_Exception); GHS3DPlugin::GHS3DEnforcedVertexList* GetEnforcedVertices(); void ClearEnforcedVertices(); + /*! + * To set an enforced mesh + */ + void SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType elementType) throw (SALOME::SALOME_Exception); + void SetEnforcedMeshSize(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType elementType, double size) throw (SALOME::SALOME_Exception); + void _SetEnforcedMesh(SMESH::SMESH_IDSource_ptr theSource, SMESH::ElementType elementType, double size) throw (SALOME::SALOME_Exception); + void ClearEnforcedMeshes(); // Get implementation ::GHS3DPlugin_Hypothesis* GetImpl(); diff --git a/src/GHS3DPlugin/Makefile.am b/src/GHS3DPlugin/Makefile.am index 2875199..4e7dda2 100644 --- a/src/GHS3DPlugin/Makefile.am +++ b/src/GHS3DPlugin/Makefile.am @@ -27,6 +27,7 @@ include $(top_srcdir)/adm_local/unix/make_common_starter.am # header files salomeinclude_HEADERS = \ + libmesh5.h \ GHS3DPlugin_Defs.hxx \ GHS3DPlugin_GHS3D.hxx \ GHS3DPlugin_GHS3D_i.hxx \ @@ -37,6 +38,7 @@ salomeinclude_HEADERS = \ lib_LTLIBRARIES = libGHS3DEngine.la dist_libGHS3DEngine_la_SOURCES = \ + libmesh5.c \ GHS3DPlugin_GHS3D.cxx \ GHS3DPlugin_GHS3D_i.cxx \ GHS3DPlugin_i.cxx \ diff --git a/src/GHS3DPlugin/libmesh5.c b/src/GHS3DPlugin/libmesh5.c new file mode 100644 index 0000000..4277b97 --- /dev/null +++ b/src/GHS3DPlugin/libmesh5.c @@ -0,0 +1,1192 @@ + + +/*----------------------------------------------------------*/ +/* */ +/* LIBMESH V 5.45 */ +/* */ +/*----------------------------------------------------------*/ +/* */ +/* Description: handle .meshb file format I/O */ +/* Author: Loic MARECHAL */ +/* Creation date: feb 16 2007 */ +/* Last modification: feb 08 2011 */ +/* */ +/*----------------------------------------------------------*/ + + +/*----------------------------------------------------------*/ +/* Includes */ +/*----------------------------------------------------------*/ + +#include +#include +#include +#include +#include +#include +#include +#include "libmesh5.h" + + +/*----------------------------------------------------------*/ +/* Defines */ +/*----------------------------------------------------------*/ + +#define Asc 1 +#define Bin 2 +#define MshFil 4 +#define SolFil 8 +#define MaxMsh 100 +#define InfKwd 1 +#define RegKwd 2 +#define SolKwd 3 +#define WrdSiz 4 +#define BufSiz 10000 + + +/*----------------------------------------------------------*/ +/* Structures */ +/*----------------------------------------------------------*/ + +typedef struct +{ + int typ, SolSiz, NmbWrd, NmbLin, NmbTyp, TypTab[ GmfMaxTyp ]; + long pos; + char fmt[ GmfMaxTyp*9 ]; +}KwdSct; + +typedef struct +{ + int dim, ver, mod, typ, cod, pos; + long NexKwdPos; + KwdSct KwdTab[ GmfMaxKwd + 1 ]; + FILE *hdl; + int *IntBuf; + float *FltBuf; + unsigned char *buf; + char FilNam[ GmfStrSiz ]; + double DblBuf[1000/8]; + unsigned char blk[ BufSiz + 1000 ]; +}GmfMshSct; + + +/*----------------------------------------------------------*/ +/* Global variables */ +/*----------------------------------------------------------*/ + +int GmfIniFlg=0; +GmfMshSct *GmfMshTab[ MaxMsh + 1 ]; +char *GmfKwdFmt[ GmfMaxKwd + 1 ][4] = +{ {"Reserved", "", "", ""}, + {"MeshVersionFormatted", "", "", "i"}, + {"Reserved", "", "", ""}, + {"Dimension", "", "", "i"}, + {"Vertices", "Vertex", "i", "dri"}, + {"Edges", "Edge", "i", "iii"}, + {"Triangles", "Triangle", "i", "iiii"}, + {"Quadrilaterals", "Quadrilateral", "i", "iiiii"}, + {"Tetrahedra", "Tetrahedron", "i", "iiiii"}, + {"Prisms", "Prism", "i", "iiiiiii"}, + {"Hexahedra", "Hexahedron", "i", "iiiiiiiii"}, + {"IterationsAll", "IterationAll","","i"}, + {"TimesAll", "TimeAll","","r"}, + {"Corners", "Corner", "i", "i"}, + {"Ridges", "Ridge", "i", "i"}, + {"RequiredVertices", "RequiredVertex", "i", "i"}, + {"RequiredEdges", "RequiredEdge", "i", "i"}, + {"RequiredTriangles", "RequiredTriangle", "i", "i"}, + {"RequiredQuadrilaterals", "RequiredQuadrilateral", "i", "i"}, + {"TangentAtEdgeVertices", "TangentAtEdgeVertex", "i", "iii"}, + {"NormalAtVertices", "NormalAtVertex", "i", "ii"}, + {"NormalAtTriangleVertices", "NormalAtTriangleVertex", "i", "iii"}, + {"NormalAtQuadrilateralVertices", "NormalAtQuadrilateralVertex", "i", "iiii"}, + {"AngleOfCornerBound", "", "", "r"}, + {"TrianglesP2", "TriangleP2", "i", "iiiiiii"}, + {"EdgesP2", "EdgeP2", "i", "iiii"}, + {"SolAtPyramids", "SolAtPyramid", "i", "sr"}, + {"QuadrilateralsQ2", "QuadrilateralQ2", "i", "iiiiiiiiii"}, + {"ISolAtPyramids", "ISolAtPyramid", "i", "iiiii"}, + {"Reserved", "", "", ""}, + {"TetrahedraP2", "TetrahedronP2", "i", "iiiiiiiiiii"}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"HexahedraQ2", "HexahedronQ2", "i", "iiiiiiiiiiiiiiiiiiiiiiiiiiii"}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Polyhedra", "Polyhedron", "i", "iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii"}, + {"Polygons", "Polygon", "", "iiiiiiiii"}, + {"Reserved", "", "", ""}, + {"Pyramids", "Pyramid", "i", "iiiiii"}, + {"BoundingBox", "", "", "drdr"}, + {"Body","i", "drdrdrdr"}, + {"PrivateTable", "PrivateTable", "i", "i"}, + {"Reserved", "", "", ""}, + {"End", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Reserved", "", "", ""}, + {"Tangents", "Tangent", "i", "dr"}, + {"Normals", "Normal", "i", "dr"}, + {"TangentAtVertices", "TangentAtVertex", "i", "ii"}, + {"SolAtVertices", "SolAtVertex", "i", "sr"}, + {"SolAtEdges", "SolAtEdge", "i", "sr"}, + {"SolAtTriangles", "SolAtTriangle", "i", "sr"}, + {"SolAtQuadrilaterals", "SolAtQuadrilateral", "i", "sr"}, + {"SolAtTetrahedra", "SolAtTetrahedron", "i", "sr"}, + {"SolAtPrisms", "SolAtPrism", "i", "sr"}, + {"SolAtHexahedra", "SolAtHexahedron", "i", "sr"}, + {"DSolAtVertices", "DSolAtVertex", "i", "sr"}, + {"ISolAtVertices", "ISolAtVertex", "i", "i"}, + {"ISolAtEdges", "ISolAtEdge", "i", "ii"}, + {"ISolAtTriangles", "ISolAtTriangle", "i", "iii"}, + {"ISolAtQuadrilaterals", "ISolAtQuadrilateral", "i", "iiii"}, + {"ISolAtTetrahedra", "ISolAtTetrahedron", "i", "iiii"}, + {"ISolAtPrisms", "ISolAtPrism", "i", "iiiiii"}, + {"ISolAtHexahedra", "ISolAtHexahedron", "i", "iiiiiiii"}, + {"Iterations", "","","i"}, + {"Time", "","","r"}, + {"Reserved", "","",""} + }; + + +/*----------------------------------------------------------*/ +/* Prototypes of local procedures */ +/*----------------------------------------------------------*/ + +static void ScaWrd(GmfMshSct *, unsigned char *); +static void ScaDblWrd(GmfMshSct *, unsigned char *); +static void ScaBlk(GmfMshSct *, unsigned char *, int); +static long GetPos(GmfMshSct *); +static void RecWrd(GmfMshSct *, unsigned char *); +static void RecDblWrd(GmfMshSct *, unsigned char *); +static void RecBlk(GmfMshSct *, unsigned char *, int); +static void SetPos(GmfMshSct *, long); +static int ScaKwdTab(GmfMshSct *); +static void ExpFmt(GmfMshSct *, int); +static void ScaKwdHdr(GmfMshSct *, int); + + +/*----------------------------------------------------------*/ +/* Open a mesh file in read or write mod */ +/*----------------------------------------------------------*/ + +int GmfOpenMesh(const char *FilNam, int mod, ...) +{ + int i, KwdCod, res, *PtrVer, *PtrDim, MshIdx=0; + char str[ GmfStrSiz ]; + va_list VarArg; + GmfMshSct *msh; + + if(!GmfIniFlg) + { + for(i=0;i<=MaxMsh;i++) + GmfMshTab[i] = NULL; + + GmfIniFlg = 1; + } + + /*---------------------*/ + /* MESH STRUCTURE INIT */ + /*---------------------*/ + + for(i=1;i<=MaxMsh;i++) + if(!GmfMshTab[i]) + { + MshIdx = i; + break; + } + + if( !MshIdx || !(msh = calloc(1, sizeof(GmfMshSct))) ) + return(0); + + /* Copy the FilNam into the structure */ + + if(strlen(FilNam) + 7 >= GmfStrSiz) + return(0); + + strcpy(msh->FilNam, FilNam); + + /* Store the opening mod (read or write) and guess the filetype (binary or ascii) depending on the extension */ + + msh->mod = mod; + msh->buf = (void *)msh->DblBuf; + msh->FltBuf = (void *)msh->DblBuf; + msh->IntBuf = (void *)msh->DblBuf; + + if(strstr(msh->FilNam, ".meshb")) + msh->typ |= (Bin | MshFil); + else if(strstr(msh->FilNam, ".mesh")) + msh->typ |= (Asc | MshFil); + else if(strstr(msh->FilNam, ".solb")) + msh->typ |= (Bin | SolFil); + else if(strstr(msh->FilNam, ".sol")) + msh->typ |= (Asc | SolFil); + else + return(0); + + /* Open the file in the required mod and initialyse the mesh structure */ + + if(msh->mod == GmfRead) + { + + /*-----------------------*/ + /* OPEN FILE FOR READING */ + /*-----------------------*/ + + va_start(VarArg, mod); + PtrVer = va_arg(VarArg, int *); + PtrDim = va_arg(VarArg, int *); + va_end(VarArg); + + /* Create the name string and open the file */ + + if(!(msh->hdl = fopen(msh->FilNam, "rb"))) + return(0); + + /* Read the endian coding tag, the mesh version and the mesh dimension (mandatory kwd) */ + + if(msh->typ & Bin) + { + fread((unsigned char *)&msh->cod, WrdSiz, 1, msh->hdl); + + if( (msh->cod != 1) && (msh->cod != 16777216) ) + return(0); + + ScaWrd(msh, (unsigned char *)&msh->ver); + + if( (msh->ver < 1) || (msh->ver > 3) ) + return(0); + + if( (msh->ver == 3) && (sizeof(long) == 4) ) + return(0); + + ScaWrd(msh, (unsigned char *)&KwdCod); + + if(KwdCod != GmfDimension) + return(0); + + GetPos(msh); + ScaWrd(msh, (unsigned char *)&msh->dim); + } + else + { + do + { + res = fscanf(msh->hdl, "%s", str); + }while( (res != EOF) && strcmp(str, "MeshVersionFormatted") ); + + if(res == EOF) + return(0); + + fscanf(msh->hdl, "%d", &msh->ver); + + if( (msh->ver < 1) || (msh->ver > 3) ) + return(0); + + do + { + res = fscanf(msh->hdl, "%s", str); + }while( (res != EOF) && strcmp(str, "Dimension") ); + + if(res == EOF) + return(0); + + fscanf(msh->hdl, "%d", &msh->dim); + } + + if( (msh->dim != 2) && (msh->dim != 3) ) + return(0); + + (*PtrVer) = msh->ver; + (*PtrDim) = msh->dim; + + /*------------*/ + /* KW READING */ + /*------------*/ + + /* Read the list of kw present in the file */ + + if(!ScaKwdTab(msh)) + return(0); + + GmfMshTab[ MshIdx ] = msh; + + return(MshIdx); + } + else if(msh->mod == GmfWrite) + { + + /*-----------------------*/ + /* OPEN FILE FOR WRITING */ + /*-----------------------*/ + + msh->cod = 1; + + /* Check if the user provided a valid version number and dimension */ + + va_start(VarArg, mod); + msh->ver = va_arg(VarArg, int); + msh->dim = va_arg(VarArg, int); + va_end(VarArg); + + if( (msh->ver < 1) || (msh->ver > 3) ) + return(0); + + if( (msh->ver == 3) && (sizeof(long) == 4) ) + return(0); + + if( (msh->dim != 2) && (msh->dim != 3) ) + return(0); + + /* Create the mesh file */ + + if(!(msh->hdl = fopen(msh->FilNam, "wb"))) + return(0); + + GmfMshTab[ MshIdx ] = msh; + + + /*------------*/ + /* KW WRITING */ + /*------------*/ + + /* Write the mesh version and dimension */ + + if(msh->typ & Asc) + { + fprintf(msh->hdl, "%s %d\n\n", GmfKwdFmt[ GmfVersionFormatted ][0], msh->ver); + fprintf(msh->hdl, "%s %d\n", GmfKwdFmt[ GmfDimension ][0], msh->dim); + } + else + { + RecWrd(msh, (unsigned char *)&msh->cod); + RecWrd(msh, (unsigned char *)&msh->ver); + GmfSetKwd(MshIdx, GmfDimension, 0); + RecWrd(msh, (unsigned char *)&msh->dim); + } + + return(MshIdx); + } + else + return(0); +} + + +/*----------------------------------------------------------*/ +/* Close a meshfile in the right way */ +/*----------------------------------------------------------*/ + +int GmfCloseMesh(int MshIdx) +{ + int res = 1; + GmfMshSct *msh; + + if( (MshIdx < 1) || (MshIdx > MaxMsh) ) + return(0); + + msh = GmfMshTab[ MshIdx ]; + RecBlk(msh, msh->buf, 0); + + /* In write down the "End" kw in write mode */ + + if(msh->mod == GmfWrite) { + if(msh->typ & Asc) + fprintf(msh->hdl, "\n%s\n", GmfKwdFmt[ GmfEnd ][0]); + else + GmfSetKwd(MshIdx, GmfEnd, 0); + } + + /* Close the file and free the mesh structure */ + + if(fclose(msh->hdl)) + res = 0; + + free(msh); + GmfMshTab[ MshIdx ] = NULL; + + return(res); +} + + +/*----------------------------------------------------------*/ +/* Read the number of lines and set the position to this kwd*/ +/*----------------------------------------------------------*/ + +int GmfStatKwd(int MshIdx, int KwdCod, ...) +{ + int i, *PtrNmbTyp, *PtrSolSiz, *TypTab; + GmfMshSct *msh; + KwdSct *kwd; + va_list VarArg; + + if( (MshIdx < 1) || (MshIdx > MaxMsh) ) + return(0); + + msh = GmfMshTab[ MshIdx ]; + + if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) ) + return(0); + + kwd = &msh->KwdTab[ KwdCod ]; + + if(!kwd->NmbLin) + return(0); + + /* Read further arguments if this kw is a sol */ + + if(kwd->typ == SolKwd) + { + va_start(VarArg, KwdCod); + + PtrNmbTyp = va_arg(VarArg, int *); + *PtrNmbTyp = kwd->NmbTyp; + + PtrSolSiz = va_arg(VarArg, int *); + *PtrSolSiz = kwd->SolSiz; + + TypTab = va_arg(VarArg, int *); + + for(i=0;iNmbTyp;i++) + TypTab[i] = kwd->TypTab[i]; + + va_end(VarArg); + } + + return(kwd->NmbLin); +} + + +/*----------------------------------------------------------*/ +/* Set the current file position to a given kwd */ +/*----------------------------------------------------------*/ + +int GmfGotoKwd(int MshIdx, int KwdCod) +{ + GmfMshSct *msh; + KwdSct *kwd; + + if( (MshIdx < 1) || (MshIdx > MaxMsh) ) + return(0); + + msh = GmfMshTab[ MshIdx ]; + + if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) ) + return(0); + + kwd = &msh->KwdTab[ KwdCod ]; + + if(!kwd->NmbLin) + return(0); + + return(fseek(msh->hdl, kwd->pos, SEEK_SET)); +} + + +/*----------------------------------------------------------*/ +/* Write the kwd and set the number of lines */ +/*----------------------------------------------------------*/ + +int GmfSetKwd(int MshIdx, int KwdCod, ...) +{ + int i, NmbLin=0, *TypTab; + long CurPos; + va_list VarArg; + GmfMshSct *msh; + KwdSct *kwd; + + if( (MshIdx < 1) || (MshIdx > MaxMsh) ) + return(0); + + msh = GmfMshTab[ MshIdx ]; + RecBlk(msh, msh->buf, 0); + + if( (KwdCod < 1) || (KwdCod > GmfMaxKwd) ) + return(0); + + kwd = &msh->KwdTab[ KwdCod ]; + + /* Read further arguments if this kw has a header */ + + if(strlen(GmfKwdFmt[ KwdCod ][2])) + { + va_start(VarArg, KwdCod); + NmbLin = va_arg(VarArg, int); + + if(!strcmp(GmfKwdFmt[ KwdCod ][3], "sr")) + { + kwd->NmbTyp = va_arg(VarArg, int); + TypTab = va_arg(VarArg, int *); + + for(i=0;iNmbTyp;i++) + kwd->TypTab[i] = TypTab[i]; + } + + va_end(VarArg); + } + + /* Setup the kwd info */ + + ExpFmt(msh, KwdCod); + + if(!kwd->typ) + return(0); + else if(kwd->typ == InfKwd) + kwd->NmbLin = 1; + else + kwd->NmbLin = NmbLin; + + /* Store the next kwd position in binary file */ + + if( (msh->typ & Bin) && msh->NexKwdPos ) + { + CurPos = ftell(msh->hdl); + fseek(msh->hdl, msh->NexKwdPos, SEEK_SET); + SetPos(msh, CurPos); + fseek(msh->hdl, CurPos, SEEK_SET); + } + + /* Write the header */ + + if(msh->typ & Asc) + { + fprintf(msh->hdl, "\n%s\n", GmfKwdFmt[ KwdCod ][0]); + + if(kwd->typ != InfKwd) + fprintf(msh->hdl, "%d\n", kwd->NmbLin); + + /* In case of solution field, write the extended header */ + + if(kwd->typ == SolKwd) + { + fprintf(msh->hdl, "%d ", kwd->NmbTyp); + + for(i=0;iNmbTyp;i++) + fprintf(msh->hdl, "%d ", kwd->TypTab[i]); + + fprintf(msh->hdl, "\n\n"); + } + } + else + { + RecWrd(msh, (unsigned char *)&KwdCod); + msh->NexKwdPos = ftell(msh->hdl); + SetPos(msh, 0); + + if(kwd->typ != InfKwd) + RecWrd(msh, (unsigned char *)&kwd->NmbLin); + + /* In case of solution field, write the extended header at once */ + + if(kwd->typ == SolKwd) + { + RecWrd(msh, (unsigned char *)&kwd->NmbTyp); + + for(i=0;iNmbTyp;i++) + RecWrd(msh, (unsigned char *)&kwd->TypTab[i]); + } + } + msh->pos = 0; + return(kwd->NmbLin); +} + + +/*----------------------------------------------------------*/ +/* Read a full line from the current kwd */ +/*----------------------------------------------------------*/ + +void GmfGetLin(int MshIdx, int KwdCod, ...) +{ + int i, j; + float *FltSolTab; + double *DblSolTab; + va_list VarArg; + GmfMshSct *msh = GmfMshTab[ MshIdx ]; + KwdSct *kwd = &msh->KwdTab[ KwdCod ]; + + /* Start decoding the arguments */ + + va_start(VarArg, KwdCod); + + if(kwd->typ != SolKwd) + { + if(msh->ver == 1) + { + if(msh->typ & Asc) + { + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + fscanf(msh->hdl, "%f", va_arg(VarArg, float *)); + else + fscanf(msh->hdl, "%d", va_arg(VarArg, int *)); + } + else + { + ScaBlk(msh, msh->buf, kwd->SolSiz); + + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + *(va_arg(VarArg, float *)) = msh->FltBuf[i]; + else + *(va_arg(VarArg, int *)) = msh->IntBuf[i]; + } + } + else + { + if(msh->typ & Asc) + { + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + fscanf(msh->hdl, "%lf", va_arg(VarArg, double *)); + else + fscanf(msh->hdl, "%d", va_arg(VarArg, int *)); + } + else + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + ScaDblWrd(msh, (unsigned char *)va_arg(VarArg, double *)); + else + ScaWrd(msh, (unsigned char *)va_arg(VarArg, int *)); + } + } + else + { + if(msh->ver == 1) + { + FltSolTab = va_arg(VarArg, float *); + + if(msh->typ & Asc) + for(j=0;jSolSiz;j++) + fscanf(msh->hdl, "%f", &FltSolTab[j]); + else + ScaBlk(msh, (unsigned char *)FltSolTab, kwd->NmbWrd); + } + else + { + DblSolTab = va_arg(VarArg, double *); + + if(msh->typ & Asc) + for(j=0;jSolSiz;j++) + fscanf(msh->hdl, "%lf", &DblSolTab[j]); + else + for(j=0;jSolSiz;j++) + ScaDblWrd(msh, (unsigned char *)&DblSolTab[j]); + } + } + + va_end(VarArg); +} + + +/*----------------------------------------------------------*/ +/* Write a full line from the current kwd */ +/*----------------------------------------------------------*/ + +void GmfSetLin(int MshIdx, int KwdCod, ...) +{ + int i, j, pos, *IntBuf; + float *FltSolTab; + double *DblSolTab, *DblBuf; + va_list VarArg; + GmfMshSct *msh = GmfMshTab[ MshIdx ]; + KwdSct *kwd = &msh->KwdTab[ KwdCod ]; + + /* Start decoding the arguments */ + + va_start(VarArg, KwdCod); + + if(kwd->typ != SolKwd) + { + if(msh->ver == 1) + { + if(msh->typ & Asc) + { + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + fprintf(msh->hdl, "%g ", (float)va_arg(VarArg, double)); + else + fprintf(msh->hdl, "%d ", va_arg(VarArg, int)); + } + else + { + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + msh->FltBuf[i] = va_arg(VarArg, double); + else + msh->IntBuf[i] = va_arg(VarArg, int); + + RecBlk(msh, msh->buf, kwd->SolSiz); + } + } + else + { + if(msh->typ & Asc) + { + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + fprintf(msh->hdl, "%.15lg ", va_arg(VarArg, double)); + else + fprintf(msh->hdl, "%d ", va_arg(VarArg, int)); + } + else + { + pos = 0; + + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'r') + { + DblBuf = (void *)&msh->buf[ pos ]; + *DblBuf = va_arg(VarArg, double); + pos += 8; + } + else + { + IntBuf = (void *)&msh->buf[ pos ]; + *IntBuf = va_arg(VarArg, int); + pos += 4; + } + RecBlk(msh, msh->buf, kwd->NmbWrd); + } + } + } + else + { + if(msh->ver == 1) + { + FltSolTab = va_arg(VarArg, float *); + + if(msh->typ & Asc) + for(j=0;jSolSiz;j++) + fprintf(msh->hdl, "%g ", FltSolTab[j]); + else + RecBlk(msh, (unsigned char *)FltSolTab, kwd->NmbWrd); + } + else + { + DblSolTab = va_arg(VarArg, double *); + + if(msh->typ & Asc) + for(j=0;jSolSiz;j++) + fprintf(msh->hdl, "%.15lg ", DblSolTab[j]); + else + RecBlk(msh, (unsigned char *)DblSolTab, kwd->NmbWrd); + } + } + + va_end(VarArg); + + if(msh->typ & Asc) + fprintf(msh->hdl, "\n"); +} + + +/*----------------------------------------------------------*/ +/* Private procedure for transmesh : copy a whole line */ +/*----------------------------------------------------------*/ + +void GmfCpyLin(int InpIdx, int OutIdx, int KwdCod) +{ + double d; + float f; + int i, a; + GmfMshSct *InpMsh = GmfMshTab[ InpIdx ], *OutMsh = GmfMshTab[ OutIdx ]; + KwdSct *kwd = &InpMsh->KwdTab[ KwdCod ]; + + for(i=0;iSolSiz;i++) + { + if(kwd->fmt[i] == 'r') + { + if(InpMsh->ver == 1) + { + if(InpMsh->typ & Asc) + fscanf(InpMsh->hdl, "%f", &f); + else + ScaWrd(InpMsh, (unsigned char *)&f); + + d = f; + } + else + { + if(InpMsh->typ & Asc) + fscanf(InpMsh->hdl, "%lf", &d); + else + ScaDblWrd(InpMsh, (unsigned char *)&d); + + f = (float)d; + } + + if(OutMsh->ver == 1) + if(OutMsh->typ & Asc) + fprintf(OutMsh->hdl, "%g ", f); + else + RecWrd(OutMsh, (unsigned char *)&f); + else + if(OutMsh->typ & Asc) + fprintf(OutMsh->hdl, "%.15g ", d); + else + RecDblWrd(OutMsh, (unsigned char *)&d); + } + else + { + if(InpMsh->typ & Asc) + fscanf(InpMsh->hdl, "%d", &a); + else + ScaWrd(InpMsh, (unsigned char *)&a); + + if(OutMsh->typ & Asc) + fprintf(OutMsh->hdl, "%d ", a); + else + RecWrd(OutMsh, (unsigned char *)&a); + } + } + + if(OutMsh->typ & Asc) + fprintf(OutMsh->hdl, "\n"); +} + + +/*----------------------------------------------------------*/ +/* Find every kw present in a meshfile */ +/*----------------------------------------------------------*/ + +static int ScaKwdTab(GmfMshSct *msh) +{ + int KwdCod; + long NexPos, CurPos, EndPos; + char str[ GmfStrSiz ]; + + if(msh->typ & Asc) + { + /* Scan each string in the file until the end */ + + while(fscanf(msh->hdl, "%s", str) != EOF) + { + /* Fast test in order to reject quickly the numeric values */ + + if(isalpha(str[0])) + { + /* Search which kwd code this string is associated with, + then get its header and save the curent position in file (just before the data) */ + + for(KwdCod=1; KwdCod<= GmfMaxKwd; KwdCod++) + if(!strcmp(str, GmfKwdFmt[ KwdCod ][0])) + { + ScaKwdHdr(msh, KwdCod); + break; + } + } + else if(str[0] == '#') + while(fgetc(msh->hdl) != '\n'); + } + } + else + { + /* Get file size */ + + CurPos = ftell(msh->hdl); + fseek(msh->hdl, 0, SEEK_END); + EndPos = ftell(msh->hdl); + fseek(msh->hdl, CurPos, SEEK_SET); + + /* Jump through kwd positions in the file */ + + do + { + /* Get the kwd code and the next kwd position */ + + ScaWrd(msh, (unsigned char *)&KwdCod); + NexPos = GetPos(msh); + + if(NexPos > EndPos) + return(0); + + /* Check if this kwd belongs to this mesh version */ + + if( (KwdCod >= 1) && (KwdCod <= GmfMaxKwd) ) + ScaKwdHdr(msh, KwdCod); + + /* Go to the next kwd */ + + if(NexPos) + fseek(msh->hdl, NexPos, SEEK_SET); + }while(NexPos && (KwdCod != GmfEnd)); + } + + return(1); +} + + +/*----------------------------------------------------------*/ +/* Read and setup the keyword's header */ +/*----------------------------------------------------------*/ + +static void ScaKwdHdr(GmfMshSct *msh, int KwdCod) +{ + int i; + KwdSct *kwd = &msh->KwdTab[ KwdCod ]; + + if(!strcmp("i", GmfKwdFmt[ KwdCod ][2])) + { + if(msh->typ & Asc) + fscanf(msh->hdl, "%d", &kwd->NmbLin); + else + ScaWrd(msh, (unsigned char *)&kwd->NmbLin); + } + else + kwd->NmbLin = 1; + + if(!strcmp("sr", GmfKwdFmt[ KwdCod ][3])) + { + if(msh->typ & Asc) + { + fscanf(msh->hdl, "%d", &kwd->NmbTyp); + + for(i=0;iNmbTyp;i++) + fscanf(msh->hdl, "%d", &kwd->TypTab[i]); + } + else + { + ScaWrd(msh, (unsigned char *)&kwd->NmbTyp); + + for(i=0;iNmbTyp;i++) + ScaWrd(msh, (unsigned char *)&kwd->TypTab[i]); + } + } + + ExpFmt(msh, KwdCod); + kwd->pos = ftell(msh->hdl); +} + + +/*----------------------------------------------------------*/ +/* Expand the compacted format and compute the line size */ +/*----------------------------------------------------------*/ + +static void ExpFmt(GmfMshSct *msh, int KwdCod) +{ + int i, j, TmpSiz=0; + char chr, *InpFmt = GmfKwdFmt[ KwdCod ][3]; + KwdSct *kwd = &msh->KwdTab[ KwdCod ]; + + /* Set the kwd's type */ + + if(!strlen(GmfKwdFmt[ KwdCod ][2])) + kwd->typ = InfKwd; + else if(!strcmp(InpFmt, "sr")) + kwd->typ = SolKwd; + else + kwd->typ = RegKwd; + + /* Get the solution-field's size */ + + if(kwd->typ == SolKwd) + for(i=0;iNmbTyp;i++) + switch(kwd->TypTab[i]) + { + case GmfSca : TmpSiz += 1; break; + case GmfVec : TmpSiz += msh->dim; break; + case GmfSymMat : TmpSiz += (msh->dim * (msh->dim+1)) / 2; break; + case GmfMat : TmpSiz += msh->dim * msh->dim; break; + } + + /* Scan each character from the format string */ + + i = kwd->SolSiz = kwd->NmbWrd = 0; + + while(i < strlen(InpFmt)) + { + chr = InpFmt[ i++ ]; + + if(chr == 'd') + { + chr = InpFmt[i++]; + + for(j=0;jdim;j++) + kwd->fmt[ kwd->SolSiz++ ] = chr; + } + else if(chr == 's') + { + chr = InpFmt[i++]; + + for(j=0;jfmt[ kwd->SolSiz++ ] = chr; + } + else + kwd->fmt[ kwd->SolSiz++ ] = chr; + } + + for(i=0;iSolSiz;i++) + if(kwd->fmt[i] == 'i') + kwd->NmbWrd++; + else if(msh->ver >= 2) + kwd->NmbWrd += 2; + else + kwd->NmbWrd++; +} + + +/*----------------------------------------------------------*/ +/* Read a four bytes word from a mesh file */ +/*----------------------------------------------------------*/ + +static void ScaWrd(GmfMshSct *msh, unsigned char *wrd) +{ + unsigned char swp; + + fread(wrd, WrdSiz, 1, msh->hdl); + + if(msh->cod == 1) + return; + + swp = wrd[3]; + wrd[3] = wrd[0]; + wrd[0] = swp; + + swp = wrd[2]; + wrd[2] = wrd[1]; + wrd[1] = swp; +} + + +/*----------------------------------------------------------*/ +/* Read an eight bytes word from a mesh file */ +/*----------------------------------------------------------*/ + +static void ScaDblWrd(GmfMshSct *msh, unsigned char *wrd) +{ + int i; + unsigned char swp; + + fread(wrd, WrdSiz, 2, msh->hdl); + + if(msh->cod == 1) + return; + + for(i=0;i<4;i++) + { + swp = wrd[7-i]; + wrd[7-i] = wrd[i]; + wrd[i] = swp; + } +} + + +/*----------------------------------------------------------*/ +/* Read ablock of four bytes word from a mesh file */ +/*----------------------------------------------------------*/ + +static void ScaBlk(GmfMshSct *msh, unsigned char *blk, int siz) +{ + int i, j; + unsigned char swp, *wrd; + + fread(blk, WrdSiz, siz, msh->hdl); + + if(msh->cod == 1) + return; + + for(i=0;iver >= 3) + ScaDblWrd(msh, (unsigned char*)&pos); + else + { + ScaWrd(msh, (unsigned char*)&IntVal); + pos = IntVal; + } + + return(pos); +} + + +/*----------------------------------------------------------*/ +/* Write a four bytes word to a mesh file */ +/*----------------------------------------------------------*/ + +static void RecWrd(GmfMshSct *msh, unsigned char *wrd) +{ + fwrite(wrd, WrdSiz, 1, msh->hdl); +} + + +/*----------------------------------------------------------*/ +/* Write an eight bytes word to a mesh file */ +/*----------------------------------------------------------*/ + +static void RecDblWrd(GmfMshSct *msh, unsigned char *wrd) +{ + fwrite(wrd, WrdSiz, 2, msh->hdl); +} + + +/*----------------------------------------------------------*/ +/* Write a block of four bytes word to a mesh file */ +/*----------------------------------------------------------*/ + +static void RecBlk(GmfMshSct *msh, unsigned char *blk, int siz) +{ + /* Copy this line-block into the main mesh buffer */ + + if(siz) + { + memcpy(&msh->blk[ msh->pos ], blk, siz * WrdSiz); + msh->pos += siz * WrdSiz; + } + + /* When the buffer is full or this procedure is called with a 0 size, flush the cache on disk */ + + if( (msh->pos > BufSiz) || (!siz && msh->pos) ) + { + fwrite(msh->blk, 1, msh->pos, msh->hdl); + msh->pos = 0; + } +} + + +/*----------------------------------------------------------*/ +/* Read a 4 or 8 bytes position in mesh file */ +/*----------------------------------------------------------*/ + +static void SetPos(GmfMshSct *msh, long pos) +{ + int IntVal; + + if(msh->ver >= 3) + RecDblWrd(msh, (unsigned char*)&pos); + else + { + IntVal = pos; + RecWrd(msh, (unsigned char*)&IntVal); + } +} diff --git a/src/GHS3DPlugin/libmesh5.h b/src/GHS3DPlugin/libmesh5.h new file mode 100755 index 0000000..cfaf9fb --- /dev/null +++ b/src/GHS3DPlugin/libmesh5.h @@ -0,0 +1,152 @@ + + +/*----------------------------------------------------------*/ +/* */ +/* LIBMESH V 5.45 */ +/* */ +/*----------------------------------------------------------*/ +/* */ +/* Description: handle .meshb file format I/O */ +/* Author: Loic MARECHAL */ +/* Creation date: feb 16 2007 */ +/* Last modification: sep 27 2010 */ +/* */ +/*----------------------------------------------------------*/ + + +/*----------------------------------------------------------*/ +/* Defines */ +/*----------------------------------------------------------*/ + +#define GmfStrSiz 1024 +#define GmfMaxTyp 1000 +#define GmfMaxKwd 79 +#define GmfMshVer 1 +#define GmfRead 1 +#define GmfWrite 2 +#define GmfSca 1 +#define GmfVec 2 +#define GmfSymMat 3 +#define GmfMat 4 +#define GmfFloat 1 +#define GmfDouble 2 + +enum GmfKwdCod +{ + GmfReserved1, \ + GmfVersionFormatted, \ + GmfReserved2, \ + GmfDimension, \ + GmfVertices, \ + GmfEdges, \ + GmfTriangles, \ + GmfQuadrilaterals, \ + GmfTetrahedra, \ + GmfPrisms, \ + GmfHexahedra, \ + GmfIterationsAll, \ + GmfTimesAll, \ + GmfCorners, \ + GmfRidges, \ + GmfRequiredVertices, \ + GmfRequiredEdges, \ + GmfRequiredTriangles, \ + GmfRequiredQuadrilaterals, \ + GmfTangentAtEdgeVertices, \ + GmfNormalAtVertices, \ + GmfNormalAtTriangleVertices, \ + GmfNormalAtQuadrilateralVertices, \ + GmfAngleOfCornerBound, \ + GmfTrianglesP2, \ + GmfEdgesP2, \ + GmfSolAtPyramids, \ + GmfQuadrilateralsQ2, \ + GmfISolAtPyramids, \ + GmfReserved6, \ + GmfTetrahedraP2, \ + GmfReserved7, \ + GmfReserved8, \ + GmfHexahedraQ2, \ + GmfReserved9, \ + GmfReserved10, \ + GmfReserved17, \ + GmfReserved18, \ + GmfReserved19, \ + GmfReserved20, \ + GmfReserved21, \ + GmfReserved22, \ + GmfReserved23, \ + GmfReserved24, \ + GmfReserved25, \ + GmfReserved26, \ + GmfPolyhedra, \ + GmfPolygons, \ + GmfReserved29, \ + GmfPyramids, \ + GmfBoundingBox, \ + GmfBody, \ + GmfPrivateTable, \ + GmfReserved33, \ + GmfEnd, \ + GmfReserved34, \ + GmfReserved35, \ + GmfReserved36, \ + GmfReserved37, \ + GmfTangents, \ + GmfNormals, \ + GmfTangentAtVertices, \ + GmfSolAtVertices, \ + GmfSolAtEdges, \ + GmfSolAtTriangles, \ + GmfSolAtQuadrilaterals, \ + GmfSolAtTetrahedra, \ + GmfSolAtPrisms, \ + GmfSolAtHexahedra, \ + GmfDSolAtVertices, \ + GmfISolAtVertices, \ + GmfISolAtEdges, \ + GmfISolAtTriangles, \ + GmfISolAtQuadrilaterals, \ + GmfISolAtTetrahedra, \ + GmfISolAtPrisms, \ + GmfISolAtHexahedra, \ + GmfIterations, \ + GmfTime, \ + GmfReserved38 +}; + + +/*----------------------------------------------------------*/ +/* External procedures */ +/*----------------------------------------------------------*/ + +extern int GmfOpenMesh(const char *, int, ...); +extern int GmfCloseMesh(int); +extern int GmfStatKwd(int, int, ...); +extern int GmfGotoKwd(int, int); +extern int GmfSetKwd(int, int, ...); +extern void GmfGetLin(int, int, ...); +extern void GmfSetLin(int, int, ...); + + +/*----------------------------------------------------------*/ +/* Fortran 77 API */ +/*----------------------------------------------------------*/ + +#if defined(F77_NO_UNDER_SCORE) +#define call(x) x +#else +#define call(x) x ## _ +#endif + + +/*----------------------------------------------------------*/ +/* Transmesh private API */ +/*----------------------------------------------------------*/ + +#ifdef TRANSMESH + +extern char *GmfKwdFmt[ GmfMaxKwd + 1 ][4]; +extern int GmfCpyLin(int, int, int); + +#endif