X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_MEFISTO_2D.cxx;h=75db563136dd7a0ab154bcb6785bd8771646fc61;hp=7b76b35455e95acc1fba4b4cfcb831cb2a266415;hb=b0a908c0d20341651771d0249fb10882f54b2aad;hpb=79b1ac2b6df9117f16f11d444b1f165d477a1813 diff --git a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx index 7b76b3545..75db56313 100644 --- a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx +++ b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx @@ -1,59 +1,59 @@ -// SMESH SMESH : implementaion of SMESH idl descriptions +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE // -// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, -// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS -// -// 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. -// -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -// -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// 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, 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 +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // + +// SMESH SMESH : implementaion of SMESH idl descriptions // File : StdMeshers_MEFISTO_2D.cxx // Moved here from SMESH_MEFISTO_2D.cxx // Author : Paul RASCLE, EDF // Module : SMESH -// $Header$ - +// #include "StdMeshers_MEFISTO_2D.hxx" +#include "SMDS_EdgePosition.hxx" +#include "SMDS_MeshElement.hxx" +#include "SMDS_MeshNode.hxx" +#include "SMESH_Comment.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Mesh.hxx" -#include "SMESH_subMesh.hxx" -#include "SMESH_Block.hxx" #include "SMESH_MesherHelper.hxx" -#include "SMESH_Comment.hxx" - +#include "SMESH_subMesh.hxx" #include "StdMeshers_FaceSide.hxx" -#include "StdMeshers_MaxElementArea.hxx" #include "StdMeshers_LengthFromEdges.hxx" +#include "StdMeshers_MaxElementArea.hxx" +#include "StdMeshers_ViscousLayers2D.hxx" + +#include "utilities.h" #include "Rn.h" #include "aptrte.h" -#include "SMDS_MeshElement.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_EdgePosition.hxx" -#include "SMDS_FacePosition.hxx" - -#include "utilities.h" - +#include #include #include +#include #include +#include #include +#include #include #include #include @@ -63,10 +63,15 @@ #include #include #include +#include #include using namespace std; +#ifdef _DEBUG_ +//#define DUMP_POINTS // to print coordinates of MEFISTO input +#endif + //============================================================================= /*! * @@ -81,12 +86,13 @@ StdMeshers_MEFISTO_2D::StdMeshers_MEFISTO_2D(int hypId, int studyId, SMESH_Gen * _shapeType = (1 << TopAbs_FACE); _compatibleHypothesis.push_back("MaxElementArea"); _compatibleHypothesis.push_back("LengthFromEdges"); + _compatibleHypothesis.push_back("ViscousLayers2D"); _edgeLength = 0; _maxElementArea = 0; _hypMaxElementArea = NULL; _hypLengthFromEdges = NULL; - myTool = 0; + _helper = 0; } //============================================================================= @@ -113,6 +119,8 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis { _hypMaxElementArea = NULL; _hypLengthFromEdges = NULL; + _edgeLength = 0; + _maxElementArea = 0; list ::const_iterator itl; const SMESHDS_Hypothesis *theHyp; @@ -121,8 +129,8 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis int nbHyp = hyps.size(); if (!nbHyp) { - aStatus = SMESH_Hypothesis::HYP_MISSING; - return false; // can't work with no hypothesis + aStatus = SMESH_Hypothesis::HYP_OK; //SMESH_Hypothesis::HYP_MISSING; + return true; // (PAL13464) can work with no hypothesis, LengthFromEdges is default one } itl = hyps.begin(); @@ -137,7 +145,6 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis _hypMaxElementArea = static_cast(theHyp); ASSERT(_hypMaxElementArea); _maxElementArea = _hypMaxElementArea->GetMaxArea(); - _edgeLength = 0; isOk = true; aStatus = SMESH_Hypothesis::HYP_OK; } @@ -146,8 +153,6 @@ bool StdMeshers_MEFISTO_2D::CheckHypothesis { _hypLengthFromEdges = static_cast(theHyp); ASSERT(_hypLengthFromEdges); - _edgeLength = 0; - _maxElementArea = 0; isOk = true; aStatus = SMESH_Hypothesis::HYP_OK; } @@ -186,13 +191,19 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh // helper builds quadratic mesh if necessary SMESH_MesherHelper helper(aMesh); - myTool = &helper; - _quadraticMesh = myTool->IsQuadraticSubMesh(aShape); - const bool ignoreMediumNodes = _quadraticMesh; + _helper = &helper; + _quadraticMesh = _helper->IsQuadraticSubMesh(aShape); + const bool skipMediumNodes = _quadraticMesh; + + // build viscous layers if required + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F ); + if ( !proxyMesh ) + return false; // get all edges of a face TError problem; - TWireVector wires = StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, problem ); + TWireVector wires = + StdMeshers_FaceSide::GetFaceWires( F, aMesh, skipMediumNodes, problem, proxyMesh ); int nbWires = wires.size(); if ( problem && !problem->IsOK() ) return error( problem ); if ( nbWires == 0 ) return error( "Problem in StdMeshers_FaceSide::GetFaceWires()"); @@ -201,7 +212,7 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh SMESH_Comment("Too few segments: ")<NbSegments()); // compute average edge length - if (_hypLengthFromEdges) + if (!_hypMaxElementArea) { _edgeLength = 0; int nbSegments = 0; @@ -215,7 +226,7 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh _edgeLength /= nbSegments; } - if (_hypLengthFromEdges && _edgeLength < DBL_MIN ) + if (/*_hypLengthFromEdges &&*/ _edgeLength < DBL_MIN ) _edgeLength = 100; Z nblf; //nombre de lignes fermees (enveloppe en tete) @@ -233,6 +244,8 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh Z nutysu = 1; // 1: il existe un fonction areteideale_() // Z nutysu=0; // 0: on utilise aretmx R aretmx = _edgeLength; // longueur max aretes future triangulation + if ( _hypMaxElementArea ) + aretmx *= 1.5; nblf = nbWires; @@ -285,6 +298,90 @@ bool StdMeshers_MEFISTO_2D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh return isOk; } + +//============================================================================= +/*! + * + */ +//============================================================================= + +bool StdMeshers_MEFISTO_2D::Evaluate(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + MapShapeNbElems& aResMap) +{ + MESSAGE("StdMeshers_MEFISTO_2D::Evaluate"); + + TopoDS_Face F = TopoDS::Face(aShape.Oriented(TopAbs_FORWARD)); + + double aLen = 0.0; + int NbSeg = 0; + bool IsQuadratic = false; + bool IsFirst = true; + TopExp_Explorer exp(F,TopAbs_EDGE); + for(; exp.More(); exp.Next()) { + TopoDS_Edge E = TopoDS::Edge(exp.Current()); + MapShapeNbElemsItr anIt = aResMap.find( aMesh.GetSubMesh(E) ); + if( anIt == aResMap.end() ) continue; + std::vector aVec = (*anIt).second; + int nbe = Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]); + NbSeg += nbe; + if(IsFirst) { + IsQuadratic = ( aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge] ); + IsFirst = false; + } + double a,b; + TopLoc_Location L; + Handle(Geom_Curve) C = BRep_Tool::Curve(E,L,a,b); + gp_Pnt P1; + C->D0(a,P1); + double dp = (b-a)/nbe; + for(int i=1; i<=nbe; i++) { + gp_Pnt P2; + C->D0(a+i*dp,P2); + aLen += P1.Distance(P2); + P1 = P2; + } + } + if(NbSeg<1) { + std::vector aResVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; iGetComputeError(); + smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED, + "Submesh can not be evaluated",this)); + return false; + } + aLen = aLen/NbSeg; // middle length + + _edgeLength = Precision::Infinite(); + double tmpLength = Min( _edgeLength, aLen ); + + GProp_GProps G; + BRepGProp::SurfaceProperties(aShape,G); + double anArea = G.Mass(); + + int nbFaces = Precision::IsInfinite( tmpLength ) ? 0 : + (int)( anArea/(tmpLength*tmpLength*sqrt(3.)/4) ); + int nbNodes = (int) ( nbFaces*3 - (NbSeg-1)*2 ) / 6; + + std::vector aVec(SMDSEntity_Last); + for(int i=SMDSEntity_Node; i(node->GetPosition().get()); + static_cast(node->GetPosition()); double u = epos->GetUParameter(); if ( u < umin ) umin = u; @@ -484,7 +581,7 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector & wires, TopTools_IndexedDataMapOfShapeListOfShape VWMap; if ( wires.size() > 1 ) { - F = TopoDS::Face( myTool->GetSubShape() ); + F = TopoDS::Face( _helper->GetSubShape() ); TopExp::MapShapesAndAncestors( F, TopAbs_VERTEX, TopAbs_WIRE, VWMap ); int nbVertices = 0; for ( int iW = 0; iW < wires.size(); ++iW ) @@ -494,20 +591,21 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector & wires, } int m = 0; - list< int > mOnVertex; for ( int iW = 0; iW < wires.size(); ++iW ) { const vector& uvPtVec = wires[ iW ]->GetUVPtStruct(); if ( uvPtVec.size() != wires[ iW ]->NbPoints() ) { return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Unexpected nb of points on wire ") - << iW << uvPtVec.size()<<" != "<NbPoints()); + << iW << ": " << uvPtVec.size()<<" != "<NbPoints() + << ", probably because of invalid node parameters on geom edges"); } if ( m + uvPtVec.size()-1 > mefistoToDS.size() ) { MESSAGE("Wrong mefistoToDS.size: "< mOnVertex; vector::const_iterator uvPt = uvPtVec.begin(); for ( ++uvPt; uvPt != uvPtVec.end(); ++uvPt ) { @@ -516,8 +614,31 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector & wires, // set UV uvslf[m].x = uvPt->u * scalex; uvslf[m].y = uvPt->v * scaley; - if ( uvPt->node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) + switch ( uvPt->node->GetPosition()->GetTypeOfPosition()) + { + case SMDS_TOP_VERTEX: mOnVertex.push_back( m ); + break; + case SMDS_TOP_EDGE: + // In order to detect degenerated faces easily, we replace + // nodes on a degenerated edge by node on the vertex of that edge + if ( _helper->IsDegenShape( uvPt->node->getshapeId() )) + { + int edgeID = uvPt->node->getshapeId(); + SMESH_subMesh* edgeSM = _helper->GetMesh()->GetSubMeshContaining( edgeID ); + SMESH_subMeshIteratorPtr smIt = edgeSM->getDependsOnIterator( /*includeSelf=*/0, + /*complexShapeFirst=*/0); + if ( smIt->more() ) + { + SMESH_subMesh* vertexSM = smIt->next(); + SMDS_NodeIteratorPtr nIt = vertexSM->GetSubMeshDS()->GetNodes(); + if ( nIt->more() ) + mefistoToDS[m] = nIt->next(); + } + } + break; + default:; + } m++; } @@ -528,9 +649,9 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector & wires, int m = *mIt; if ( iW && !VWMap.IsEmpty()) { // except outer wire // avoid passing same uv point for a vertex common to 2 wires - int vID = mefistoToDS[m]->GetPosition()->GetShapeId(); - TopoDS_Vertex V = TopoDS::Vertex( myTool->GetMeshDS()->IndexToShape( vID )); - if ( fixCommonVertexUV( uvslf[m], V, F, VWMap, *myTool->GetMesh(), + int vID = mefistoToDS[m]->getshapeId(); + TopoDS_Vertex V = TopoDS::Vertex( _helper->GetMeshDS()->IndexToShape( vID )); + if ( fixCommonVertexUV( uvslf[m], V, F, VWMap, *_helper->GetMesh(), scalex, scaley, _quadraticMesh )) { myNodesOnCommonV.push_back( mefistoToDS[m] ); continue; @@ -545,6 +666,13 @@ bool StdMeshers_MEFISTO_2D::LoadPoints(TWireVector & wires, } } +#ifdef DUMP_POINTS + cout << "MEFISTO INPUT************" << endl; + for ( int i =0; i < m; ++i ) + cout << i << ": \t" << uvslf[i].x << ", " << uvslf[i].y + << " Node " << mefistoToDS[i]->GetID()<< endl; +#endif + return true; } @@ -642,6 +770,23 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, ASSERT(scaley); } +// namespace +// { +// bool isDegenTria( const SMDS_MeshNode * nn[3] ) +// { +// SMESH_TNodeXYZ p1( nn[0] ); +// SMESH_TNodeXYZ p2( nn[1] ); +// SMESH_TNodeXYZ p3( nn[2] ); +// gp_XYZ vec1 = p2 - p1; +// gp_XYZ vec2 = p3 - p1; +// gp_XYZ cross = vec1 ^ vec2; +// const double eps = 1e-100; +// return ( fabs( cross.X() ) < eps && +// fabs( cross.Y() ) < eps && +// fabs( cross.Z() ) < eps ); +// } +// } + //============================================================================= /*! * @@ -652,12 +797,13 @@ void StdMeshers_MEFISTO_2D::StoreResult(Z nbst, R2 * uvst, Z nbt, Z * nust, vector< const SMDS_MeshNode*>&mefistoToDS, double scalex, double scaley) { - SMESHDS_Mesh * meshDS = myTool->GetMeshDS(); - int faceID = myTool->GetSubShapeID(); + _helper->SetElementsOnShape( true ); - TopoDS_Face F = TopoDS::Face( myTool->GetSubShape() ); + TopoDS_Face F = TopoDS::Face( _helper->GetSubShape() ); Handle(Geom_Surface) S = BRep_Tool::Surface( F ); + //const size_t nbInputNodes = mefistoToDS.size(); + Z n = mefistoToDS.size(); // nb input points mefistoToDS.resize( nbst ); for ( ; n < nbst; n++) @@ -668,12 +814,7 @@ void StdMeshers_MEFISTO_2D::StoreResult(Z nbst, R2 * uvst, Z nbt, Z * nust, double v = uvst[n][1] / scaley; gp_Pnt P = S->Value(u, v); - SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); - meshDS->SetNodeOnFace(node, faceID, u, v); - - //MESSAGE(P.X()<<" "<AddNode( P.X(), P.Y(), P.Z(), 0, u, v ); } } @@ -682,22 +823,32 @@ void StdMeshers_MEFISTO_2D::StoreResult(Z nbst, R2 * uvst, Z nbt, Z * nust, // triangle points must be in trigonometric order if face is Forward // else they must be put clockwise - bool triangleIsWellOriented = ( F.Orientation() == TopAbs_FORWARD ); + int i1 = 1, i2 = 2; + if ( F.Orientation() != TopAbs_FORWARD ) + std::swap( i1, i2 ); + const SMDS_MeshNode * nn[3]; for (n = 1; n <= nbt; n++) { - const SMDS_MeshNode * n1 = mefistoToDS[ nust[m++] - 1 ]; - const SMDS_MeshNode * n2 = mefistoToDS[ nust[m++] - 1 ]; - const SMDS_MeshNode * n3 = mefistoToDS[ nust[m++] - 1 ]; + // const bool allNodesAreOld = ( nust[m + 0] <= nbInputNodes && + // nust[m + 1] <= nbInputNodes && + // nust[m + 2] <= nbInputNodes ); + nn[ 0 ] = mefistoToDS[ nust[m++] - 1 ]; + nn[ 1 ] = mefistoToDS[ nust[m++] - 1 ]; + nn[ 2 ] = mefistoToDS[ nust[m++] - 1 ]; + m++; - SMDS_MeshElement * elt; - if (triangleIsWellOriented) - elt = myTool->AddFace(n1, n2, n3); - else - elt = myTool->AddFace(n1, n3, n2); + // avoid creating degenetrated faces + bool isDegen = ( _helper->HasDegeneratedEdges() && + ( nn[0] == nn[1] || nn[1] == nn[2] || nn[2] == nn[0] )); - meshDS->SetMeshElementOnShape(elt, faceID); - m++; + // It was an attemp to fix a problem of a zero area face whose all nodes + // are on one staight EDGE. But omitting this face makes a hole in the mesh :( + // if ( !isDegen && allNodesAreOld ) + // isDegen = isDegenTria( nn ); + + if ( !isDegen ) + _helper->AddFace( nn[0], nn[i1], nn[i2] ); } // remove bad elements built on vertices shared by wires @@ -717,7 +868,7 @@ void StdMeshers_MEFISTO_2D::StoreResult(Z nbst, R2 * uvst, Z nbt, Z * nust, nbSame++; if (nbSame > 1) { MESSAGE( "RM bad element " << elem->GetID()); - meshDS->RemoveElement( elem ); + _helper->GetMeshDS()->RemoveElement( elem ); } } }