X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FGHS3DPlugin%2FGHS3DPlugin_GHS3D.cxx;h=a4498543c8fd5c794a5ea680f4c0cb3dd789a9bc;hb=6a9f3ba14d636dfcaf9d00dd044573beeabeaf8a;hp=b8a777e4a478602f2ef0ca9386fd692f1991cdb1;hpb=09d65a759f378f5eb9b23c0fd7d281fecb882e36;p=plugins%2Fghs3dplugin.git diff --git a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx index b8a777e..a449854 100644 --- a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx +++ b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx @@ -1,9 +1,9 @@ -// Copyright (C) 2004-2013 CEA/DEN, EDF R&D +// Copyright (C) 2004-2016 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -26,34 +26,32 @@ // #include "GHS3DPlugin_GHS3D.hxx" #include "GHS3DPlugin_Hypothesis.hxx" +#include "MG_Tetra_API.hxx" -#include - -//#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - +#include +#include #include #include -#include #include - #include - +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include #include #include +#include #include #include #include @@ -62,9 +60,7 @@ #include #include #include -#include #include -#include #include #include #include @@ -74,50 +70,58 @@ #include #include #include -#include +#include #include +#include #include -#ifdef WIN32 -#include -#else -#include -#endif #include +#include -//#include - +#ifdef _DEBUG_ +//#define _MY_DEBUG_ +#endif #define castToNode(n) static_cast( n ); -#ifdef _DEBUG_ -#define DUMP(txt) \ -// std::cout << txt -#else -#define DUMP(txt) -#endif +using namespace std; -extern "C" +#define HOLE_ID -1 + +// flags returning state of enforced entities, returned from writeGMFFile +enum InvalidEnforcedFlags { FLAG_BAD_ENF_VERT = 1, + FLAG_BAD_ENF_NODE = 2, + FLAG_BAD_ENF_EDGE = 4, + FLAG_BAD_ENF_TRIA = 8 +}; +static std::string flagsToErrorStr( int anInvalidEnforcedFlags ) { -#ifndef WNT -#include -#include -#endif -#include -#include + std::string str; + if ( anInvalidEnforcedFlags != 0 ) + { + if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_VERT ) + str = "There are enforced vertices incorrectly defined.\n"; + if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_NODE ) + str += "There are enforced nodes incorrectly defined.\n"; + if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_EDGE ) + str += "There are enforced edge incorrectly defined.\n"; + if ( anInvalidEnforcedFlags & FLAG_BAD_ENF_TRIA ) + str += "There are enforced triangles incorrectly defined.\n"; + } + return str; } -#define HOLE_ID -1 - typedef const list TTriaList; +static const char theDomainGroupNamePrefix[] = "Domain_"; + static void removeFile( const TCollection_AsciiString& fileName ) { try { - OSD_File( fileName ).Remove(); + SMESH_File( fileName.ToCString() ).remove(); } - catch ( Standard_ProgramError ) { + catch ( ... ) { MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied"); } } @@ -128,49 +132,35 @@ static void removeFile( const TCollection_AsciiString& fileName ) */ //============================================================================= -GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen) - : SMESH_3D_Algo(hypId, studyId, gen) +GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, SMESH_Gen* gen) + : SMESH_3D_Algo(hypId, gen), _isLibUsed( false ) { - MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D"); - _name = "GHS3D_3D"; + _name = Name(); _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type _onlyUnaryInput = false; // Compute() will be called on a compound of solids _iShape=0; _nbShape=0; _compatibleHypothesis.push_back( GHS3DPlugin_Hypothesis::GetHypType()); _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() ); - _requireShape = false; // can work without shape_studyId - - smeshGen_i = SMESH_Gen_i::GetSMESHGen(); - CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager"); - SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject); - - MESSAGE("studyid = " << _studyId); - - myStudy = NULL; - myStudy = aStudyMgr->GetStudyByID(_studyId); - if (myStudy) - MESSAGE("myStudy->StudyId() = " << myStudy->StudyId()); + _requireShape = false; // can work without shape -#ifdef WITH_SMESH_CANCEL_COMPUTE - _compute_canceled = false; -#endif + _computeCanceled = false; + _progressAdvance = 1e-4; } //============================================================================= /*! - * + * */ //============================================================================= GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D() { - MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D"); } //============================================================================= /*! - * + * */ //============================================================================= @@ -183,6 +173,8 @@ bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh& aMesh, _hyp = 0; _viscousLayersHyp = 0; _keepFiles = false; + _removeLogOnSuccess = true; + _logInStandardOutput = false; const list & hyps = GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false); @@ -195,9 +187,16 @@ bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh& aMesh, _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h ); } if ( _hyp ) - _keepFiles = _hyp->GetKeepFiles(); + { + _keepFiles = _hyp->GetKeepFiles(); + _removeLogOnSuccess = _hyp->GetRemoveLogOnSuccess(); + _logInStandardOutput = _hyp->GetStandardOutputLog(); + } - return true; + if ( _viscousLayersHyp ) + error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus )); + + return aStatus == HYP_OK; } @@ -208,124 +207,76 @@ bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh& aMesh, TopoDS_Shape GHS3DPlugin_GHS3D::entryToShape(std::string entry) { - MESSAGE("GHS3DPlugin_GHS3D::entryToShape "<_is_nil() ) + throw SALOME_Exception("MG-Tetra plugin can't work w/o publishing in the study"); + GEOM::GEOM_Object_var aGeomObj; TopoDS_Shape S = TopoDS_Shape(); - SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() ); + SALOMEDS::SObject_var aSObj = SMESH_Gen_i::getStudyServant()->FindObjectID( entry.c_str() ); if (!aSObj->_is_nil() ) { CORBA::Object_var obj = aSObj->GetObject(); aGeomObj = GEOM::GEOM_Object::_narrow(obj); aSObj->UnRegister(); } if ( !aGeomObj->_is_nil() ) - S = smeshGen_i->GeomObjectToShape( aGeomObj.in() ); + S = SMESH_Gen_i::GetSMESHGen()->GeomObjectToShape( aGeomObj.in() ); return S; } -//======================================================================= -//function : findShape -//purpose : -//======================================================================= - -static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[], - TopoDS_Shape aShape, - const TopoDS_Shape shape[], - double** box, - const int nShape, - TopAbs_State * state = 0) -{ - gp_XYZ aPnt(0,0,0); - int j, iShape, nbNode = 4; - - for ( j=0; jX(), aNode[j]->Y(), aNode[j]->Z() ); - if ( aNode[j]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE ) { - aPnt = p; - break; - } - aPnt += p / nbNode; - } - - BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion()); - if (state) *state = SC.State(); - if ( SC.State() != TopAbs_IN || aShape.IsNull() || aShape.ShapeType() != TopAbs_SOLID) { - for (iShape = 0; iShape < nShape; iShape++) { - aShape = shape[iShape]; - if ( !( aPnt.X() < box[iShape][0] || box[iShape][1] < aPnt.X() || - aPnt.Y() < box[iShape][2] || box[iShape][3] < aPnt.Y() || - aPnt.Z() < box[iShape][4] || box[iShape][5] < aPnt.Z()) ) { - BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion()); - if (state) *state = SC.State(); - if (SC.State() == TopAbs_IN) - break; - } - } - } - return aShape; -} - -//======================================================================= -//function : readMapIntLine -//purpose : -//======================================================================= - -static char* readMapIntLine(char* ptr, int tab[]) { - long int intVal; - std::cout << std::endl; - - for ( int i=0; i<17; i++ ) { - intVal = strtol(ptr, &ptr, 10); - if ( i < 3 ) - tab[i] = intVal; - } - return ptr; -} - //================================================================================ /*! - * \brief returns true if a triangle defined by the nodes is a temporary face on a - * side facet of pyramid and defines sub-domian inside the pyramid + * \brief returns id of a solid if a triangle defined by the nodes is a temporary face + * either on a side facet of pyramid or a top of pentahedron and defines sub-domian + * outside the volume; else returns HOLE_ID */ //================================================================================ -static bool isTmpFace(const SMDS_MeshNode* node1, - const SMDS_MeshNode* node2, - const SMDS_MeshNode* node3) +static int checkTmpFace(const SMDS_MeshNode* node1, + const SMDS_MeshNode* node2, + const SMDS_MeshNode* node3) { // find a pyramid sharing the 3 nodes - //const SMDS_MeshElement* pyram = 0; SMDS_ElemIteratorPtr vIt1 = node1->GetInverseElementIterator(SMDSAbs_Volume); while ( vIt1->more() ) { - const SMDS_MeshElement* pyram = vIt1->next(); - if ( pyram->NbCornerNodes() != 5 ) continue; + const SMDS_MeshElement* vol = vIt1->next(); + const int nbNodes = vol->NbCornerNodes(); + if ( nbNodes != 5 && nbNodes != 6 ) continue; int i2, i3; - if ( (i2 = pyram->GetNodeIndex( node2 )) >= 0 && - (i3 = pyram->GetNodeIndex( node3 )) >= 0 ) + if ( (i2 = vol->GetNodeIndex( node2 )) >= 0 && + (i3 = vol->GetNodeIndex( node3 )) >= 0 ) { - // Triangle defines sub-domian inside the pyramid if it's - // normal points out of the pyram - - // make i2 and i3 hold indices of base nodes of the pyram while - // keeping the nodes order in the triangle - const int iApex = 4; - if ( i2 == iApex ) - i2 = i3, i3 = pyram->GetNodeIndex( node1 ); - else if ( i3 == iApex ) - i3 = i2, i2 = pyram->GetNodeIndex( node1 ); - - int i3base = (i2+1) % 4; // next index after i2 within the pyramid base - return ( i3base != i3 ); + if ( nbNodes == 5 ) + { + // Triangle defines sub-domian inside the pyramid if it's + // normal points out of the vol + + // make i2 and i3 hold indices of base nodes of the vol while + // keeping the nodes order in the triangle + const int iApex = 4; + if ( i2 == iApex ) + i2 = i3, i3 = vol->GetNodeIndex( node1 ); + else if ( i3 == iApex ) + i3 = i2, i2 = vol->GetNodeIndex( node1 ); + + int i3base = (i2+1) % 4; // next index after i2 within the pyramid base + bool isDomainInPyramid = ( i3base != i3 ); + return isDomainInPyramid ? HOLE_ID : vol->getshapeId(); + } + else + { + return vol->getshapeId(); // triangle is a prism top + } } } - return false; + return HOLE_ID; } //======================================================================= //function : findShapeID -//purpose : find the solid corresponding to GHS3D sub-domain following -// the technique proposed in GHS3D manual (available within -// ghs3d installation) in chapter "B.4 Subdomain (sub-region) assignment". +//purpose : find the solid corresponding to MG-Tetra sub-domain following +// the technique proposed in MG-Tetra manual (available within +// MG-Tetra installation) in chapter "B.4 Subdomain (sub-region) assignment". // In brief: normal of the triangle defined by the given nodes // points out of the domain it is associated to //======================================================================= @@ -340,17 +291,21 @@ static int findShapeID(SMESH_Mesh& mesh, SMESHDS_Mesh* meshDS = mesh.GetMeshDS(); // face the nodes belong to - const SMDS_MeshElement * face = meshDS->FindFace(node1,node2,node3); + vector nodes(3); + nodes[0] = node1; + nodes[1] = node2; + nodes[2] = node3; + const SMDS_MeshElement * face = meshDS->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/true); if ( !face ) - return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID; -#ifdef _DEBUG_ + return checkTmpFace(node1, node2, node3); +#ifdef _MY_DEBUG_ std::cout << "bnd face " << face->GetID() << " - "; #endif // geom face the face assigned to SMESH_MeshEditor editor(&mesh); int geomFaceID = editor.FindShape( face ); if ( !geomFaceID ) - return isTmpFace(node1, node2, node3) ? HOLE_ID : invalidID; + return checkTmpFace(node1, node2, node3); TopoDS_Shape shape = meshDS->IndexToShape( geomFaceID ); if ( shape.IsNull() || shape.ShapeType() != TopAbs_FACE ) return invalidID; @@ -378,25 +333,22 @@ static int findShapeID(SMESH_Mesh& mesh, if ( toMeshHoles ) return meshDS->ShapeToIndex( solid1 ); - //////////// UNCOMMENT AS SOON AS - //////////// http://tracker.dev.opencascade.org/view.php?id=23129 - //////////// IS SOLVED // - Are we at a hole boundary face? - // if ( shells(1).IsSame( BRepTools::OuterShell( solid1 )) ) - // { // - No, but maybe a hole is bound by two shapes? Does shells(1) touches another shell? - // bool touch = false; - // TopExp_Explorer eExp( shells(1), TopAbs_EDGE ); - // // check if any edge of shells(1) belongs to another shell - // for ( ; eExp.More() && !touch; eExp.Next() ) { - // ansIt = mesh.GetAncestors( eExp.Current() ); - // for ( ; ansIt.More() && !touch; ansIt.Next() ) { - // if ( ansIt.Value().ShapeType() == TopAbs_SHELL ) - // touch = ( !ansIt.Value().IsSame( shells(1) )); - // } - // } - // if (!touch) - // return meshDS->ShapeToIndex( solid1 ); - // } + if ( shells(1).IsSame( BRepClass3d::OuterShell( solid1 )) ) + { // - No, but maybe a hole is bound by two shapes? Does shells(1) touch another shell? + bool touch = false; + TopExp_Explorer eExp( shells(1), TopAbs_EDGE ); + // check if any edge of shells(1) belongs to another shell + for ( ; eExp.More() && !touch; eExp.Next() ) { + ansIt = mesh.GetAncestors( eExp.Current() ); + for ( ; ansIt.More() && !touch; ansIt.Next() ) { + if ( ansIt.Value().ShapeType() == TopAbs_SHELL ) + touch = ( !ansIt.Value().IsSame( shells(1) )); + } + } + if (!touch) + return meshDS->ShapeToIndex( solid1 ); + } } // find orientation of geom face within the first solid TopExp_Explorer fExp( solid1, TopAbs_FACE ); @@ -421,7 +373,6 @@ static int findShapeID(SMESH_Mesh& mesh, // get normale to geomFace at any node bool geomNormalOK = false; gp_Vec geomNormal; - const SMDS_MeshNode* nodes[3] = { node1, node2, node3 }; SMESH_MesherHelper helper( mesh ); helper.SetSubShape( geomFace ); for ( int i = 0; !geomNormalOK && i < 3; ++i ) { @@ -433,7 +384,7 @@ static int findShapeID(SMESH_Mesh& mesh, else nNotOnSeamEdge = nodes[(i+1)%3]; } - bool uvOK; + bool uvOK = true; gp_XY uv = helper.GetNodeUV( geomFace, nodes[i], nNotOnSeamEdge, &uvOK ); // check that uv is correct if (uvOK) { @@ -475,521 +426,6 @@ static int findShapeID(SMESH_Mesh& mesh, 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 -// // ======================================================================= -// -// 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) -// { -// 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 nbSubdomains = 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; -// const SMDS_MeshNode** node; -// // std::vector< SMDS_MeshNode* > enfNode( nbRef ); -// SMDS_MeshElement * aGMFElement; -// -// node = new const 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 : addElemInMeshGroup //purpose : Update or create groups in mesh @@ -1013,7 +449,6 @@ static void addElemInMeshGroup(SMESH_Mesh* theMesh, SMESHDS_Group* aGroupDS = static_cast( groupDS ); aGroupDS->SMDSGroup().Add(anElem); groupDone = true; -// MESSAGE("Successfully added enforced element to existing group " << groupName); break; } } @@ -1025,7 +460,6 @@ static void addElemInMeshGroup(SMESH_Mesh* theMesh, aGroup->SetName( groupName.c_str() ); SMESHDS_Group* aGroupDS = static_cast( aGroup->GetGroupDS() ); aGroupDS->SMDSGroup().Add(anElem); -// MESSAGE("Successfully created enforced vertex group " << groupName); groupDone = true; } if (!groupDone) @@ -1049,63 +483,140 @@ static void updateMeshGroups(SMESH_Mesh* theMesh, std::set groupsTo std::string currentGroupName = (string)group->GetName(); if (groupDS->IsEmpty() && groupsToRemove.find(currentGroupName) != groupsToRemove.end()) { // Previous group created by enforced elements - MESSAGE("Delete previous group created by removed enforced elements: " << group->GetName()) theMesh->RemoveGroup(groupDS->GetID()); } } } +//======================================================================= +//function : removeEmptyGroupsOfDomains +//purpose : remove empty groups named "Domain_nb" created due to +// "To make groups of domains" option. +//======================================================================= + +static void removeEmptyGroupsOfDomains(SMESH_Mesh* mesh, + bool notEmptyAsWell = false) +{ + const char* refName = theDomainGroupNamePrefix; + const size_t refLen = strlen( theDomainGroupNamePrefix ); + + std::list groupIDs = mesh->GetGroupIds(); + std::list::const_iterator id = groupIDs.begin(); + for ( ; id != groupIDs.end(); ++id ) + { + SMESH_Group* group = mesh->GetGroup( *id ); + if ( !group || ( !group->GetGroupDS()->IsEmpty() && !notEmptyAsWell )) + continue; + const char* name = group->GetName(); + char* end; + // check the name + if ( strncmp( name, refName, refLen ) == 0 && // starts from refName; + isdigit( *( name + refLen )) && // refName is followed by a digit; + strtol( name + refLen, &end, 10) >= 0 && // there are only digits ... + *end == '\0') // ... till a string end. + { + mesh->RemoveGroup( *id ); + } + } +} + +//================================================================================ +/*! + * \brief Create the groups corresponding to domains + */ +//================================================================================ + +static void makeDomainGroups( std::vector< std::vector< const SMDS_MeshElement* > >& elemsOfDomain, + SMESH_MesherHelper* theHelper) +{ + // int nbDomains = 0; + // for ( size_t i = 0; i < elemsOfDomain.size(); ++i ) + // nbDomains += ( elemsOfDomain[i].size() > 0 ); + + // if ( nbDomains > 1 ) + for ( size_t iDomain = 0; iDomain < elemsOfDomain.size(); ++iDomain ) + { + std::vector< const SMDS_MeshElement* > & elems = elemsOfDomain[ iDomain ]; + if ( elems.empty() ) continue; + + // find existing groups + std::vector< SMESH_Group* > groupOfType( SMDSAbs_NbElementTypes, (SMESH_Group*)NULL ); + const std::string domainName = ( SMESH_Comment( theDomainGroupNamePrefix ) << iDomain ); + SMESH_Mesh::GroupIteratorPtr groupIt = theHelper->GetMesh()->GetGroups(); + while ( groupIt->more() ) + { + SMESH_Group* group = groupIt->next(); + if ( domainName == group->GetName() && + dynamic_cast< SMESHDS_Group* >( group->GetGroupDS()) ) + groupOfType[ group->GetGroupDS()->GetType() ] = group; + } + // create and fill the groups + size_t iElem = 0; + int groupID; + do + { + SMESH_Group* group = groupOfType[ elems[ iElem ]->GetType() ]; + if ( !group ) + group = theHelper->GetMesh()->AddGroup( elems[ iElem ]->GetType(), + domainName.c_str(), groupID ); + SMDS_MeshGroup& groupDS = + static_cast< SMESHDS_Group* >( group->GetGroupDS() )->SMDSGroup(); + + while ( iElem < elems.size() && groupDS.Add( elems[iElem] )) + ++iElem; + + } while ( iElem < elems.size() ); + } +} + //======================================================================= //function : readGMFFile //purpose : read GMF file w/o geometry associated to mesh //======================================================================= -static bool readGMFFile(const char* theFile, -#ifdef WITH_SMESH_CANCEL_COMPUTE +static bool readGMFFile(MG_Tetra_API* MGOutput, + const char* theFile, GHS3DPlugin_GHS3D* theAlgo, -#endif SMESH_MesherHelper* theHelper, - TopoDS_Shape theSolid, - vector & theNodeByGhs3dId, + std::vector & theNodeByGhs3dId, + std::vector & theFaceByGhs3dId, map & theNodeToGhs3dIdMap, std::vector & aNodeGroupByGhs3dId, std::vector & anEdgeGroupByGhs3dId, std::vector & aFaceGroupByGhs3dId, - std::set & groupsToRemove - ) + std::set & groupsToRemove, + bool toMakeGroupsOfDomains=false, + bool toMeshHoles=true) { std::string tmpStr; SMESHDS_Mesh* theMeshDS = theHelper->GetMeshDS(); + const bool hasGeom = ( theHelper->GetMesh()->HasShapeToMesh() ); int nbInitialNodes = theNodeByGhs3dId.size(); - int nbMeshNodes = theMeshDS->NbNodes(); +#ifdef _MY_DEBUG_ const bool isQuadMesh = theHelper->GetMesh()->NbEdges( ORDER_QUADRATIC ) || theHelper->GetMesh()->NbFaces( ORDER_QUADRATIC ) || theHelper->GetMesh()->NbVolumes( ORDER_QUADRATIC ); - -#ifdef _DEBUG_ std::cout << "theNodeByGhs3dId.size(): " << nbInitialNodes << std::endl; - std::cout << "theHelper->GetMesh()->NbNodes(): " << nbMeshNodes << std::endl; + std::cout << "theHelper->GetMesh()->NbNodes(): " << theMeshDS->NbNodes() << std::endl; std::cout << "isQuadMesh: " << isQuadMesh << std::endl; #endif - if (theHelper->GetSubShapeID() != 0) - theHelper->IsQuadraticSubMesh( theHelper->GetSubShape() ); - // --------------------------------- // Read generated elements and nodes // --------------------------------- int nbElem = 0, nbRef = 0; - int aGMFNodeID = 0/*, shapeID*/; - //int *nodeAssigne; - const SMDS_MeshNode** GMFNode; -#ifdef _DEBUG_ + int aGMFNodeID = 0; + std::vector< const SMDS_MeshNode*> GMFNode; +#ifdef _MY_DEBUG_ std::map > subdomainId2tetraId; #endif std::map tabRef; + const bool force3d = !hasGeom; + const int noID = 0; tabRef[GmfVertices] = 3; // for new nodes and enforced nodes tabRef[GmfCorners] = 1; @@ -1117,47 +628,95 @@ static bool readGMFFile(const char* theFile, tabRef[GmfHexahedra] = 8; int ver, dim; - MESSAGE("Read " << theFile << " file"); - int InpMsh = GmfOpenMesh(theFile, GmfRead, &ver, &dim); + int InpMsh = MGOutput->GmfOpenMesh( theFile, GmfRead, &ver, &dim); if (!InpMsh) return false; - MESSAGE("Done "); + + // Read ids of domains + vector< int > solidIDByDomain; + if ( hasGeom ) + { + int solid1; // id used in case of 1 domain or some reading failure + if ( theHelper->GetSubShape().ShapeType() == TopAbs_SOLID ) + solid1 = theHelper->GetSubShapeID(); + else + solid1 = theMeshDS->ShapeToIndex + ( TopExp_Explorer( theHelper->GetSubShape(), TopAbs_SOLID ).Current() ); + + int nbDomains = MGOutput->GmfStatKwd( InpMsh, GmfSubDomainFromGeom ); + if ( nbDomains > 1 ) + { + solidIDByDomain.resize( nbDomains+1, theHelper->GetSubShapeID() ); + int faceNbNodes, faceIndex, orientation, domainNb; + MGOutput->GmfGotoKwd( InpMsh, GmfSubDomainFromGeom ); + for ( int i = 0; i < nbDomains; ++i ) + { + faceIndex = 0; + MGOutput->GmfGetLin( InpMsh, GmfSubDomainFromGeom, + &faceNbNodes, &faceIndex, &orientation, &domainNb, i); + solidIDByDomain[ domainNb ] = 1; + if ( 0 < faceIndex && faceIndex-1 < (int)theFaceByGhs3dId.size() ) + { + const SMDS_MeshElement* face = theFaceByGhs3dId[ faceIndex-1 ]; + const SMDS_MeshNode* nn[3] = { face->GetNode(0), + face->GetNode(1), + face->GetNode(2) }; + if ( orientation < 0 ) + std::swap( nn[1], nn[2] ); + solidIDByDomain[ domainNb ] = + findShapeID( *theHelper->GetMesh(), nn[0], nn[1], nn[2], toMeshHoles ); + if ( solidIDByDomain[ domainNb ] > 0 ) + { +#ifdef _MY_DEBUG_ + std::cout << "solid " << solidIDByDomain[ domainNb ] << std::endl; +#endif + const TopoDS_Shape& foundShape = theMeshDS->IndexToShape( solidIDByDomain[ domainNb ] ); + if ( ! theHelper->IsSubShape( foundShape, theHelper->GetSubShape() )) + solidIDByDomain[ domainNb ] = HOLE_ID; + } + } + } + } + if ( solidIDByDomain.size() < 2 ) + solidIDByDomain.resize( 2, solid1 ); + } // 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_MeshAlgos::GetElementSearcher( *theHelper->GetMeshDS() ); + std::vector< const SMDS_MeshElement* > foundVolumes; + if ( !hasGeom && theHelper->GetMesh()->NbVolumes() > 0 ) + elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *theMeshDS ); + auto_ptr< SMESH_ElementSearcher > elemSearcherDeleter( elemSearcher ); - int nbVertices = GmfStatKwd(InpMsh, GmfVertices) - nbInitialNodes; - GMFNode = new const SMDS_MeshNode*[ nbVertices + 1 ]; - //nodeAssigne = new int[ nbVertices + 1 ]; + // IMP 0022172: [CEA 790] create the groups corresponding to domains + std::vector< std::vector< const SMDS_MeshElement* > > elemsOfDomain; + + int nbVertices = MGOutput->GmfStatKwd( InpMsh, GmfVertices ) - nbInitialNodes; + if ( nbVertices < 0 ) + return false; + GMFNode.resize( 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; + int dummy, solidID; GmfKwdCod token = it->first; - nbRef = it->second; + nbRef = it->second; - nbElem = GmfStatKwd(InpMsh, token); + nbElem = MGOutput->GmfStatKwd( InpMsh, token); if (nbElem > 0) { - GmfGotoKwd(InpMsh, token); + MGOutput->GmfGotoKwd( InpMsh, token); std::cout << "Read " << nbElem; } else continue; std::vector id (nbElem*tabRef[token]); // node ids + std::vector domainID( nbElem ); // domain if (token == GmfVertices) { (nbElem <= 1) ? tmpStr = " vertex" : tmpStr = " vertices"; @@ -1182,24 +741,18 @@ static bool readGMFFile(const char* theFile, double x, y, z; const SMDS_MeshNode * aGMFNode; - //shapeID = theMeshDS->ShapeToIndex( theSolid ); for ( int iElem = 0; iElem < nbElem; iElem++ ) { -#ifdef WITH_SMESH_CANCEL_COMPUTE if(theAlgo->computeCanceled()) { - GmfCloseMesh(InpMsh); - delete [] GMFNode; - //delete [] nodeAssigne; return false; } -#endif if (ver == GmfFloat) { - GmfGetLin(InpMsh, token, &VerTab_f[0], &VerTab_f[1], &VerTab_f[2], &dummy); + MGOutput->GmfGetLin( InpMsh, token, &VerTab_f[0], &VerTab_f[1], &VerTab_f[2], &dummy); x = VerTab_f[0]; y = VerTab_f[1]; z = VerTab_f[2]; } else { - GmfGetLin(InpMsh, token, &x, &y, &z, &dummy); + MGOutput->GmfGetLin( InpMsh, token, &x, &y, &z, &dummy); } if (iElem >= nbInitialNodes) { if ( elemSearcher && @@ -1210,8 +763,7 @@ static bool readGMFFile(const char* theFile, aGMFID = iElem -nbInitialNodes +1; GMFNode[ aGMFID ] = aGMFNode; - //nodeAssigne[ aGMFID ] = 0; - if (aGMFID-1 < aNodeGroupByGhs3dId.size() && !aNodeGroupByGhs3dId.at(aGMFID-1).empty()) + if (aGMFID-1 < (int)aNodeGroupByGhs3dId.size() && !aNodeGroupByGhs3dId.at(aGMFID-1).empty()) addElemInMeshGroup(theHelper->GetMesh(), aGMFNode, aNodeGroupByGhs3dId.at(aGMFID-1), groupsToRemove); } } @@ -1219,43 +771,42 @@ static bool readGMFFile(const char* theFile, else if (token == GmfCorners && nbElem > 0) { (nbElem <= 1) ? tmpStr = " corner" : tmpStr = " corners"; for ( int iElem = 0; iElem < nbElem; iElem++ ) - GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]); + MGOutput->GmfGetLin( InpMsh, token, &id[iElem*tabRef[token]]); } else if (token == GmfRidges && nbElem > 0) { (nbElem <= 1) ? tmpStr = " ridge" : tmpStr = " ridges"; for ( int iElem = 0; iElem < nbElem; iElem++ ) - GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]]); + MGOutput->GmfGetLin( InpMsh, token, &id[iElem*tabRef[token]]); } else if (token == GmfEdges && nbElem > 0) { (nbElem <= 1) ? tmpStr = " edge" : tmpStr = " edges"; for ( int iElem = 0; iElem < nbElem; iElem++ ) - GmfGetLin(InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &dummy); + MGOutput->GmfGetLin( InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &domainID[iElem]); } else if (token == GmfTriangles && nbElem > 0) { (nbElem <= 1) ? tmpStr = " triangle" : tmpStr = " triangles"; 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); + MGOutput->GmfGetLin( InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &domainID[iElem]); } else if (token == GmfQuadrilaterals && nbElem > 0) { (nbElem <= 1) ? tmpStr = " Quadrilateral" : tmpStr = " Quadrilaterals"; 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); + MGOutput->GmfGetLin( InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &domainID[iElem]); } else if (token == GmfTetrahedra && nbElem > 0) { (nbElem <= 1) ? tmpStr = " Tetrahedron" : tmpStr = " Tetrahedra"; 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); -#ifdef _DEBUG_ + MGOutput->GmfGetLin( InpMsh, token, &id[iElem*tabRef[token]], &id[iElem*tabRef[token]+1], &id[iElem*tabRef[token]+2], &id[iElem*tabRef[token]+3], &domainID[iElem]); +#ifdef _MY_DEBUG_ subdomainId2tetraId[dummy].insert(iElem+1); -// MESSAGE("subdomainId2tetraId["< 0) { (nbElem <= 1) ? tmpStr = " Hexahedron" : tmpStr = " Hexahedra"; 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); + MGOutput->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], &domainID[iElem]); } std::cout << tmpStr << std::endl; std::cout << std::endl; @@ -1276,14 +827,9 @@ static bool readGMFFile(const char* theFile, for ( int iElem = 0; iElem < nbElem; iElem++ ) { -#ifdef WITH_SMESH_CANCEL_COMPUTE if(theAlgo->computeCanceled()) { - GmfCloseMesh(InpMsh); - delete [] GMFNode; - //delete [] nodeAssigne; return false; } -#endif // Check if elem is already in input mesh. If yes => skip bool fullyCreatedElement = false; // if at least one of the nodes was created for ( int iRef = 0; iRef < nbRef; iRef++ ) @@ -1302,87 +848,129 @@ static bool readGMFFile(const char* theFile, node [ iRef ] = GMFNode[ aGMFNodeID ]; } } - + aCreatedElem = 0; switch (token) { case GmfEdges: if (fullyCreatedElement) { - aCreatedElem = theHelper->AddEdge( node[0], node[1], /*id =*/0, /*force3d =*/false ); + aCreatedElem = theHelper->AddEdge( node[0], node[1], noID, force3d ); if (anEdgeGroupByGhs3dId.size() && !anEdgeGroupByGhs3dId[iElem].empty()) addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, anEdgeGroupByGhs3dId[iElem], groupsToRemove); } break; case GmfTriangles: if (fullyCreatedElement) { - aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], /*id =*/0, /*force3d =*/false ); - // for ( int iRef = 0; iRef < nbRef; iRef++ ) - // nodeAssigne[ nodeID[ iRef ]] = 1; + aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], noID, force3d ); if (aFaceGroupByGhs3dId.size() && !aFaceGroupByGhs3dId[iElem].empty()) addElemInMeshGroup(theHelper->GetMesh(), aCreatedElem, aFaceGroupByGhs3dId[iElem], groupsToRemove); } break; case GmfQuadrilaterals: if (fullyCreatedElement) { - theHelper->AddFace( node[0], node[1], node[2], node[3], /*id =*/0, /*force3d =*/false ); - // for ( int iRef = 0; iRef < nbRef; iRef++ ) - // nodeAssigne[ nodeID[ iRef ]] = 1; + aCreatedElem = theHelper->AddFace( node[0], node[1], node[2], node[3], noID, force3d ); } break; case GmfTetrahedra: - 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 )) - break; + if ( hasGeom ) + { + solidID = solidIDByDomain[ domainID[iElem]]; + if ( solidID != HOLE_ID ) + { + aCreatedElem = theHelper->AddVolume( node[1], node[0], node[2], node[3], + noID, force3d ); + theMeshDS->SetMeshElementOnShape( aCreatedElem, solidID ); + for ( int iN = 0; iN < 4; ++iN ) + if ( node[iN]->getshapeId() < 1 ) + theMeshDS->SetNodeInVolume( node[iN], solidID ); + } + } + else + { + 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 )) + break; + } + aCreatedElem = theHelper->AddVolume( node[1], node[0], node[2], node[3], + noID, force3d ); } - theHelper->AddVolume( node[1], node[0], node[2], node[3], /*id =*/0, /*force3d =*/false ); -// theMeshDS->SetMeshElementOnShape( aTet, shapeID ); break; case GmfHexahedra: - 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] || !node[4] || !node[5] || !node[6] || !node[7]) - continue; - if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) + - SMESH_TNodeXYZ(node[1]) + - SMESH_TNodeXYZ(node[2]) + - SMESH_TNodeXYZ(node[3]) + - SMESH_TNodeXYZ(node[4]) + - SMESH_TNodeXYZ(node[5]) + - SMESH_TNodeXYZ(node[6]) + - SMESH_TNodeXYZ(node[7])) / 8., - SMDSAbs_Volume, foundVolumes )) - break; + if ( hasGeom ) + { + solidID = solidIDByDomain[ domainID[iElem]]; + if ( solidID != HOLE_ID ) + { + aCreatedElem = theHelper->AddVolume( node[0], node[3], node[2], node[1], + node[4], node[7], node[6], node[5], + noID, force3d ); + theMeshDS->SetMeshElementOnShape( aCreatedElem, solidID ); + for ( int iN = 0; iN < 8; ++iN ) + if ( node[iN]->getshapeId() < 1 ) + theMeshDS->SetNodeInVolume( node[iN], solidID ); + } + } + else + { + 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] || !node[4] || !node[5] || !node[6] || !node[7]) + continue; + if ( elemSearcher->FindElementsByPoint((SMESH_TNodeXYZ(node[0]) + + SMESH_TNodeXYZ(node[1]) + + SMESH_TNodeXYZ(node[2]) + + SMESH_TNodeXYZ(node[3]) + + SMESH_TNodeXYZ(node[4]) + + SMESH_TNodeXYZ(node[5]) + + SMESH_TNodeXYZ(node[6]) + + SMESH_TNodeXYZ(node[7])) / 8., + SMDSAbs_Volume, foundVolumes )) + break; + } + aCreatedElem = theHelper->AddVolume( node[0], node[3], node[2], node[1], + node[4], node[7], node[6], node[5], + noID, force3d ); } - theHelper->AddVolume( node[0], node[3], node[2], node[1], - node[4], node[7], node[6], node[5], /*id =*/0, /*force3d =*/false ); -// theMeshDS->SetMeshElementOnShape( aTet, shapeID ); break; default: continue; + } // switch (token) + + if ( aCreatedElem && toMakeGroupsOfDomains ) + { + if ( domainID[iElem] >= (int) elemsOfDomain.size() ) + elemsOfDomain.resize( domainID[iElem] + 1 ); + elemsOfDomain[ domainID[iElem] ].push_back( aCreatedElem ); } - } + } // loop on elements of one type break; - } - } + } // case ... + default:; + } // switch (token) + } // loop on tabRef + + // remove nodes in holes + if ( hasGeom ) + { + for ( int i = 1; i <= nbVertices; ++i ) + if ( GMFNode[i]->NbInverseElements() == 0 ) + theMeshDS->RemoveFreeNode( GMFNode[i], /*sm=*/0, /*fromGroups=*/false ); } - // for ( int i = 0; i < nbVertices; ++i ) { - // if ( !nodeAssigne[ i+1 ]) - // theMeshDS->SetNodeInVolume( GMFNode[ i+1 ], shapeID ); - // } + MGOutput->GmfCloseMesh( InpMsh); - GmfCloseMesh(InpMsh); - delete [] GMFNode; - //delete [] nodeAssigne; -#ifdef _DEBUG_ - MESSAGE("Nb subdomains " << subdomainId2tetraId.size()); + // 0022172: [CEA 790] create the groups corresponding to domains + if ( toMakeGroupsOfDomains ) + makeDomainGroups( elemsOfDomain, theHelper ); + +#ifdef _MY_DEBUG_ std::map >::const_iterator subdomainIt = subdomainId2tetraId.begin(); TCollection_AsciiString aSubdomainFileName = theFile; aSubdomainFileName = aSubdomainFileName + ".subdomain"; @@ -1392,7 +980,6 @@ static bool readGMFFile(const char* theFile, for(;subdomainIt != subdomainId2tetraId.end() ; ++subdomainIt) { int subdomainId = subdomainIt->first; std::set tetraIds = subdomainIt->second; - MESSAGE("Subdomain #"<::const_iterator tetraIdsIt = tetraIds.begin(); aSubdomainFile << subdomainId << std::endl; for(;tetraIdsIt != tetraIds.end() ; ++tetraIdsIt) { @@ -1406,12 +993,15 @@ static bool readGMFFile(const char* theFile, return true; } -static bool writeGMFFile(const char* theMeshFileName, + +static bool writeGMFFile(MG_Tetra_API* MGInput, + const char* theMeshFileName, const char* theRequiredFileName, const char* theSolFileName, const SMESH_ProxyMesh& theProxyMesh, - SMESH_Mesh * theMesh, + SMESH_MesherHelper& theHelper, std::vector & theNodeByGhs3dId, + std::vector & theFaceByGhs3dId, std::map & aNodeToGhs3dIdMap, std::vector & aNodeGroupByGhs3dId, std::vector & anEdgeGroupByGhs3dId, @@ -1420,9 +1010,9 @@ static bool writeGMFFile(const char* theMesh GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges, GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles, std::map, std::string> & enfVerticesWithGroup, - GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices) + GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices, + int & theInvalidEnforcedFlags) { - MESSAGE("writeGMFFile w/o geometry"); std::string tmpStr; int idx, idxRequired = 0, idxSol = 0; const int dummyint = 0; @@ -1439,14 +1029,18 @@ static bool writeGMFFile(const char* theMesh GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap::iterator elemIt; TIDSortedElemSet::iterator elemSetIt; bool isOK; - auto_ptr< SMESH_ElementSearcher > pntCls + SMESH_Mesh* theMesh = theHelper.GetMesh(); + const bool hasGeom = theMesh->HasShapeToMesh(); + SMESHUtils::Deleter< SMESH_ElementSearcher > pntCls ( SMESH_MeshAlgos::GetElementSearcher(*theMesh->GetMeshDS())); int nbEnforcedVertices = theEnforcedVertices.size(); - + theInvalidEnforcedFlags = 0; + // count faces int nbFaces = theProxyMesh.NbFaces(); int nbNodes; + theFaceByGhs3dId.reserve( nbFaces ); // groups management int usedEnforcedNodes = 0; @@ -1455,13 +1049,14 @@ static bool writeGMFFile(const char* theMesh if ( nbFaces == 0 ) return false; - idx = GmfOpenMesh(theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + idx = MGInput->GmfOpenMesh( theMeshFileName, GmfWrite, GMFVERSION, GMFDIMENSION); if (!idx) return false; /* ========================== FACES ========================== */ /* TRIANGLES ========================== */ - SMDS_ElemIteratorPtr eIt = theProxyMesh.GetFaces(); + SMDS_ElemIteratorPtr eIt = + hasGeom ? theProxyMesh.GetFaces( theHelper.GetSubShape()) : theProxyMesh.GetFaces(); while ( eIt->more() ) { elem = eIt->next(); @@ -1470,9 +1065,9 @@ static bool writeGMFFile(const char* theMesh nbNodes = elem->NbCornerNodes(); while ( nodeIt->more() && nbNodes--) { - // find GHS3D ID + // find MG-Tetra ID const SMDS_MeshNode* node = castToNode( nodeIt->next() ); - int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1 + int newId = aNodeToGhs3dIdMap.size() + 1; // MG-Tetra ids count from 1 aNodeToGhs3dIdMap.insert( make_pair( node, newId )); } } @@ -1486,13 +1081,14 @@ static bool writeGMFFile(const char* theMesh nodeIt = elem->nodesIterator(); nbNodes = 2; while ( nodeIt->more() && nbNodes-- ) { - // find GHS3D ID + // find MG-Tetra 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_OUT ) { isOK = false; + theInvalidEnforcedFlags |= FLAG_BAD_ENF_EDGE; break; } aNodeToTopAbs_StateMap.insert( make_pair( node, result )); @@ -1502,17 +1098,17 @@ static bool writeGMFFile(const char* theMesh nbNodes = 2; int newId = -1; while ( nodeIt->more() && nbNodes-- ) { - // find GHS3D ID + // find MG-Tetra ID const SMDS_MeshNode* node = castToNode( nodeIt->next() ); gp_Pnt myPoint(node->X(),node->Y(),node->Z()); nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems); -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << "Node at "<X()<<", "<Y()<<", "<Z()<nodesIterator(); nbNodes = 3; while ( nodeIt->more() && nbNodes--) { - // find GHS3D ID + // find MG-Tetra 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_OUT ) { isOK = false; + theInvalidEnforcedFlags |= FLAG_BAD_ENF_TRIA; break; } aNodeToTopAbs_StateMap.insert( make_pair( node, result )); @@ -1557,16 +1156,16 @@ static bool writeGMFFile(const char* theMesh nbNodes = 3; int newId = -1; while ( nodeIt->more() && nbNodes--) { - // find GHS3D ID + // find MG-Tetra ID const SMDS_MeshNode* node = castToNode( nodeIt->next() ); gp_Pnt myPoint(node->X(),node->Y(),node->Z()); nbFoundElems = pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems); -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << "Nb nodes found : "<first: "<first<second - 1 ] = n2id->first; // ghs3d ids count from 1 + theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // MG-Tetra ids count from 1 } // put nodes to anEnforcedNodeToGhs3dIdMap vector -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << "anEnforcedNodeToGhs3dIdMap.size(): "<second > aNodeToGhs3dIdMap.size()) { - theEnforcedNodeByGhs3dId[ n2id->second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // ghs3d ids count from 1 + if (n2id->second > (int)aNodeToGhs3dIdMap.size()) { + theEnforcedNodeByGhs3dId[ n2id->second - aNodeToGhs3dIdMap.size() - 1 ] = n2id->first; // MG-Tetra ids count from 1 } } @@ -1643,13 +1244,13 @@ static bool writeGMFFile(const char* theMesh coords.push_back(node->X()); coords.push_back(node->Y()); coords.push_back(node->Z()); -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << "Node at " << node->X()<<", " <Y()<<", " <Z(); #endif if (nodesCoords.find(coords) != nodesCoords.end()) { // node already exists in original mesh -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " found" << std::endl; #endif continue; @@ -1657,7 +1258,7 @@ static bool writeGMFFile(const char* theMesh if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) { // node already exists in enforced vertices -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " found" << std::endl; #endif continue; @@ -1679,7 +1280,7 @@ static bool writeGMFFile(const char* theMesh // theOrderedNodes.push_back(existingNode); // } -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " not found" << std::endl; #endif @@ -1700,7 +1301,7 @@ static bool writeGMFFile(const char* theMesh coords.push_back(node->X()); coords.push_back(node->Y()); coords.push_back(node->Z()); -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << "Node at " << node->X()<<", " <Y()<<", " <Z(); #endif @@ -1708,14 +1309,15 @@ static bool writeGMFFile(const char* theMesh gp_Pnt myPoint(node->X(),node->Y(),node->Z()); TopAbs_State result = pntCls->GetPointState( myPoint ); if ( result == TopAbs_OUT ) { -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " out of volume" << std::endl; #endif + theInvalidEnforcedFlags |= FLAG_BAD_ENF_NODE; continue; } if (nodesCoords.find(coords) != nodesCoords.end()) { -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " found in nodesCoords" << std::endl; #endif // theRequiredNodes.push_back(node); @@ -1723,7 +1325,7 @@ static bool writeGMFFile(const char* theMesh } if (theEnforcedVertices.find(coords) != theEnforcedVertices.end()) { -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " found in theEnforcedVertices" << std::endl; #endif continue; @@ -1752,7 +1354,7 @@ static bool writeGMFFile(const char* theMesh // if ( result != TopAbs_IN ) // continue; -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << " not found" << std::endl; #endif nodesCoords.insert(coords); @@ -1776,12 +1378,11 @@ static bool writeGMFFile(const char* theMesh gp_Pnt myPoint(x,y,z); TopAbs_State result = pntCls->GetPointState( myPoint ); if ( result == TopAbs_OUT ) - continue; - //if (pntCls->FindElementsByPoint(myPoint, SMDSAbs_Node, foundElems) == 0) - //continue; - -// if ( result != TopAbs_IN ) -// continue; + { + std::cout << "Warning: enforced vertex at ( " << x << "," << y << "," << z << " ) is out of the meshed domain!!!" << std::endl; + theInvalidEnforcedFlags |= FLAG_BAD_ENF_VERT; + //continue; + } std::vector coords; coords.push_back(x); coords.push_back(y); @@ -1796,9 +1397,9 @@ static bool writeGMFFile(const char* theMesh // GmfVertices std::cout << "Begin writting required nodes in GmfVertices" << std::endl; std::cout << "Nb vertices: " << theOrderedNodes.size() << std::endl; - GmfSetKwd(idx, GmfVertices, theOrderedNodes.size()/*+solSize*/); + MGInput->GmfSetKwd( idx, GmfVertices, theOrderedNodes.size()/*+solSize*/); for (ghs3dNodeIt = theOrderedNodes.begin();ghs3dNodeIt != theOrderedNodes.end();++ghs3dNodeIt) { - GmfSetLin(idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint); + MGInput->GmfSetLin( idx, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint); } std::cout << "End writting required nodes in GmfVertices" << std::endl; @@ -1806,27 +1407,23 @@ static bool writeGMFFile(const char* theMesh if (requiredNodes + solSize) { std::cout << "Begin writting in req and sol file" << std::endl; aNodeGroupByGhs3dId.resize( requiredNodes + solSize ); - idxRequired = GmfOpenMesh(theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + idxRequired = MGInput->GmfOpenMesh( theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION); if (!idxRequired) { - GmfCloseMesh(idx); return false; } - idxSol = GmfOpenMesh(theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION); + idxSol = MGInput->GmfOpenMesh( theSolFileName, GmfWrite, GMFVERSION, GMFDIMENSION); if (!idxSol) { - GmfCloseMesh(idx); - if (idxRequired) - GmfCloseMesh(idxRequired); return false; } int TypTab[] = {GmfSca}; double ValTab[] = {0.0}; - GmfSetKwd(idxRequired, GmfVertices, requiredNodes + solSize); - GmfSetKwd(idxSol, GmfSolAtVertices, requiredNodes + solSize, 1, TypTab); + MGInput->GmfSetKwd( idxRequired, GmfVertices, requiredNodes + solSize); + MGInput->GmfSetKwd( idxSol, GmfSolAtVertices, requiredNodes + solSize, 1, TypTab); // int usedEnforcedNodes = 0; // std::string gn = ""; for (ghs3dNodeIt = theRequiredNodes.begin();ghs3dNodeIt != theRequiredNodes.end();++ghs3dNodeIt) { - GmfSetLin(idxRequired, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint); - GmfSetLin(idxSol, GmfSolAtVertices, ValTab); + MGInput->GmfSetLin( idxRequired, GmfVertices, (*ghs3dNodeIt)->X(), (*ghs3dNodeIt)->Y(), (*ghs3dNodeIt)->Z(), dummyint); + MGInput->GmfSetLin( idxSol, GmfSolAtVertices, ValTab); if (theEnforcedNodes.find((*ghs3dNodeIt)) != theEnforcedNodes.end()) gn = theEnforcedNodes.find((*ghs3dNodeIt))->second; aNodeGroupByGhs3dId[usedEnforcedNodes] = gn; @@ -1835,14 +1432,14 @@ static bool writeGMFFile(const char* theMesh for (int i=0;iGmfSetLin( idxRequired, GmfVertices, ReqVerTab[i][0], ReqVerTab[i][1], ReqVerTab[i][2], dummyint); + MGInput->GmfSetLin( idxSol, GmfSolAtVertices, solTab); aNodeGroupByGhs3dId[usedEnforcedNodes] = enfVerticesWithGroup.find(ReqVerTab[i])->second; -#ifdef _DEBUG_ +#ifdef _MY_DEBUG_ std::cout << "aNodeGroupByGhs3dId["<GmfOpenMesh( theRequiredFileName, GmfWrite, GMFVERSION, GMFDIMENSION); // if (!idxRequired) // return false; - GmfSetKwd(idx, GmfEdges, theKeptEnforcedEdges.size()); -// GmfSetKwd(idxRequired, GmfEdges, theKeptEnforcedEdges.size()); + MGInput->GmfSetKwd( idx, GmfEdges, theKeptEnforcedEdges.size()); +// MGInput->GmfSetKwd( idxRequired, GmfEdges, theKeptEnforcedEdges.size()); for(elemSetIt = theKeptEnforcedEdges.begin() ; elemSetIt != theKeptEnforcedEdges.end() ; ++elemSetIt) { elem = (*elemSetIt); nodeIt = elem->nodesIterator(); int index=0; while ( nodeIt->more() ) { - // find GHS3D ID + // find MG-Tetra ID const SMDS_MeshNode* node = castToNode( nodeIt->next() ); map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node); if (it == anEnforcedNodeToGhs3dIdMap.end()) { @@ -1877,19 +1474,18 @@ static bool writeGMFFile(const char* theMesh nedge[index] = it->second; index++; } - GmfSetLin(idx, GmfEdges, nedge[0], nedge[1], dummyint); + MGInput->GmfSetLin( idx, GmfEdges, nedge[0], nedge[1], dummyint); anEdgeGroupByGhs3dId[usedEnforcedEdges] = theEnforcedEdges.find(elem)->second; -// GmfSetLin(idxRequired, GmfEdges, nedge[0], nedge[1], dummyint); +// MGInput->GmfSetLin( idxRequired, GmfEdges, nedge[0], nedge[1], dummyint); usedEnforcedEdges++; } -// GmfCloseMesh(idxRequired); } if (usedEnforcedEdges) { - GmfSetKwd(idx, GmfRequiredEdges, usedEnforcedEdges); + MGInput->GmfSetKwd( idx, GmfRequiredEdges, usedEnforcedEdges); for (int enfID=1;enfID<=usedEnforcedEdges;enfID++) { - GmfSetLin(idx, GmfRequiredEdges, enfID); + MGInput->GmfSetLin( idx, GmfRequiredEdges, enfID); } } @@ -1897,14 +1493,15 @@ static bool writeGMFFile(const char* theMesh int usedEnforcedTriangles = 0; if (anElemSet.size()+theKeptEnforcedTriangles.size()) { aFaceGroupByGhs3dId.resize( anElemSet.size()+theKeptEnforcedTriangles.size() ); - GmfSetKwd(idx, GmfTriangles, anElemSet.size()+theKeptEnforcedTriangles.size()); + MGInput->GmfSetKwd( idx, GmfTriangles, anElemSet.size()+theKeptEnforcedTriangles.size()); int k=0; for(elemSetIt = anElemSet.begin() ; elemSetIt != anElemSet.end() ; ++elemSetIt,++k) { elem = (*elemSetIt); + theFaceByGhs3dId.push_back( elem ); nodeIt = elem->nodesIterator(); int index=0; for ( int j = 0; j < 3; ++j ) { - // find GHS3D ID + // find MG-Tetra ID const SMDS_MeshNode* node = castToNode( nodeIt->next() ); map< const SMDS_MeshNode*,int >::iterator it = aNodeToGhs3dIdMap.find(node); if (it == aNodeToGhs3dIdMap.end()) @@ -1912,16 +1509,18 @@ static bool writeGMFFile(const char* theMesh ntri[index] = it->second; index++; } - GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint); + MGInput->GmfSetLin( idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint); aFaceGroupByGhs3dId[k] = ""; } + if ( !theHelper.GetMesh()->HasShapeToMesh() ) + SMESHUtils::FreeVector( theFaceByGhs3dId ); if (theKeptEnforcedTriangles.size()) { for(elemSetIt = theKeptEnforcedTriangles.begin() ; elemSetIt != theKeptEnforcedTriangles.end() ; ++elemSetIt,++k) { elem = (*elemSetIt); nodeIt = elem->nodesIterator(); int index=0; for ( int j = 0; j < 3; ++j ) { - // find GHS3D ID + // find MG-Tetra ID const SMDS_MeshNode* node = castToNode( nodeIt->next() ); map< const SMDS_MeshNode*,int >::iterator it = anEnforcedNodeToGhs3dIdMap.find(node); if (it == anEnforcedNodeToGhs3dIdMap.end()) { @@ -1932,7 +1531,7 @@ static bool writeGMFFile(const char* theMesh ntri[index] = it->second; index++; } - GmfSetLin(idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint); + MGInput->GmfSetLin( idx, GmfTriangles, ntri[0], ntri[1], ntri[2], dummyint); aFaceGroupByGhs3dId[k] = theEnforcedTriangles.find(elem)->second; usedEnforcedTriangles++; } @@ -1941,1230 +1540,23 @@ static bool writeGMFFile(const char* theMesh if (usedEnforcedTriangles) { - GmfSetKwd(idx, GmfRequiredTriangles, usedEnforcedTriangles); + MGInput->GmfSetKwd( idx, GmfRequiredTriangles, usedEnforcedTriangles); for (int enfID=1;enfID<=usedEnforcedTriangles;enfID++) - GmfSetLin(idx, GmfRequiredTriangles, anElemSet.size()+enfID); + MGInput->GmfSetLin( idx, GmfRequiredTriangles, anElemSet.size()+enfID); } - GmfCloseMesh(idx); + MGInput->GmfCloseMesh(idx); if (idxRequired) - GmfCloseMesh(idxRequired); + MGInput->GmfCloseMesh(idxRequired); if (idxSol) - GmfCloseMesh(idxSol); - + MGInput->GmfCloseMesh(idxSol); + return true; - } -// static bool writeGMFFile(const char* theMeshFileName, -// const char* theRequiredFileName, -// const char* theSolFileName, -// SMESH_MesherHelper& theHelper, -// const SMESH_ProxyMesh& theProxyMesh, -// std::map & theNodeId2NodeIndexMap, -// std::map & theSmdsToGhs3dIdMap, -// std::map & theGhs3dIdToNodeMap, -// TIDSortedNodeSet & theEnforcedNodes, -// TIDSortedElemSet & theEnforcedEdges, -// TIDSortedElemSet & theEnforcedTriangles, -// // TIDSortedElemSet & theEnforcedQuadrangles, -// GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices) -// { -// MESSAGE("writeGMFFile with geometry"); -// int idx, idxRequired, idxSol; -// int nbv, nbev, nben, aGhs3dID = 0; -// const int dummyint = 0; -// GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt; -// std::vector enfVertexSizes; -// TIDSortedNodeSet::const_iterator enfNodeIt; -// const SMDS_MeshNode* node; -// SMDS_NodeIteratorPtr nodeIt; -// -// 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(); -// -// // 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 nodes which are free -// // and replace not-free nodes on edges by the node on vertex -// TNodeNodeMap n2nDegen; // map a node on degenerated edge to a node on vertex -// TNodeNodeMap::iterator n2nDegenIt; -// if ( theHelper.HasDegeneratedEdges() ) -// { -// 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 ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() ) -// { -// TopoDS_Shape vertex = TopoDS_Iterator( e.Current() ).Value(); -// const SMDS_MeshNode* vNode = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), theMeshDS); -// { -// SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); -// while ( nIt->more() ) -// n2nDegen.insert( make_pair( nIt->next(), vNode )); -// } -// } -// } -// } -// } -// -// const bool isQuadMesh = -// theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) || -// theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) || -// theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC ); -// -// std::vector > VerTab; -// std::set > VerMap; -// VerTab.clear(); -// std::vector aVerTab; -// // Loop from 1 to NB_NODES -// -// nodeIt = theMeshDS->nodesIterator(); -// -// while ( nodeIt->more() ) -// { -// node = nodeIt->next(); -// if ( isQuadMesh && theHelper.IsMedium( node )) // Issue 0021238 -// continue; -// if ( n2nDegen.count( node ) ) // Issue 0020674 -// continue; -// -// std::vector coords; -// coords.push_back(node->X()); -// coords.push_back(node->Y()); -// coords.push_back(node->Z()); -// if (VerMap.find(coords) != VerMap.end()) { -// aGhs3dID = theSmdsToGhs3dIdMap[node->GetID()]; -// theGhs3dIdToNodeMap[theSmdsToGhs3dIdMap[node->GetID()]] = node; -// continue; -// } -// VerTab.push_back(coords); -// VerMap.insert(coords); -// aGhs3dID++; -// theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID )); -// theGhs3dIdToNodeMap.insert( make_pair( aGhs3dID, node )); -// } -// -// -// /* ENFORCED NODES ========================== */ -// if (nben) { -// std::cout << "Add " << nben << " enforced nodes to input .mesh file" << std::endl; -// 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; -// if (VerMap.find(coords) != VerMap.end()) -// continue; -// VerTab.push_back(coords); -// VerMap.insert(coords); -// aGhs3dID++; -// theNodeId2NodeIndexMap.insert( make_pair( (*enfNodeIt)->GetID(), aGhs3dID )); -// } -// } -// -// -// /* ENFORCED VERTICES ========================== */ -// int solSize = 0; -// std::vector > ReqVerTab; -// ReqVerTab.clear(); -// if (nbev) { -// std::cout << "Add " << nbev << " enforced vertices to input .mesh file" << std::endl; -// 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); -// std::vector coords; -// coords.push_back(x); -// coords.push_back(y); -// coords.push_back(z); -// if (VerMap.find(coords) != VerMap.end()) -// continue; -// ReqVerTab.push_back(coords); -// VerMap.insert(coords); -// solSize++; -// } -// } -// -// -// /* ========================== FACES ========================== */ -// -// int nbTriangles = 0/*, nbQuadrangles = 0*/, aSmdsID; -// TopTools_IndexedMapOfShape facesMap, trianglesMap/*, quadranglesMap*/; -// TIDSortedElemSet::const_iterator elemIt; -// const SMESHDS_SubMesh* theSubMesh; -// TopoDS_Shape aShape; -// SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace; -// const SMDS_MeshElement* aFace; -// map::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(); -// int nbCornerNodes = elem->NbCornerNodes(); -// if (nbCornerNodes == 3) -// { -// trianglesMap.Add(facesMap(i)); -// nbTriangles ++; -// } -// // else if (nbCornerNodes == 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(); -// for ( int j = 0; j < 3; ++j ) { -// // find GHS3D ID -// node = castToNode( itOnSubFace->next() ); -// if (( n2nDegenIt = n2nDegen.find( node )) != n2nDegen.end() ) -// node = n2nDegenIt->second; -// aSmdsID = node->GetID(); -// itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID ); -// ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() ); -// att.push_back((*itOnMap).second); -// } -// tt.push_back(att); -// } -// } -// } -// -// if (theEnforcedTriangles.size()) { -// std::cout << "Add " << theEnforcedTriangles.size() << " enforced triangles to input .mesh file" << std::endl; -// // Iterate over the enforced triangles -// for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) { -// aFace = (*elemIt); -// itOnSubFace = aFace->nodesIterator(); -// bool isOK = true; -// att.clear(); -// -// for ( int j = 0; j < 3; ++j ) { -// node = castToNode( itOnSubFace->next() ); -// if (( n2nDegenIt = n2nDegen.find( node )) != n2nDegen.end() ) -// node = n2nDegenIt->second; -// // std::cout << node; -// double x = node->X(); -// double y = node->Y(); -// double z = node->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 ) { -// isOK = false; -// theEnforcedTriangles.erase(elemIt); -// continue; -// } -// std::vector coords; -// coords.push_back(x); -// coords.push_back(y); -// coords.push_back(z); -// if (VerMap.find(coords) != VerMap.end()) { -// att.push_back(theNodeId2NodeIndexMap[node->GetID()]); -// continue; -// } -// VerTab.push_back(coords); -// VerMap.insert(coords); -// aGhs3dID++; -// theNodeId2NodeIndexMap.insert( make_pair( node->GetID(), aGhs3dID )); -// att.push_back(aGhs3dID); -// } -// if (isOK) -// tt.push_back(att); -// } -// } -// -// -// /* ========================== EDGES ========================== */ -// -// if (theEnforcedEdges.size()) { -// // Iterate over the enforced edges -// std::cout << "Add " << theEnforcedEdges.size() << " enforced edges to input .mesh file" << std::endl; -// for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) { -// aFace = (*elemIt); -// bool isOK = true; -// itOnSubFace = aFace->nodesIterator(); -// aet.clear(); -// for ( int j = 0; j < 2; ++j ) { -// node = castToNode( itOnSubFace->next() ); -// if (( n2nDegenIt = n2nDegen.find( node )) != n2nDegen.end() ) -// node = n2nDegenIt->second; -// double x = node->X(); -// double y = node->Y(); -// double z = node->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 ) { -// isOK = false; -// theEnforcedEdges.erase(elemIt); -// continue; -// } -// std::vector coords; -// coords.push_back(x); -// coords.push_back(y); -// coords.push_back(z); -// if (VerMap.find(coords) != VerMap.end()) { -// aet.push_back(theNodeId2NodeIndexMap[node->GetID()]); -// continue; -// } -// VerTab.push_back(coords); -// VerMap.insert(coords); -// -// aGhs3dID++; -// theNodeId2NodeIndexMap.insert( make_pair( node->GetID(), aGhs3dID )); -// aet.push_back(aGhs3dID); -// } -// if (isOK) -// et.push_back(aet); -// } -// } -// -// -// /* Write vertices number */ -// MESSAGE("Number of vertices: "<GetElements(); -// // for ( int j = 0; j < 4; ++j ) -// // { -// // 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(); -// // for ( int j = 0; j < 4; ++j ) { -// // 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); -// // } -// // } -// // -// -// // /* Write quadrilaterals number */ -// // if (qt.size()) { -// // GmfSetKwd(idx, GmfQuadrilaterals, qt.size()); -// // for (int i=0;i & theSmdsToGhs3dIdMap, - const map & theEnforcedNodeIdToGhs3dIdMap, - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges, - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles) -{ - // 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(); - - // count triangles bound to geometry - int nbTriangles = 0; - - TopTools_IndexedMapOfShape facesMap, trianglesMap; - TopExp::MapShapes( theShape, TopAbs_FACE, facesMap ); - - int nbFaces = facesMap.Extent(); - - for ( int i = 1; i <= nbFaces; ++i ) - if (( theSubMesh = theMesh.GetSubMesh( facesMap(i)))) - nbTriangles += theSubMesh->NbElements(); - std::string tmpStr; - (nbFaces == 0 || nbFaces == 1) ? tmpStr = " shape " : tmpStr = " shapes " ; - std::cout << " " << nbFaces << tmpStr << "of 2D dimension"; - int nbEnforcedElements = nbEnforcedEdges+nbEnforcedTriangles; - if (nbEnforcedElements > 0) { - (nbEnforcedElements == 1) ? tmpStr = "shape:" : tmpStr = "shapes:"; - std::cout << " and" << std::endl; - std::cout << " " << nbEnforcedElements - << " enforced " << tmpStr << std::endl; - } - else - std::cout << std::endl; - if (nbEnforcedEdges) { - (nbEnforcedEdges == 1) ? tmpStr = "edge" : tmpStr = "edges"; - std::cout << " " << nbEnforcedEdges << " enforced " << tmpStr << std::endl; - } - if (nbEnforcedTriangles) { - (nbEnforcedTriangles == 1) ? tmpStr = "triangle" : tmpStr = "triangles"; - std::cout << " " << nbEnforcedTriangles << " enforced " << tmpStr << std::endl; - } - std::cout << std::endl; - -// theFile << space << nbTriangles << space << dummyint << std::endl; - std::ostringstream globalStream, localStream, aStream; - - 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(); - - localStream << nbNodes << space; - - itOnSubFace = aFace->nodesIterator(); - for ( int j = 0; j < 3; ++j ) { - // 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 ; -// for ( int j = 0; j < 2; ++j ) { -// 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); -// nbNodes = aFace->NbCornerNodes(); -// isOK = true; -// itOnSubFace = aFace->nodesIterator(); -// aStream.str(""); -// aStream << nbNodes << space ; -// for ( int j = 0; j < 3; ++j ) { -// 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 -// // - - theFile - << nbTriangles/*+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::TGHS3DEnforcedVertexCoordsValues & theEnforcedVertices, - GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap & theEnforcedNodes, - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges, - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles) -{ - // record structure: - // - // NB_NODES - // Loop from 1 to NB_NODES - // X Y Z DUMMY_INT - - SMESHDS_Mesh * theMeshDS = theHelper.GetMeshDS(); - int nbNodes = theMeshDS->NbNodes(); - if ( nbNodes == 0 ) - return false; - - int nbEnforcedVertices = theEnforcedVertices.size(); - int nbEnforcedNodes = theEnforcedNodes.size(); - - const TopoDS_Shape shapeToMesh = theMeshDS->ShapeToMesh(); - - int aGhs3dID = 1; - SMDS_NodeIteratorPtr nodeIt = theMeshDS->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 nodes which are free - // and replace not-free nodes on degenerated edges by the node on vertex - TNodeNodeMap n2nDegen; // map a node on degenerated edge to a node on vertex - TNodeNodeMap::iterator n2nDegenIt; - if ( theHelper.HasDegeneratedEdges() ) - { - 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 ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() ) - { - TopoDS_Shape vertex = TopoDS_Iterator( e.Current() ).Value(); - const SMDS_MeshNode* vNode = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), theMeshDS); - { - SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); - while ( nIt->more() ) - n2nDegen.insert( make_pair( nIt->next(), vNode )); - } - } - } - } - nbNodes -= n2nDegen.size(); - } - - 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() )) - nbNodes -= int( theHelper.IsMedium( node )); - } - nodeIt = theMeshDS->nodesIterator(); - } - - const char* space = " "; - const int dummyint = 0; - - std::string tmpStr; - (nbNodes == 0 || nbNodes == 1) ? tmpStr = " node" : tmpStr = " nodes"; - // NB_NODES - std::cout << std::endl; - std::cout << "The initial 2D mesh contains :" << std::endl; - std::cout << " " << nbNodes << tmpStr << std::endl; - if (nbEnforcedVertices > 0) { - (nbEnforcedVertices == 1) ? tmpStr = "vertex" : tmpStr = "vertices"; - std::cout << " " << nbEnforcedVertices << " enforced " << tmpStr << std::endl; - } - if (nbEnforcedNodes > 0) { - (nbEnforcedNodes == 1) ? tmpStr = "node" : tmpStr = "nodes"; - std::cout << " " << nbEnforcedNodes << " enforced " << tmpStr << std::endl; - } - std::cout << std::endl; - std::cout << "Start writing in 'points' file ..." << std::endl; - - theFile << nbNodes << std::endl; - - // Loop from 1 to NB_NODES - - while ( nodeIt->more() ) - { - node = nodeIt->next(); - if ( isQuadMesh && theHelper.IsMedium( node )) // Issue 0021238 - continue; - if ( n2nDegen.count( node ) ) // 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; - - theFile << std::endl; - - } - - // Iterate over the enforced nodes - std::map enfVertexIndexSizeMap; - if (nbEnforcedNodes) { - GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator nodeIt = theEnforcedNodes.begin(); - for( ; nodeIt != theEnforcedNodes.end() ; ++nodeIt) { - double x = nodeIt->first->X(); - double y = nodeIt->first->Y(); - double z = nodeIt->first->Z(); - // 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 ) - 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->first->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 - << -1 << space - << dummyint << space; - theFile << std::endl; - theEnforcedNodeIdToGhs3dIdMap.insert( make_pair( nodeIt->first->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::TGHS3DEnforcedVertexCoordsValues::const_iterator vertexIt = theEnforcedVertices.begin(); - for( ; 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 ) - 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++; - } - } - - - std::cout << std::endl; - std::cout << "End writing in 'points' file." << std::endl; - - return true; -} - -//======================================================================= -//function : readResultFile -//purpose : readResultFile with geometry -//======================================================================= - -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 tabShape[], - double** tabBox, - const int nbShape, - map & theGhs3dIdToNodeMap, - std::map & theNodeId2NodeIndexMap, - bool toMeshHoles, - int nbEnforcedVertices, - int nbEnforcedNodes, - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedEdges, - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap & theEnforcedTriangles) -{ - MESSAGE("GHS3DPlugin_GHS3D::readResultFile()"); - Kernel_Utils::Localizer loc; - struct stat status; - size_t length; - - std::string tmpStr; - - char *ptr, *mapPtr; - char *tetraPtr; - char *shapePtr; - - SMESHDS_Mesh* theMeshDS = theHelper.GetMeshDS(); - - int nbElems, nbNodes, nbInputNodes; - int nbTriangle; - int ID, shapeID, ghs3dShapeID; - int IdShapeRef = 1; - int compoundID = - nbShape ? theMeshDS->ShapeToIndex( tabShape[0] ) : theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() ); - - int *tab, *tabID, *nodeID, *nodeAssigne; - double *coord; - const SMDS_MeshNode **node; - - tab = new int[3]; - nodeID = new int[4]; - coord = new double[3]; - node = new const SMDS_MeshNode*[4]; - - TopoDS_Shape aSolid; - SMDS_MeshNode * aNewNode; - map ::iterator itOnNode; - SMDS_MeshElement* aTet; -#ifdef _DEBUG_ - set shapeIDs; -#endif - - // Read the file state - 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]; - - nodeAssigne = new int[ nbNodes+1 ]; - - if (nbShape > 0) - aSolid = tabShape[0]; - - // Reading the nodeId - for (int i=0; i < 4*nbElems; i++) - strtol(ptr, &ptr, 10); - - MESSAGE("nbInputNodes: "<computeCanceled()) - return false; -#endif - for (int iCoor=0; iCoor < 3; iCoor++) - coord[ iCoor ] = strtod(ptr, &ptr); - nodeAssigne[ iNode ] = 1; - if ( iNode > (nbInputNodes-(nbEnforcedVertices+nbEnforcedNodes)) ) { - // Creating SMESH nodes - // - for enforced vertices - // - for vertices of forced edges - // - for ghs3d nodes - nodeAssigne[ iNode ] = 0; - aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] ); - theGhs3dIdToNodeMap.insert(theGhs3dIdToNodeMap.end(), make_pair( iNode, aNewNode )); - } - } - - // Reading the number of triangles which corresponds to the number of sub-domains - nbTriangle = strtol(ptr, &ptr, 10); - - tabID = new int[nbTriangle]; - for (int i=0; i < nbTriangle; i++) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if(theAlgo->computeCanceled()) - return false; -#endif - tabID[i] = 0; - // find the solid corresponding to GHS3D sub-domain following - // the technique proposed in GHS3D manual in chapter - // "B.4 Subdomain (sub-region) assignment" - int nodeId1 = strtol(ptr, &ptr, 10); - int nodeId2 = strtol(ptr, &ptr, 10); - int nodeId3 = strtol(ptr, &ptr, 10); - if ( nbTriangle > 1 ) { - const SMDS_MeshNode* n1 = theGhs3dIdToNodeMap[ nodeId1 ]; - const SMDS_MeshNode* n2 = theGhs3dIdToNodeMap[ nodeId2 ]; - const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ]; - if (!n1 || !n2 || !n3) { - tabID[i] = HOLE_ID; - continue; - } - 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 - 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 - } - } - } - - shapePtr = ptr; - - if ( nbTriangle <= nbShape ) // no holes - toMeshHoles = true; // not avoid creating tetras in holes - - // Associating the tetrahedrons to the shapes - shapeID = compoundID; - for (int iElem = 0; iElem < nbElems; iElem++) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if(theAlgo->computeCanceled()) - return false; -#endif - for (int iNode = 0; iNode < 4; iNode++) { - ID = strtol(tetraPtr, &tetraPtr, 10); - itOnNode = theGhs3dIdToNodeMap.find(ID); - node[ iNode ] = itOnNode->second; - nodeID[ iNode ] = ID; - } - // We always run GHS3D with "to mesh holes"==TRUE but we must not create - // tetras within holes depending on hypo option, - // so we first check if aTet is inside a hole and then create it - //aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] ); - if ( nbTriangle > 1 ) { - shapeID = HOLE_ID; // negative shapeID means not to create tetras if !toMeshHoles - ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef; - if ( tabID[ ghs3dShapeID ] == 0 ) { - TopAbs_State state; - aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape, &state); - if ( toMeshHoles || state == TopAbs_IN ) - shapeID = theMeshDS->ShapeToIndex( aSolid ); - tabID[ ghs3dShapeID ] = shapeID; - } - else - shapeID = tabID[ ghs3dShapeID ]; - } - else if ( nbShape > 1 ) { - // Case where nbTriangle == 1 while nbShape == 2 encountered - // with compound of 2 boxes and "To mesh holes"==False, - // so there are no subdomains specified for each tetrahedron. - // Try to guess a solid by a node already bound to shape - shapeID = 0; - for ( int i=0; i<4 && shapeID==0; i++ ) { - if ( nodeAssigne[ nodeID[i] ] == 1 && - node[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE && - node[i]->getshapeId() > 1 ) - { - shapeID = node[i]->getshapeId(); - } - } - if ( shapeID==0 ) { - aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape); - shapeID = theMeshDS->ShapeToIndex( aSolid ); - } - } - // set new nodes and tetrahedron onto the shape - for ( int i=0; i<4; i++ ) { - if ( nodeAssigne[ nodeID[i] ] == 0 ) { - if ( shapeID != HOLE_ID ) - theMeshDS->SetNodeInVolume( node[i], shapeID ); - nodeAssigne[ nodeID[i] ] = shapeID; - } - } - if ( toMeshHoles || shapeID != HOLE_ID ) { - aTet = theHelper.AddVolume( node[1], node[0], node[2], node[3], - /*id=*/0, /*force3d=*/false); - theMeshDS->SetMeshElementOnShape( aTet, shapeID ); - } -#ifdef _DEBUG_ - shapeIDs.insert( shapeID ); -#endif - } - - // Add enforced elements - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap::const_iterator elemIt; - const SMDS_MeshElement* anElem; - SMDS_ElemIteratorPtr itOnEnfElem; - map::const_iterator itOnMap; - shapeID = compoundID; - // Enforced edges - if (theEnforcedEdges.size()) { - (theEnforcedEdges.size() <= 1) ? tmpStr = " enforced edge" : " enforced edges"; - std::cout << "Add " << theEnforcedEdges.size() << tmpStr << std::endl; - std::vector< const SMDS_MeshNode* > node( 2 ); - // Iterate over the enforced edges - for(elemIt = theEnforcedEdges.begin() ; elemIt != theEnforcedEdges.end() ; ++elemIt) { - anElem = elemIt->first; - bool addElem = true; - itOnEnfElem = anElem->nodesIterator(); - for ( int j = 0; j < 2; ++j ) { - int aNodeID = itOnEnfElem->next()->GetID(); - itOnMap = theNodeId2NodeIndexMap.find(aNodeID); - if (itOnMap != theNodeId2NodeIndexMap.end()) { - itOnNode = theGhs3dIdToNodeMap.find((*itOnMap).second); - if (itOnNode != theGhs3dIdToNodeMap.end()) { - node.push_back((*itOnNode).second); -// shapeID =(*itOnNode).second->getshapeId(); - } - else - addElem = false; - } - else - addElem = false; - } - if (addElem) { - aTet = theHelper.AddEdge( node[0], node[1], 0, false); - theMeshDS->SetMeshElementOnShape( aTet, shapeID ); - } - } - } - // Enforced faces - if (theEnforcedTriangles.size()) { - (theEnforcedTriangles.size() <= 1) ? tmpStr = " enforced triangle" : " enforced triangles"; - std::cout << "Add " << theEnforcedTriangles.size() << " enforced triangles" << std::endl; - std::vector< const SMDS_MeshNode* > node( 3 ); - // Iterate over the enforced triangles - for(elemIt = theEnforcedTriangles.begin() ; elemIt != theEnforcedTriangles.end() ; ++elemIt) { - anElem = elemIt->first; - bool addElem = true; - itOnEnfElem = anElem->nodesIterator(); - for ( int j = 0; j < 3; ++j ) { - int aNodeID = itOnEnfElem->next()->GetID(); - itOnMap = theNodeId2NodeIndexMap.find(aNodeID); - if (itOnMap != theNodeId2NodeIndexMap.end()) { - itOnNode = theGhs3dIdToNodeMap.find((*itOnMap).second); - if (itOnNode != theGhs3dIdToNodeMap.end()) { - node.push_back((*itOnNode).second); -// shapeID =(*itOnNode).second->getshapeId(); - } - else - addElem = false; - } - else - addElem = false; - } - if (addElem) { - aTet = theHelper.AddFace( node[0], node[1], node[2], 0, false); - theMeshDS->SetMeshElementOnShape( aTet, shapeID ); - } - } - } - - // Remove nodes of tetras inside holes if !toMeshHoles - if ( !toMeshHoles ) { - itOnNode = theGhs3dIdToNodeMap.find( nbInputNodes ); - for ( ; itOnNode != theGhs3dIdToNodeMap.end(); ++itOnNode) { - ID = itOnNode->first; - if ( nodeAssigne[ ID ] == HOLE_ID ) - theMeshDS->RemoveFreeNode( itOnNode->second, 0 ); - } - } - - - if ( nbElems ) { - (nbElems <= 1) ? tmpStr = " tetrahedra" : " tetrahedrons"; - cout << nbElems << tmpStr << " have been associated to " << nbShape; - (nbShape <= 1) ? tmpStr = " shape" : " shapes"; - cout << tmpStr << endl; - } -#ifdef WNT - UnmapViewOfFile(mapPtr); - CloseHandle(hMapObject); - CloseHandle(fd); -#else - munmap(mapPtr, length); -#endif - close(fileOpen); - - delete [] tab; - delete [] tabID; - delete [] nodeID; - delete [] coord; - delete [] node; - delete [] nodeAssigne; - -#ifdef _DEBUG_ - shapeIDs.erase(-1); - if ( shapeIDs.size() != nbShape ) { - (shapeIDs.size() <= 1) ? tmpStr = " solid" : " solids"; - std::cout << "Only " << shapeIDs.size() << tmpStr << " of " << nbShape << " found" << std::endl; - for (int i=0; iShapeToIndex( tabShape[i] ); - if ( shapeIDs.find( shapeID ) == shapeIDs.end() ) - std::cout << " Solid #" << shapeID << " not found" << std::endl; - } - } -#endif - - return true; -} - - //============================================================================= /*! - *Here we are going to use the GHS3D mesher with geometry + *Here we are going to use the MG-Tetra mesher with geometry */ //============================================================================= @@ -3172,146 +1564,131 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape) { bool Ok(false); - //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS(); - - // we count the number of shapes - // _nbShape = countShape( meshDS, TopAbs_SOLID ); -- 0020330: Pb with ghs3d as a submesh - _nbShape = 0; TopExp_Explorer expBox ( theShape, TopAbs_SOLID ); - for ( ; expBox.More(); expBox.Next() ) - _nbShape++; - - // create bounding box for every shape inside the compound - - int iShape = 0; - TopoDS_Shape* tabShape; - double** tabBox; - tabShape = new TopoDS_Shape[_nbShape]; - tabBox = new double*[_nbShape]; - for (int i=0; i<_nbShape; i++) - tabBox[i] = new double[6]; - Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax; - - for (expBox.ReInit(); expBox.More(); expBox.Next()) { - tabShape[iShape] = expBox.Current(); - Bnd_Box BoundingBox; - BRepBndLib::Add(expBox.Current(), BoundingBox); - BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax); - tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax; - tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax; - tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax; - iShape++; - } // a unique working file name // to avoid access to the same files by eg different users - TCollection_AsciiString aGenericName - = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str(); + _genericName = GHS3DPlugin_Hypothesis::GetFileName(_hyp); + TCollection_AsciiString aGenericName((char*) _genericName.c_str() ); + TCollection_AsciiString aGenericNameRequired = aGenericName + "_required"; - TCollection_AsciiString aResultFileName; TCollection_AsciiString aLogFileName = aGenericName + ".log"; // log - // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) -// TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName; -// TCollection_AsciiString aGenericNameRequired = aGenericName + "_required"; -// #ifdef _DEBUG_ -// aGMFFileName = aGenericName + ".mesh"; // GMF mesh file -// aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file -// aRequiredVerticesFileName = aGenericNameRequired + ".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 = aGenericNameRequired + ".meshb"; // GMF required vertices mesh file -// aSolFileName = aGenericName + "_required.solb"; // GMF solution file -// #endif - - TCollection_AsciiString aFacesFileName, aPointsFileName, 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); - } + TCollection_AsciiString aResultFileName; + TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName, aResSolFileName; + aGMFFileName = aGenericName + ".mesh"; // GMF mesh file + aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file + aResSolFileName = aGenericName + "Vol.sol"; // GMF mesh file + aRequiredVerticesFileName = aGenericNameRequired + ".mesh"; // GMF required vertices mesh file + aSolFileName = aGenericNameRequired + ".sol"; // GMF solution file + std::map aNodeId2NodeIndexMap, aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap; - std::map aGhs3dIdToNodeMap; std::map nodeID2nodeIndexMap; + std::map, std::string> enfVerticesWithGroup; GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues coordsSizeMap = GHS3DPlugin_Hypothesis::GetEnforcedVerticesCoordsSize(_hyp); GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp); GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp); GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp); -// TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp); GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp); + GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList enfVertices = GHS3DPlugin_Hypothesis::GetEnforcedVertices(_hyp); + GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexList::const_iterator enfVerIt = enfVertices.begin(); + std::vector coords; + + for ( ; enfVerIt != enfVertices.end() ; ++enfVerIt) + { + GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertex* enfVertex = (*enfVerIt); + if (enfVertex->coords.size()) { + coordsSizeMap.insert(make_pair(enfVertex->coords,enfVertex->size)); + enfVerticesWithGroup.insert(make_pair(enfVertex->coords,enfVertex->groupName)); + } + else { + TopoDS_Shape GeomShape = entryToShape(enfVertex->geomEntry); + for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){ + coords.clear(); + if (it.Value().ShapeType() == TopAbs_VERTEX){ + gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value())); + coords.push_back(aPnt.X()); + coords.push_back(aPnt.Y()); + coords.push_back(aPnt.Z()); + if (coordsSizeMap.find(coords) == coordsSizeMap.end()) { + coordsSizeMap.insert(make_pair(coords,enfVertex->size)); + enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName)); + } + } + } + } + } int nbEnforcedVertices = coordsSizeMap.size(); int nbEnforcedNodes = enforcedNodes.size(); - + std::string tmpStr; (nbEnforcedNodes <= 1) ? tmpStr = "node" : "nodes"; std::cout << nbEnforcedNodes << " enforced " << tmpStr << " from hypo" << std::endl; (nbEnforcedVertices <= 1) ? tmpStr = "vertex" : "vertices"; std::cout << nbEnforcedVertices << " enforced " << tmpStr << " from hypo" << std::endl; - + SMESH_MesherHelper helper( theMesh ); helper.SetSubShape( theShape ); - { - SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh )); + std::vector aNodeByGhs3dId, anEnforcedNodeByGhs3dId; + std::vector aFaceByGhs3dId; + std::map aNodeToGhs3dIdMap; + std::vector aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId; - // make prisms on quadrangles - if ( theMesh.NbQuadrangles() > 0 ) + MG_Tetra_API mgTetra( _computeCanceled, _progress ); + + _isLibUsed = mgTetra.IsLibrary(); + if ( theMesh.NbQuadrangles() > 0 ) + _progressAdvance /= 10; + if ( _viscousLayersHyp ) + _progressAdvance /= 10; + + // proxyMesh must live till readGMFFile() as a proxy face can be used by + // MG-Tetra for domain indication + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh )); + + // make prisms on quadrangles and viscous layers + if ( theMesh.NbQuadrangles() > 0 || _viscousLayersHyp ) + { + vector components; + for (expBox.ReInit(); expBox.More(); expBox.Next()) { - vector components; - for (expBox.ReInit(); expBox.More(); expBox.Next()) + if ( _viscousLayersHyp ) + { + proxyMesh = _viscousLayersHyp->Compute( theMesh, expBox.Current() ); + if ( !proxyMesh ) + return false; + } + if ( theMesh.NbQuadrangles() > 0 ) { - if ( _viscousLayersHyp ) - { - proxyMesh = _viscousLayersHyp->Compute( theMesh, expBox.Current() ); - if ( !proxyMesh ) - return false; - } StdMeshers_QuadToTriaAdaptor* q2t = new StdMeshers_QuadToTriaAdaptor; - q2t->Compute( theMesh, expBox.Current(), proxyMesh.get() ); + Ok = q2t->Compute( theMesh, expBox.Current(), proxyMesh.get() ); components.push_back( SMESH_ProxyMesh::Ptr( q2t )); + if ( !Ok ) + return false; } - proxyMesh.reset( new SMESH_ProxyMesh( components )); - } - // build viscous layers - else if ( _viscousLayersHyp ) - { - proxyMesh = _viscousLayersHyp->Compute( theMesh, theShape ); - if ( !proxyMesh ) - return false; } - - Ok = (writePoints( aPointsFile, helper, - aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, aGhs3dIdToNodeMap, - nodeIDToSizeMap, - coordsSizeMap, enforcedNodes, enforcedEdges, enforcedTriangles) - && - writeFaces ( aFacesFile, *proxyMesh, theShape, - aSmdsToGhs3dIdMap, anEnforcedNodeIdToGhs3dIdMap, - enforcedEdges, enforcedTriangles )); -// Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(), -// helper, *proxyMesh, -// aNodeId2NodeIndexMap, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap, -// enforcedNodes, enforcedEdges, enforcedTriangles, /*enforcedQuadrangles,*/ -// coordsSizeMap); + proxyMesh.reset( new SMESH_ProxyMesh( components )); } + // build viscous layers + // else if ( _viscousLayersHyp ) + // { + // proxyMesh = _viscousLayersHyp->Compute( theMesh, theShape ); + // if ( !proxyMesh ) + // return false; + // } + + int anInvalidEnforcedFlags = 0; + Ok = writeGMFFile(&mgTetra, + aGMFFileName.ToCString(), + aRequiredVerticesFileName.ToCString(), + aSolFileName.ToCString(), + *proxyMesh, helper, + aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap, + aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId, + enforcedNodes, enforcedEdges, enforcedTriangles, + enfVerticesWithGroup, coordsSizeMap, anInvalidEnforcedFlags); // Write aSmdsToGhs3dIdMap to temp file TCollection_AsciiString aSmdsToGhs3dIdMapFileName; @@ -3323,23 +1700,19 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, return error(SMESH_Comment("Can't write into ") << aSmdsToGhs3dIdMapFileName); } INFOS( "Writing ids relation into " << aSmdsToGhs3dIdMapFileName); - aIdsFile << "Smds Ghs3d" << std::endl; + aIdsFile << "Smds MG-Tetra" << std::endl; map ::const_iterator myit; for (myit=aSmdsToGhs3dIdMap.begin() ; myit != aSmdsToGhs3dIdMap.end() ; ++myit) { aIdsFile << myit->first << " " << myit->second << std::endl; } aIdsFile.close(); - aFacesFile.close(); - aPointsFile.close(); - + if ( ! Ok ) { if ( !_keepFiles ) { -// removeFile( aGMFFileName ); -// removeFile( aRequiredVerticesFileName ); -// removeFile( aSolFileName ); - removeFile( aFacesFileName ); - removeFile( aPointsFileName ); + removeFile( aGMFFileName ); + removeFile( aRequiredVerticesFileName ); + removeFile( aSolFileName ); removeFile( aSmdsToGhs3dIdMapFileName ); } return error(COMPERR_BAD_INPUT_MESH); @@ -3347,75 +1720,57 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, removeFile( aResultFileName ); // needed for boundary recovery module usage // ----------------- - // run ghs3d mesher + // run MG-Tetra mesher // ----------------- - TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp ).c_str() ); - cmd += TCollection_AsciiString(" -f ") + aGenericName; // file to read - cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file - // The output .mesh file does not contain yet the subdomain-info (Ghs3D 4.2) -// cmd += TCollection_AsciiString(" --in ") + aGenericName; -// cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired; -// cmd += TCollection_AsciiString(" --out ") + aResultGMFFileName; -// cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file + TCollection_AsciiString cmd = GHS3DPlugin_Hypothesis::CommandToRun( _hyp, true, mgTetra.IsExecutable() ).c_str(); + if ( mgTetra.IsExecutable() ) + { + cmd += TCollection_AsciiString(" --in ") + aGMFFileName; + if ( nbEnforcedVertices + nbEnforcedNodes) + cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired; + cmd += TCollection_AsciiString(" --out ") + aResultFileName; + } + if ( !_logInStandardOutput ) + { + mgTetra.SetLogFile( aLogFileName.ToCString() ); + cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file + } std::cout << std::endl; - std::cout << "Ghs3d execution..." << std::endl; + std::cout << "MG-Tetra execution..." << std::endl; std::cout << cmd << std::endl; -#ifdef WITH_SMESH_CANCEL_COMPUTE - _compute_canceled = false; -#endif + _computeCanceled = false; - system( cmd.ToCString() ); // run + std::string errStr; + Ok = mgTetra.Compute( cmd.ToCString(), errStr ); // run - std::cout << std::endl; - std::cout << "End of Ghs3d execution !" << std::endl; + if ( _logInStandardOutput && mgTetra.IsLibrary() ) + std::cout << std::endl << mgTetra.GetLog() << std::endl; + if ( Ok ) + std::cout << std::endl << "End of MG-Tetra execution !" << std::endl; // -------------- // read a result // -------------- - // Mapping the result file - - int fileOpen; - 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 { - bool toMeshHoles = - _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles(); + GHS3DPlugin_Hypothesis::TSetStrings groupsToRemove = GHS3DPlugin_Hypothesis::GetGroupsToRemove(_hyp); + bool toMeshHoles = + _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles(); + const bool toMakeGroupsOfDomains = GHS3DPlugin_Hypothesis::GetToMakeGroupsOfDomains( _hyp ); - helper.IsQuadraticSubMesh( theShape ); - helper.SetElementsOnShape( false ); + helper.IsQuadraticSubMesh( theShape ); + helper.SetElementsOnShape( false ); - Ok = readResultFile( fileOpen, -#ifdef WNT - aResultFileName.ToCString(), -#endif -#ifdef WITH_SMESH_CANCEL_COMPUTE - this, -#endif - /*theMesh, */helper, tabShape, tabBox, _nbShape, - aGhs3dIdToNodeMap, aNodeId2NodeIndexMap, - toMeshHoles, - nbEnforcedVertices, nbEnforcedNodes, - enforcedEdges, enforcedTriangles ); - -// Ok = readGMFFile( -// #ifndef GMF_HAS_SUBDOMAIN_INFO -// fileOpen, -// #endif -// aGenericName.ToCString(), theMesh, -// _nbShape, tabShape, tabBox, -// aGhs3dIdToNodeMap, toMeshHoles, -// nbEnforcedVertices, nbEnforcedNodes); - } + Ok = readGMFFile(&mgTetra, + aResultFileName.ToCString(), + this, + &helper, aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap, + aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId, + groupsToRemove, toMakeGroupsOfDomains, toMeshHoles); + removeEmptyGroupsOfDomains( helper.GetMesh(), /*notEmptyAsWell =*/ !toMakeGroupsOfDomains ); @@ -3425,90 +1780,89 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, if ( Ok ) { - if ( !_keepFiles ) + if ( anInvalidEnforcedFlags ) + error( COMPERR_WARNING, flagsToErrorStr( anInvalidEnforcedFlags )); + if ( _removeLogOnSuccess ) removeFile( aLogFileName ); + // if ( _hyp && _hyp->GetToMakeGroupsOfDomains() ) + // error( COMPERR_WARNING, "'toMakeGroupsOfDomains' is ignored since the mesh is on shape" ); } - else if ( OSD_File( aLogFileName ).Size() > 0 ) + else if ( mgTetra.HasLog() ) { - // get problem description from the log file - _Ghs2smdsConvertor conv( aGhs3dIdToNodeMap ); - storeErrorDescription( aLogFileName, conv ); + if( _computeCanceled ) + error( "interruption initiated by user" ); + else + { + // get problem description from the log file + _Ghs2smdsConvertor conv( aNodeByGhs3dId, proxyMesh ); + error( getErrorDescription( _logInStandardOutput ? 0 : aLogFileName.ToCString(), + mgTetra.GetLog(), conv )); + } } - else + else if ( !errStr.empty() ) { // the log file is empty removeFile( aLogFileName ); - INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" ); - error(COMPERR_ALGO_FAILED, "ghs3d: command not found" ); + INFOS( "MG-Tetra Error, " << errStr); + error(COMPERR_ALGO_FAILED, errStr); } if ( !_keepFiles ) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if (! Ok) - if(_compute_canceled) - removeFile( aLogFileName ); -#endif - removeFile( aFacesFileName ); - removeFile( aPointsFileName ); + if (! Ok && _computeCanceled ) + removeFile( aLogFileName ); + removeFile( aGMFFileName ); + removeFile( aRequiredVerticesFileName ); + removeFile( aSolFileName ); + removeFile( aResSolFileName ); removeFile( aResultFileName ); - removeFile( aBadResFileName ); - removeFile( aBbResFileName ); removeFile( aSmdsToGhs3dIdMapFileName ); } - std::cout << "<" << aResultFileName.ToCString() << "> GHS3D output file "; - if ( !Ok ) - std::cout << "not "; - std::cout << "treated !" << std::endl; - std::cout << std::endl; - - _nbShape = 0; // re-initializing _nbShape for the next Compute() method call - delete [] tabShape; - delete [] tabBox; - + if ( mgTetra.IsExecutable() ) + { + std::cout << "<" << aResultFileName.ToCString() << "> MG-Tetra output file "; + if ( !Ok ) + std::cout << "not "; + std::cout << "treated !" << std::endl; + std::cout << std::endl; + } + else + { + std::cout << "MG-Tetra " << ( Ok ? "succeeded" : "failed") << std::endl; + } return Ok; } //============================================================================= /*! - *Here we are going to use the GHS3D mesher w/o geometry + *Here we are going to use the MG-Tetra mesher w/o geometry */ //============================================================================= bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, SMESH_MesherHelper* theHelper) { - MESSAGE("GHS3DPlugin_GHS3D::Compute()"); - - //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS(); - TopoDS_Shape theShape = theHelper->GetSubShape(); + theHelper->IsQuadraticSubMesh( theHelper->GetSubShape() ); // a unique working file name // to avoid access to the same files by eg different users - TCollection_AsciiString aGenericName - = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str(); + _genericName = GHS3DPlugin_Hypothesis::GetFileName(_hyp); + TCollection_AsciiString aGenericName((char*) _genericName.c_str() ); TCollection_AsciiString aGenericNameRequired = aGenericName + "_required"; TCollection_AsciiString aLogFileName = aGenericName + ".log"; // log TCollection_AsciiString aResultFileName; bool Ok; - TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName; -//#ifdef _DEBUG_ + TCollection_AsciiString aGMFFileName, aRequiredVerticesFileName, aSolFileName, aResSolFileName; aGMFFileName = aGenericName + ".mesh"; // GMF mesh file aResultFileName = aGenericName + "Vol.mesh"; // GMF mesh file + aResSolFileName = aGenericName + "Vol.sol"; // GMF mesh file aRequiredVerticesFileName = aGenericNameRequired + ".mesh"; // GMF required vertices mesh file aSolFileName = aGenericNameRequired + ".sol"; // GMF solution file -//#else -// aGMFFileName = aGenericName + ".meshb"; // GMF mesh file -// aResultFileName = aGenericName + "Vol.meshb"; // GMF mesh file -// aRequiredVerticesFileName = aGenericNameRequired + ".meshb"; // GMF required vertices mesh file -// aSolFileName = aGenericNameRequired + ".solb"; // GMF solution file -//#endif std::map nodeID2nodeIndexMap; std::map, std::string> enfVerticesWithGroup; GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertexCoordsValues coordsSizeMap; TopoDS_Shape GeomShape; -// TopAbs_ShapeEnum GeomType; std::vector coords; gp_Pnt aPnt; GHS3DPlugin_Hypothesis::TGHS3DEnforcedVertex* enfVertex; @@ -3519,70 +1873,32 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, for ( ; enfVerIt != enfVertices.end() ; ++enfVerIt) { enfVertex = (*enfVerIt); -// if (enfVertex->geomEntry.empty() && enfVertex->coords.size()) { if (enfVertex->coords.size()) { coordsSizeMap.insert(make_pair(enfVertex->coords,enfVertex->size)); enfVerticesWithGroup.insert(make_pair(enfVertex->coords,enfVertex->groupName)); -// MESSAGE("enfVerticesWithGroup.insert(make_pair(("<coords[0]<<","<coords[1]<<","<coords[2]<<"),\""<groupName<<"\"))"); } else { -// if (!enfVertex->geomEntry.empty()) { GeomShape = entryToShape(enfVertex->geomEntry); -// GeomType = GeomShape.ShapeType(); - -// if (!enfVertex->isCompound) { -// // if (GeomType == TopAbs_VERTEX) { -// coords.clear(); -// aPnt = BRep_Tool::Pnt(TopoDS::Vertex(GeomShape)); -// coords.push_back(aPnt.X()); -// coords.push_back(aPnt.Y()); -// coords.push_back(aPnt.Z()); -// if (coordsSizeMap.find(coords) == coordsSizeMap.end()) { -// coordsSizeMap.insert(make_pair(coords,enfVertex->size)); -// enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName)); -// } -// } -// -// // Group Management -// else { -// if (GeomType == TopAbs_COMPOUND){ - for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){ - coords.clear(); - if (it.Value().ShapeType() == TopAbs_VERTEX){ - aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value())); - coords.push_back(aPnt.X()); - coords.push_back(aPnt.Y()); - coords.push_back(aPnt.Z()); - if (coordsSizeMap.find(coords) == coordsSizeMap.end()) { - coordsSizeMap.insert(make_pair(coords,enfVertex->size)); - enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName)); -// MESSAGE("enfVerticesWithGroup.insert(make_pair(("<groupName<<"\"))"); - } + for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()){ + coords.clear(); + if (it.Value().ShapeType() == TopAbs_VERTEX){ + aPnt = BRep_Tool::Pnt(TopoDS::Vertex(it.Value())); + coords.push_back(aPnt.X()); + coords.push_back(aPnt.Y()); + coords.push_back(aPnt.Z()); + if (coordsSizeMap.find(coords) == coordsSizeMap.end()) { + coordsSizeMap.insert(make_pair(coords,enfVertex->size)); + enfVerticesWithGroup.insert(make_pair(coords,enfVertex->groupName)); } } -// } + } } } -// const SMDS_MeshNode* enfNode; - GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp); -// GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap::const_iterator enfNodeIt = enforcedNodes.begin(); -// for ( ; enfNodeIt != enforcedNodes.end() ; ++enfNodeIt) -// { -// enfNode = enfNodeIt->first; -// coords.clear(); -// coords.push_back(enfNode->X()); -// coords.push_back(enfNode->Y()); -// coords.push_back(enfNode->Z()); -// if (enfVerticesWithGro -// enfVerticesWithGroup.insert(make_pair(coords,enfNodeIt->second)); -// } - - - GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp); + GHS3DPlugin_Hypothesis::TIDSortedNodeGroupMap enforcedNodes = GHS3DPlugin_Hypothesis::GetEnforcedNodes(_hyp); + GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedEdges = GHS3DPlugin_Hypothesis::GetEnforcedEdges(_hyp); GHS3DPlugin_Hypothesis::TIDSortedElemGroupMap enforcedTriangles = GHS3DPlugin_Hypothesis::GetEnforcedTriangles(_hyp); -// TIDSortedElemSet enforcedQuadrangles = GHS3DPlugin_Hypothesis::GetEnforcedQuadrangles(_hyp); - GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp); + GHS3DPlugin_Hypothesis::TID2SizeMap nodeIDToSizeMap = GHS3DPlugin_Hypothesis::GetNodeIDToSizeMap(_hyp); std::string tmpStr; @@ -3594,64 +1910,85 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, std::cout << nbEnforcedVertices << " enforced " << tmpStr << " from hypo" << std::endl; std::vector aNodeByGhs3dId, anEnforcedNodeByGhs3dId; + std::vector aFaceByGhs3dId; std::map aNodeToGhs3dIdMap; std::vector aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId; - { - SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh )); - if ( theMesh.NbQuadrangles() > 0 ) - { - StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor; - aQuad2Trias->Compute( theMesh ); - proxyMesh.reset( aQuad2Trias ); - } - Ok = writeGMFFile(aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(), - *proxyMesh, &theMesh, - aNodeByGhs3dId, aNodeToGhs3dIdMap, - aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId, - enforcedNodes, enforcedEdges, enforcedTriangles, - enfVerticesWithGroup, coordsSizeMap); + + MG_Tetra_API mgTetra( _computeCanceled, _progress ); + + _isLibUsed = mgTetra.IsLibrary(); + if ( theMesh.NbQuadrangles() > 0 ) + _progressAdvance /= 10; + + // proxyMesh must live till readGMFFile() as a proxy face can be used by + // MG-Tetra for domain indication + SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( theMesh )); + if ( theMesh.NbQuadrangles() > 0 ) + { + StdMeshers_QuadToTriaAdaptor* aQuad2Trias = new StdMeshers_QuadToTriaAdaptor; + Ok = aQuad2Trias->Compute( theMesh ); + proxyMesh.reset( aQuad2Trias ); + if ( !Ok ) + return false; } + int anInvalidEnforcedFlags = 0; + Ok = writeGMFFile(&mgTetra, + aGMFFileName.ToCString(), aRequiredVerticesFileName.ToCString(), aSolFileName.ToCString(), + *proxyMesh, *theHelper, + aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap, + aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId, + enforcedNodes, enforcedEdges, enforcedTriangles, + enfVerticesWithGroup, coordsSizeMap, anInvalidEnforcedFlags); + // ----------------- - // run ghs3d mesher + // run MG-Tetra mesher // ----------------- - TCollection_AsciiString cmd = TCollection_AsciiString((char*)GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false ).c_str()); - - cmd += TCollection_AsciiString(" --in ") + aGenericName; - if ( nbEnforcedVertices + nbEnforcedNodes) - cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired; - cmd += TCollection_AsciiString(" --out ") + aResultFileName; - cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file + TCollection_AsciiString cmd = GHS3DPlugin_Hypothesis::CommandToRun( _hyp, false, mgTetra.IsExecutable() ).c_str(); + if ( mgTetra.IsExecutable() ) + { + cmd += TCollection_AsciiString(" --in ") + aGMFFileName; + if ( nbEnforcedVertices + nbEnforcedNodes) + cmd += TCollection_AsciiString(" --required_vertices ") + aGenericNameRequired; + cmd += TCollection_AsciiString(" --out ") + aResultFileName; + } + if ( !_logInStandardOutput ) + { + mgTetra.SetLogFile( aLogFileName.ToCString() ); + cmd += TCollection_AsciiString(" 1>" ) + aLogFileName; // dump into file + } std::cout << std::endl; - std::cout << "Ghs3d execution..." << std::endl; + std::cout << "MG-Tetra execution..." << std::endl; std::cout << cmd << std::endl; -#ifdef WITH_SMESH_CANCEL_COMPUTE - _compute_canceled = false; -#endif + _computeCanceled = false; - system( cmd.ToCString() ); // run + std::string errStr; + Ok = mgTetra.Compute( cmd.ToCString(), errStr ); // run - std::cout << std::endl; - std::cout << "End of Ghs3d execution !" << std::endl; + if ( _logInStandardOutput && mgTetra.IsLibrary() ) + std::cout << std::endl << mgTetra.GetLog() << std::endl; + if ( Ok ) + std::cout << std::endl << "End of MG-Tetra execution !" << std::endl; // -------------- // read a result // -------------- GHS3DPlugin_Hypothesis::TSetStrings groupsToRemove = GHS3DPlugin_Hypothesis::GetGroupsToRemove(_hyp); + const bool toMakeGroupsOfDomains = GHS3DPlugin_Hypothesis::GetToMakeGroupsOfDomains( _hyp ); - Ok = readGMFFile(aResultFileName.ToCString(), -#ifdef WITH_SMESH_CANCEL_COMPUTE - this, -#endif - theHelper, theShape, aNodeByGhs3dId, aNodeToGhs3dIdMap, - aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId, - groupsToRemove); + Ok = Ok && readGMFFile(&mgTetra, + aResultFileName.ToCString(), + this, + theHelper, aNodeByGhs3dId, aFaceByGhs3dId, aNodeToGhs3dIdMap, + aNodeGroupByGhs3dId, anEdgeGroupByGhs3dId, aFaceGroupByGhs3dId, + groupsToRemove, toMakeGroupsOfDomains); updateMeshGroups(theHelper->GetMesh(), groupsToRemove); + removeEmptyGroupsOfDomains( theHelper->GetMesh(), /*notEmptyAsWell =*/ !toMakeGroupsOfDomains ); if ( Ok ) { GHS3DPlugin_Hypothesis* that = (GHS3DPlugin_Hypothesis*)this->_hyp; @@ -3664,60 +2001,64 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh& theMesh, if ( Ok ) { - if ( !_keepFiles ) + if ( anInvalidEnforcedFlags ) + error( COMPERR_WARNING, flagsToErrorStr( anInvalidEnforcedFlags )); + if ( _removeLogOnSuccess ) removeFile( aLogFileName ); + + //if ( !toMakeGroupsOfDomains && _hyp && _hyp->GetToMakeGroupsOfDomains() ) + //error( COMPERR_WARNING, "'toMakeGroupsOfDomains' is ignored since 'toMeshHoles' is OFF." ); } - else if ( OSD_File( aLogFileName ).Size() > 0 ) + else if ( mgTetra.HasLog() ) { - // get problem description from the log file - _Ghs2smdsConvertor conv( aNodeByGhs3dId ); - storeErrorDescription( aLogFileName, conv ); + if( _computeCanceled ) + error( "interruption initiated by user" ); + else + { + // get problem description from the log file + _Ghs2smdsConvertor conv( aNodeByGhs3dId, proxyMesh ); + error( getErrorDescription( _logInStandardOutput ? 0 : aLogFileName.ToCString(), + mgTetra.GetLog(), conv )); + } } else { // the log file is empty removeFile( aLogFileName ); - INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" ); - error(COMPERR_ALGO_FAILED, "ghs3d: command not found" ); + INFOS( "MG-Tetra Error, " << errStr); + error(COMPERR_ALGO_FAILED, errStr); } if ( !_keepFiles ) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if (! Ok) - if(_compute_canceled) - removeFile( aLogFileName ); -#endif + if (! Ok && _computeCanceled) + removeFile( aLogFileName ); removeFile( aGMFFileName ); removeFile( aResultFileName ); removeFile( aRequiredVerticesFileName ); removeFile( aSolFileName ); + removeFile( aResSolFileName ); } return Ok; } -#ifdef WITH_SMESH_CANCEL_COMPUTE void GHS3DPlugin_GHS3D::CancelCompute() { - _compute_canceled = true; -#ifdef WNT -#else - TCollection_AsciiString aGenericName - = (char*) GHS3DPlugin_Hypothesis::GetFileName(_hyp).c_str(); - TCollection_AsciiString cmd = - TCollection_AsciiString("ps ux | grep ") + aGenericName; - cmd += TCollection_AsciiString(" | grep -v grep | awk '{print $2}' | xargs kill -9 > /dev/null 2>&1"); - system( cmd.ToCString() ); + _computeCanceled = true; +#if !defined WIN32 && !defined __APPLE__ + std::string cmd = "ps xo pid,args | grep " + _genericName; + //cmd += " | grep -e \"^ *[0-9]\\+ \\+" + GHS3DPlugin_Hypothesis::GetExeName() + "\""; + cmd += " | awk '{print $1}' | xargs kill -9 > /dev/null 2>&1"; + system( cmd.c_str() ); #endif } -#endif //================================================================================ /*! - * \brief Provide human readable text by error code reported by ghs3d + * \brief Provide human readable text by error code reported by MG-Tetra */ //================================================================================ -static string translateError(const int errNum) +static const char* translateError(const int errNum) { switch ( errNum ) { case 0: @@ -3862,6 +2203,146 @@ static string translateError(const int errNum) "The surface mesh is probably very bad in terms of quality."; case 23602: return "Bad vertex number."; + case 1001200: + return "Cannot close mesh file NomFil."; + case 1002010: + return "There are wrong data."; + case 1002120: + return "The number of faces is negative or null."; + case 1002170: + return "The number of vertices is negative or null in the '.sol' file."; + case 1002190: + return "The number of tetrahedra is negative or null."; + case 1002210: + return "The number of vertices is negative or null."; + case 1002211: + return "A face has a vertex negative or null."; + case 1002270: + return "The field is not a size in file NomFil."; + case 1002280: + return "A count is wrong in the enclosing box in the .boite.mesh input " + "file (option '--read_boite')."; + case 1002290: + return "A tetrahedron has a vertex with a negative number."; + case 1002300: + return "the 'MeshVersionFormatted' is not 1 or 2 in the '.mesh' file or the '.sol'."; + case 1002370: + return "The number of values in the '.sol' (metric file) is incompatible with " + "the expected value of number of mesh vertices in the '.mesh' file."; + case 1003000: + return "Not enough memory."; + case 1003020: + return "Not enough memory for the face table."; + case 1003050: + return "Insufficient memory ressources detected due to a bad quality " + "surface mesh leading to too many swaps."; + case 1005010: + return "The surface coordinates of a vertex are differing from the " + "volume coordinates, probably due to a precision problem."; + case 1005050: + return "Invalid dimension. Dimension 3 expected."; + case 1005100: + return "A point has a tag 0. This point is probably outside the domain which has been meshed."; + case 1005103: + return "The vertices of an element are too close to one another or coincident."; + case 1005104: + return "There are at least two points whose distance is very small, and considered as coincident."; + case 1005105: + return "Two vertices are too close to one another or coincident."; + case 1005106: + return "A vertex cannot be inserted."; + case 1005107: + return "Two vertices are too close to one another or coincident. Note : When " + "this error occurs during the overconstrained processing phase, this is only " + "a warning which means that it is difficult to break some overconstrained facets."; + case 1005110: + return "Two surface edges are intersecting."; + case 1005120: + return "A surface edge intersects a surface face."; + case 1005150: + return "A boundary point lies within a surface face."; + case 1005160: + return "A boundary point lies within a surface edge."; + case 1005200: + return "A surface mesh appears more than once in the input surface mesh."; + case 1005210: + return "An edge appears more than once in the input surface mesh."; + case 1005225: + return "Surface with unvalid triangles."; + case 1005270: + return "The metric in the '.sol' file contains more than one field."; + case 1005300: + return "The surface mesh includes at least one hole. The domain is not well defined."; + case 1005301: + return "Presumably, the surface mesh is not compatible with the domain being processed (warning)."; + case 1005302: + return "Probable faces overlapping somewher."; + case 1005320: + return "The quadratic version does not work with prescribed free edges."; + case 1005321: + return "The quadratic version does not work with a volume mesh."; + case 1005370: + return "The metric in the '.sol' file is inadequate (values not per vertices)."; + case 1005371: + return "The number of vertices in the '.sol' is different from the one in the " + "'.mesh' file for the required vertices (option '--required_vertices')."; + case 1005372: + return "More than one type in file NomFil. The type must be equal to 1 in the '.sol'" + "for the required vertices (option '--required_vertices')."; + case 1005515: + return "Bad vertex number."; + case 1005560: + return "No guess to start the definition of the connected component(s)."; + case 1005602: + return "Some initial points cannot be inserted."; + case 1005620: + return "A too bad quality face is detected. This face is considered degenerated."; + case 1005621: + return "A too bad quality face is detected. This face is degenerated."; + case 1005622: + return "The algorithm cannot run further."; + case 1005690: + return "A too small volume element is detected."; + case 1005691: + return "A tetrahedra is suspected to be very bad shaped or wrong."; + case 1005692: + return "There is at least a null or negative volume element. The resulting mesh" + "may be inappropriate."; + case 1005693: + return "There are some null or negative volume element. The resulting mesh may" + "be inappropriate."; + case 1005820: + return "An edge is unique (i.e., bounds a hole in the surface)."; + case 1007000: + return "Abnormal or internal error."; + case 1007010: + return "Too many components with respect to too many sub-domain."; + case 1007400: + return "An internal error has been encountered or a signal has been received. " + "Current mesh will not be saved."; + case 1008491: + return "Impossible to define a component."; + case 1008410: + return "There are some overconstrained edges."; + case 1008420: + return "There are some overconstrained facets."; + case 1008422: + return "Give the number of missing faces (information given when regeneration phase failed)."; + case 1008423: + return "A constrained face cannot be enforced (information given when regeneration phase failed)."; + case 1008441: + return "A constrained edge cannot be enforced."; + case 1008460: + return "It is dramatically tedious to enforce the boundary items."; + case 1008480: + return "The surface mesh regeneration step has failed. A .boite.mesh and .boite.map files are created."; + case 1008490: + return "Invalid resulting mesh."; + case 1008495: + return "P2 correction not successful."; + case 1009000: + return "Program has received an interruption or a termination signal sent by the " + "user or the system administrator. Current mesh will not be saved."; } return ""; } @@ -3893,39 +2374,42 @@ static char* getIds( char* ptr, int nbIds, vector& ids ) */ //================================================================================ -bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& logFile, - const _Ghs2smdsConvertor & toSmdsConvertor ) +SMESH_ComputeErrorPtr +GHS3DPlugin_GHS3D::getErrorDescription(const char* logFile, + const std::string& log, + const _Ghs2smdsConvertor & toSmdsConvertor, + const bool isOk/* = false*/ ) { -#ifdef WITH_SMESH_CANCEL_COMPUTE - if(_compute_canceled) - return error(SMESH_Comment("interruption initiated by user")); -#endif - // open file -#ifdef WNT - int file = ::_open (logFile.ToCString(), _O_RDONLY|_O_BINARY); -#else - int file = ::open (logFile.ToCString(), O_RDONLY); -#endif - if ( file < 0 ) - return error( SMESH_Comment("See ") << logFile << " for problem description"); - - // get file size -// struct stat status; -// fstat(file, &status); -// size_t length = status.st_size; - off_t length = lseek( file, 0, SEEK_END); - lseek( file, 0, SEEK_SET); - - // read file - vector< char > buf( length ); - int nBytesRead = ::read (file, & buf[0], length); - ::close (file); - char* ptr = & buf[0]; - char* bufEnd = ptr + nBytesRead; + SMESH_BadInputElements* badElemsErr = + new SMESH_BadInputElements( toSmdsConvertor.getMesh(), COMPERR_ALGO_FAILED ); + SMESH_ComputeErrorPtr err( badElemsErr ); + + char* ptr = const_cast( log.c_str() ); + char* buf = ptr, * bufEnd = ptr + log.size(); + SMESH_Comment errDescription; - enum { NODE = 1, EDGE, TRIA, VOL, ID = 1 }; + enum { NODE = 1, EDGE, TRIA, VOL, SKIP_ID = 1 }; + + // look for MeshGems version + // Since "MG-TETRA -- MeshGems 1.1-3 (January, 2013)" error codes change. + // To discriminate old codes from new ones we add 1000000 to the new codes. + // This way value of the new codes is same as absolute value of codes printed + // in the log after "MGMESSAGE" string. + int versionAddition = 0; + { + char* verPtr = ptr; + while ( ++verPtr < bufEnd ) + { + if ( strncmp( verPtr, "MG-TETRA -- MeshGems ", 21 ) != 0 ) + continue; + if ( strcmp( verPtr, "MG-TETRA -- MeshGems 1.1-3 " ) >= 0 ) + versionAddition = 1000000; + ptr = verPtr; + break; + } + } // look for errors "ERR #" @@ -3936,115 +2420,113 @@ bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& log if ( strncmp( ptr, "ERR ", 4 ) != 0 ) continue; - list badElems; + list& badElems = badElemsErr->myBadElements; vector nodeIds; ptr += 4; char* errBeg = ptr; - int errNum = strtol(ptr, &ptr, 10); - switch ( errNum ) { // we treat errors enumerated in [SALOME platform 0019316] issue - case 0015: - // The face number (numfac) with vertices (f 1, f 2, f 3) has a null vertex. - ptr = getIds(ptr, NODE, nodeIds); + int errNum = strtol(ptr, &ptr, 10) + versionAddition; + // we treat errors enumerated in [SALOME platform 0019316] issue + // and all errors from a new (Release 1.1) MeshGems User Manual + switch ( errNum ) { + case 0015: // The face number (numfac) with vertices (f 1, f 2, f 3) has a null vertex. + case 1005620 : // a too bad quality face is detected. This face is considered degenerated. + ptr = getIds(ptr, SKIP_ID, nodeIds); ptr = getIds(ptr, TRIA, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 1000: // ERR 1000 : 1 3 2 - // Face (f 1, f 2, f 3) appears more than once in the input surface mesh. + case 1005621 : // a too bad quality face is detected. This face is degenerated. + // hence the is degenerated it is invisible, add its edges in addition + ptr = getIds(ptr, SKIP_ID, nodeIds); ptr = getIds(ptr, TRIA, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); + { + vector edgeNodes( nodeIds.begin(), --nodeIds.end() ); // 01 + badElems.push_back( toSmdsConvertor.getElement(edgeNodes)); + edgeNodes[1] = nodeIds[2]; // 02 + badElems.push_back( toSmdsConvertor.getElement(edgeNodes)); + edgeNodes[0] = nodeIds[1]; // 12 + } break; - case 1001: - // Edge (e1, e2) appears more than once in the input surface mesh - ptr = getIds(ptr, EDGE, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - break; - case 1002: - // Face (f 1, f 2, f 3) has a vertex negative or null + case 1000: // Face (f 1, f 2, f 3) appears more than once in the input surface mesh. + // ERR 1000 : 1 3 2 + case 1002: // Face (f 1, f 2, f 3) has a vertex negative or null + case 3019: // Constrained face (f 1, f 2, f 3) cannot be enforced + case 1002211: // a face has a vertex negative or null. + case 1005200 : // a surface mesh appears more than once in the input surface mesh. + case 1008423 : // a constrained face cannot be enforced (regeneration phase failed). ptr = getIds(ptr, TRIA, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 2004: - // Vertex v1 and vertex v2 are too close to one another or coincident (warning). - ptr = getIds(ptr, NODE, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - ptr = getIds(ptr, NODE, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - break; - case 2012: - // Vertex v1 cannot be inserted (warning). - ptr = getIds(ptr, NODE, nodeIds); + case 1001: // Edge (e1, e2) appears more than once in the input surface mesh + case 3009: // Constrained edge (e1, e2) cannot be enforced (warning). + // ERR 3109 : EDGE 5 6 UNIQUE + case 3109: // Edge (e1, e2) is unique (i.e., bounds a hole in the surface) + case 1005210 : // an edge appears more than once in the input surface mesh. + case 1005820 : // an edge is unique (i.e., bounds a hole in the surface). + case 1008441 : // a constrained edge cannot be enforced. + ptr = getIds(ptr, EDGE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 2014: - // There are at least two points whose distance is dist, i.e., considered as coincident - case 2103: // ERR 2103 : 16 WITH 3 - // Vertex v1 and vertex v2 are too close to one another or coincident (warning). + case 2004: // Vertex v1 and vertex v2 are too close to one another or coincident (warning). + case 2014: // at least two points whose distance is dist, i.e., considered as coincident + case 2103: // Vertex v1 and vertex v2 are too close to one another or coincident (warning). + // ERR 2103 : 16 WITH 3 + case 1005105 : // two vertices are too close to one another or coincident. + case 1005107: // Two vertices are too close to one another or coincident. ptr = getIds(ptr, NODE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); ptr = getIds(ptr, NODE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 3009: - // Constrained edge (e1, e2) cannot be enforced (warning). - ptr = getIds(ptr, EDGE, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - break; - case 3019: - // Constrained face (f 1, f 2, f 3) cannot be enforced - ptr = getIds(ptr, TRIA, nodeIds); + case 2012: // Vertex v1 cannot be inserted (warning). + case 1005106 : // a vertex cannot be inserted. + ptr = getIds(ptr, NODE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 3103: // ERR 3103 : 1 2 WITH 7 3 - // The surface edge (e1, e2) intersects another surface edge (e3, e4) + case 3103: // The surface edge (e1, e2) intersects another surface edge (e3, e4) + case 1005110 : // two surface edges are intersecting. + // ERR 3103 : 1 2 WITH 7 3 ptr = getIds(ptr, EDGE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); ptr = getIds(ptr, EDGE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 3104: // ERR 3104 : 9 10 WITH 1 2 3 - // The surface edge (e1, e2) intersects the surface face (f 1, f 2, f 3) + case 3104: // The surface edge (e1, e2) intersects the surface face (f 1, f 2, f 3) + // ERR 3104 : 9 10 WITH 1 2 3 + case 3106: // One surface edge (say e1, e2) intersects a surface face (f 1, f 2, f 3) + case 1005120 : // a surface edge intersects a surface face. ptr = getIds(ptr, EDGE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); ptr = getIds(ptr, TRIA, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 3105: // ERR 3105 : 8 IN 2 3 5 - // One boundary point (say p1) lies within a surface face (f 1, f 2, f 3) + case 3105: // One boundary point (say p1) lies within a surface face (f 1, f 2, f 3) + // ERR 3105 : 8 IN 2 3 5 + case 1005150 : // a boundary point lies within a surface face. ptr = getIds(ptr, NODE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); ptr = getIds(ptr, TRIA, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 3106: - // One surface edge (say e1, e2) intersects a surface face (f 1, f 2, f 3) - ptr = getIds(ptr, EDGE, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - ptr = getIds(ptr, TRIA, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - break; - case 3107: // ERR 3107 : 2 IN 4 1 - // One boundary point (say p1) lies within a surface edge (e1, e2) (stop). + case 3107: // One boundary point (say p1) lies within a surface edge (e1, e2) (stop). + // ERR 3107 : 2 IN 4 1 + case 1005160 : // a boundary point lies within a surface edge. ptr = getIds(ptr, NODE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); ptr = getIds(ptr, EDGE, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); break; - case 3109: // ERR 3109 : EDGE 5 6 UNIQUE - // Edge (e1, e2) is unique (i.e., bounds a hole in the surface) - ptr = getIds(ptr, EDGE, nodeIds); - badElems.push_back( toSmdsConvertor.getElement(nodeIds)); - break; case 9000: // ERR 9000 // ELEMENT 261 WITH VERTICES : 7 396 -8 242 // VOLUME : -1.11325045E+11 W.R.T. EPSILON 0. // A too small volume element is detected. Are reported the index of the element, // its four vertex indices, its volume and the tolerance threshold value - ptr = getIds(ptr, ID, nodeIds); + ptr = getIds(ptr, SKIP_ID, nodeIds); ptr = getIds(ptr, VOL, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); // even if all nodes found, volume it most probably invisible, - // add its faces to demenstrate it anyhow + // add its faces to demonstrate it anyhow { vector faceNodes( nodeIds.begin(), --nodeIds.end() ); // 012 badElems.push_back( toSmdsConvertor.getElement(faceNodes)); @@ -4077,7 +2559,7 @@ bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& log // SMALL INRADIUS : 0. // A too bad quality face is detected. This face is degenerated, // its index, its three vertex indices together with its inradius are reported - ptr = getIds(ptr, ID, nodeIds); + ptr = getIds(ptr, SKIP_ID, nodeIds); ptr = getIds(ptr, TRIA, nodeIds); badElems.push_back( toSmdsConvertor.getElement(nodeIds)); // add triangle edges as it most probably has zero area and hence invisible @@ -4091,6 +2573,12 @@ bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& log badElems.push_back( toSmdsConvertor.getElement(edgeNodes)); } break; + case 1005103 : // the vertices of an element are too close to one another or coincident. + ptr = getIds(ptr, TRIA, nodeIds); + if ( nodeIds.back() == 0 ) // index of the third vertex of the element (0 for an edge) + nodeIds.resize( EDGE ); + badElems.push_back( toSmdsConvertor.getElement(nodeIds)); + break; } bool isNewError = foundErrorStr.insert( string( errBeg, ptr )).second; @@ -4106,13 +2594,6 @@ bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& log // continue; // not to report different types of errors with bad elements // } - // store bad elements - //if ( allElemsOk ) { - list::iterator elem = badElems.begin(); - for ( ; elem != badElems.end(); ++elem ) - addBadInputElement( *elem ); - //} - // make error text string text = translateError( errNum ); if ( errDescription.find( text ) == text.npos ) { @@ -4133,16 +2614,24 @@ bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& log { char msg2[] = "SEGMENTATION FAULT"; if ( search( &buf[0], bufEnd, msg2, msg2 + strlen(msg2)) != bufEnd ) - errDescription << "ghs3d: SEGMENTATION FAULT. "; + errDescription << "MG-Tetra: SEGMENTATION FAULT. "; } } - if ( errDescription.empty() ) - errDescription << "See " << logFile << " for problem description"; - else - errDescription << "\nSee " << logFile << " for more information"; + if ( !isOk && logFile && logFile[0] ) + { + if ( errDescription.empty() ) + errDescription << "See " << logFile << " for problem description"; + else + errDescription << "\nSee " << logFile << " for more information"; + } - return error( errDescription ); + err->myComment = errDescription; + + if ( err->myComment.empty() && !err->HasBadElems() ) + err = SMESH_ComputeError::New(); // OK + + return err; } //================================================================================ @@ -4151,8 +2640,9 @@ bool GHS3DPlugin_GHS3D::storeErrorDescription(const TCollection_AsciiString& log */ //================================================================================ -_Ghs2smdsConvertor::_Ghs2smdsConvertor( const map & ghs2NodeMap) - :_ghs2NodeMap( & ghs2NodeMap ), _nodeByGhsId( 0 ) +_Ghs2smdsConvertor::_Ghs2smdsConvertor( const map & ghs2NodeMap, + SMESH_ProxyMesh::Ptr mesh) + :_ghs2NodeMap( & ghs2NodeMap ), _nodeByGhsId( 0 ), _mesh( mesh ) { } @@ -4162,14 +2652,15 @@ _Ghs2smdsConvertor::_Ghs2smdsConvertor( const map & g */ //================================================================================ -_Ghs2smdsConvertor::_Ghs2smdsConvertor( const vector & nodeByGhsId) - : _ghs2NodeMap( 0 ), _nodeByGhsId( &nodeByGhsId ) +_Ghs2smdsConvertor::_Ghs2smdsConvertor( const vector & nodeByGhsId, + SMESH_ProxyMesh::Ptr mesh) + : _ghs2NodeMap( 0 ), _nodeByGhsId( &nodeByGhsId ), _mesh( mesh ) { } //================================================================================ /*! - * \brief Return SMDS element by ids of GHS3D nodes + * \brief Return SMDS element by ids of MG-Tetra nodes */ //================================================================================ @@ -4186,7 +2677,7 @@ const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector& ghsNod nodes[ i ] = in->second; } else { - if ( ghsNode < 1 || ghsNode > _nodeByGhsId->size() ) + if ( ghsNode < 1 || ghsNode > (int)_nodeByGhsId->size() ) return 0; nodes[ i ] = (*_nodeByGhsId)[ ghsNode-1 ]; } @@ -4196,13 +2687,13 @@ const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector& ghsNod if ( nbNodes == 2 ) { const SMDS_MeshElement* edge= SMDS_Mesh::FindEdge( nodes[0], nodes[1] ); - if ( !edge ) + if ( !edge || edge->GetID() < 1 || _mesh->IsTemporary( edge )) edge = new SMDS_LinearEdge( nodes[0], nodes[1] ); return edge; } if ( nbNodes == 3 ) { const SMDS_MeshElement* face = SMDS_Mesh::FindFace( nodes ); - if ( !face ) + if ( !face || face->GetID() < 1 || _mesh->IsTemporary( face )) face = new SMDS_FaceOfNodes( nodes[0], nodes[1], nodes[2] ); return face; } @@ -4212,6 +2703,16 @@ const SMDS_MeshElement* _Ghs2smdsConvertor::getElement(const vector& ghsNod return 0; } +//================================================================================ +/*! + * \brief Return a mesh + */ +//================================================================================ + +const SMDS_Mesh* _Ghs2smdsConvertor::getMesh() const +{ + return _mesh->GetMeshDS(); +} //============================================================================= /*! @@ -4295,21 +2796,11 @@ bool GHS3DPlugin_GHS3D::Evaluate(SMESH_Mesh& aMesh, bool GHS3DPlugin_GHS3D::importGMFMesh(const char* theGMFFileName, SMESH_Mesh& theMesh) { - SMESH_MesherHelper* helper = new SMESH_MesherHelper(theMesh ); -// TopoDS_Shape theShape = theMesh.GetShapeToMesh(); - std::vector dummyNodeVector; - std::map dummyNodeMap; - std::map, std::string> dummyEnfVertGroup; - std::vector dummyElemGroup; - std::set dummyGroupsToRemove; - - bool ok = readGMFFile(theGMFFileName, -#ifdef WITH_SMESH_CANCEL_COMPUTE - this, -#endif - helper, theMesh.GetShapeToMesh(), dummyNodeVector, dummyNodeMap, dummyElemGroup, dummyElemGroup, dummyElemGroup, dummyGroupsToRemove); + SMESH_ComputeErrorPtr err = theMesh.GMFToMesh( theGMFFileName, /*makeRequiredGroups =*/ true ); + theMesh.GetMeshDS()->Modified(); - return ok; + + return ( !err || err->IsOK()); } namespace @@ -4357,13 +2848,41 @@ namespace */ static GHS3DPlugin_Hypothesis* GetGHSHypothesis( SMESH_subMesh* subMesh ) { - SMESH_HypoFilter ghsHypFilter( SMESH_HypoFilter::HasName( "GHS3D_Parameters" )); + SMESH_HypoFilter ghsHypFilter + ( SMESH_HypoFilter::HasName( GHS3DPlugin_Hypothesis::GetHypType() )); return (GHS3DPlugin_Hypothesis* ) subMesh->GetFather()->GetHypothesis( subMesh->GetSubShape(), ghsHypFilter, /*visitAncestors=*/true); } }; + + //================================================================================ + /*! + * \brief Sub-mesh event listener removing empty groups created due to "To make + * groups of domains". + */ + struct _GroupsOfDomainsRemover : public SMESH_subMeshEventListener + { + _GroupsOfDomainsRemover(): + SMESH_subMeshEventListener( /*isDeletable = */true, + "GHS3DPlugin_GHS3D::_GroupsOfDomainsRemover" ) {} + /*! + * \brief Treat events of the subMesh + */ + void ProcessEvent(const int event, + const int eventType, + SMESH_subMesh* subMesh, + SMESH_subMeshEventListenerData* data, + const SMESH_Hypothesis* hyp) + { + if (SMESH_subMesh::ALGO_EVENT == eventType && + !subMesh->GetAlgo() ) + { + removeEmptyGroupsOfDomains( subMesh->GetFather(), /*notEmptyAsWell=*/true ); + } + } + }; } //================================================================================ @@ -4392,3 +2911,42 @@ void GHS3DPlugin_GHS3D::SubmeshRestored(SMESH_subMesh* subMesh) } } } + +//================================================================================ +/*! + * \brief Sets an event listener removing empty groups created due to "To make + * groups of domains". + * \param subMesh - submesh where algo is set + * + * This method is called when a submesh gets HYP_OK algo_state. + * After being set, event listener is notified on each event of a submesh. + */ +//================================================================================ + +void GHS3DPlugin_GHS3D::SetEventListener(SMESH_subMesh* subMesh) +{ + subMesh->SetEventListener( new _GroupsOfDomainsRemover(), 0, subMesh ); +} + +//================================================================================ +/*! + * \brief If possible, returns progress of computation [0.,1.] + */ +//================================================================================ + +double GHS3DPlugin_GHS3D::GetProgress() const +{ + if ( _isLibUsed ) + { + // this->_progress is advanced by MG_Tetra_API according to messages from MG library + // but sharply. Advance it a bit to get smoother advancement. + GHS3DPlugin_GHS3D* me = const_cast( this ); + if ( _progress < 0.1 ) // the first message is at 10% + me->_progress = GetProgressByTic(); + else if ( _progress < 0.98 ) + me->_progress += _progressAdvance; + return _progress; + } + + return -1; +}