- MESSAGE("StdMeshers_MEFISTO_2D::Compute");
-
- if (_hypLengthFromEdges)
- _edgeLength = ComputeEdgeElementLength(aMesh, aShape);
-
- bool isOk = false;
- //const SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
- //SMESH_subMesh *theSubMesh = aMesh.GetSubMesh(aShape);
-
- const TopoDS_Face & FF = TopoDS::Face(aShape);
- bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
- TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
-
- Z nblf; //nombre de lignes fermees (enveloppe en tete)
- Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee
- R2 *uvslf = NULL;
- Z nbpti = 0; //nombre points internes futurs sommets de la triangulation
- R2 *uvpti = NULL;
-
- Z nbst;
- R2 *uvst = NULL;
- Z nbt;
- Z *nust = NULL;
- Z ierr = 0;
-
- Z nutysu = 1; // 1: il existe un fonction areteideale_()
- // Z nutysu=0; // 0: on utilise aretmx
- R aretmx = _edgeLength; // longueur max aretes future triangulation
-
- nblf = NumberOfWires(F);
-
- nudslf = new Z[1 + nblf];
- nudslf[0] = 0;
- int iw = 1;
- int nbpnt = 0;
-
- myOuterWire = BRepTools::OuterWire(F);
- nbpnt += NumberOfPoints(aMesh, myOuterWire);
- if ( nbpnt < 3 ) // ex: a circle with 2 segments
- return false;
- nudslf[iw++] = nbpnt;
-
- for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
- {
- const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
- if (!myOuterWire.IsSame(W))
- {
- nbpnt += NumberOfPoints(aMesh, W);
- nudslf[iw++] = nbpnt;
- }
- }
-
- // avoid passing same uv points for a vertex common to 2 wires
- TopTools_IndexedDataMapOfShapeListOfShape VWMap;
- if ( iw - 1 > 1 ) // nbofWires > 1
- TopExp::MapShapesAndAncestors( F , TopAbs_VERTEX, TopAbs_WIRE, VWMap );
-
- uvslf = new R2[nudslf[nblf]];
- int m = 0;
-
- double scalex, scaley;
- ComputeScaleOnFace(aMesh, F, scalex, scaley);
-
- map<int, const SMDS_MeshNode*> mefistoToDS; // correspondence mefisto index--> points IDNodes
- if ( !LoadPoints(aMesh, F, myOuterWire, uvslf, m,
- mefistoToDS, scalex, scaley, VWMap))
- return false;
-
- for (TopExp_Explorer exp(F, TopAbs_WIRE); exp.More(); exp.Next())
- {
- const TopoDS_Wire & W = TopoDS::Wire(exp.Current());
- if (!myOuterWire.IsSame(W))
- {
- if (! LoadPoints(aMesh, F, W, uvslf, m,
- mefistoToDS, scalex, scaley, VWMap ))
- return false;
- }
- }
-
- uvst = NULL;
- nust = NULL;
- aptrte(nutysu, aretmx,
- nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
-
- if (ierr == 0)
- {
- MESSAGE("... End Triangulation Generated Triangle Number " << nbt);
- MESSAGE(" Node Number " << nbst);
- StoreResult(aMesh, nbst, uvst, nbt, nust, F,
- faceIsForward, mefistoToDS, scalex, scaley);
- isOk = true;
- }
- else
- {
- MESSAGE("Error in Triangulation");
- isOk = false;
- }
- if (nudslf != NULL)
- delete[]nudslf;
- if (uvslf != NULL)
- delete[]uvslf;
- if (uvst != NULL)
- delete[]uvst;
- if (nust != NULL)
- delete[]nust;
- return isOk;
+ MESSAGE("StdMeshers_MEFISTO_2D::Compute");
+
+ TopoDS_Face F = TopoDS::Face(aShape.Oriented(TopAbs_FORWARD));
+
+ // helper builds quadratic mesh if necessary
+ SMESH_MesherHelper helper(aMesh);
+ _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, skipMediumNodes, problem, proxyMesh );
+ int nbWires = wires.size();
+ if ( problem && !problem->IsOK() ) return error( problem );
+ if ( nbWires == 0 ) return error( "Problem in StdMeshers_FaceSide::GetFaceWires()");
+ if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
+ return error(COMPERR_BAD_INPUT_MESH,
+ SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments());
+
+ // compute average edge length
+ if (!_hypMaxElementArea)
+ {
+ _edgeLength = 0;
+ int nbSegments = 0;
+ for ( int iW = 0; iW < nbWires; ++iW )
+ {
+ StdMeshers_FaceSidePtr wire = wires[ iW ];
+ _edgeLength += wire->Length();
+ nbSegments += wire->NbSegments();
+ }
+ if ( nbSegments )
+ _edgeLength /= nbSegments;
+ }
+
+ if (/*_hypLengthFromEdges &&*/ _edgeLength < DBL_MIN )
+ _edgeLength = 100;
+
+ Z nblf; //nombre de lignes fermees (enveloppe en tete)
+ Z *nudslf = NULL; //numero du dernier sommet de chaque ligne fermee
+ R2 *uvslf = NULL;
+ Z nbpti = 0; //nombre points internes futurs sommets de la triangulation
+ R2 *uvpti = NULL;
+
+ Z nbst;
+ R2 *uvst = NULL;
+ Z nbt;
+ Z *nust = NULL;
+ Z ierr = 0;
+
+ 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;
+
+ nudslf = new Z[1 + nblf];
+ nudslf[0] = 0;
+ int iw = 1;
+ int nbpnt = 0;
+
+ // count nb of input points
+ for ( int iW = 0; iW < nbWires; ++iW )
+ {
+ nbpnt += wires[iW]->NbPoints() - 1;
+ nudslf[iw++] = nbpnt;
+ }
+
+ uvslf = new R2[nudslf[nblf]];
+
+ double scalex, scaley;
+ ComputeScaleOnFace(aMesh, F, scalex, scaley);
+
+ // correspondence mefisto index --> Nodes
+ vector< const SMDS_MeshNode*> mefistoToDS(nbpnt, (const SMDS_MeshNode*)0);
+
+ bool isOk = false;
+
+ // fill input points UV
+ if ( LoadPoints(wires, uvslf, mefistoToDS, scalex, scaley) )
+ {
+ // Compute
+ aptrte(nutysu, aretmx,
+ nblf, nudslf, uvslf, nbpti, uvpti, nbst, uvst, nbt, nust, ierr);
+
+ if (ierr == 0)
+ {
+ MESSAGE("... End Triangulation Generated Triangle Number " << nbt);
+ MESSAGE(" Node Number " << nbst);
+ StoreResult(nbst, uvst, nbt, nust, mefistoToDS, scalex, scaley);
+ isOk = true;
+ }
+ else
+ {
+ error(ierr,"Error in Triangulation (aptrte())");
+ }
+ }
+ if (nudslf != NULL) delete[]nudslf;
+ if (uvslf != NULL) delete[]uvslf;
+ if (uvst != NULL) delete[]uvst;
+ if (nust != NULL) delete[]nust;
+
+ 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<int> 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<int> aResVec(SMDSEntity_Last);
+ for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ aResMap.insert(std::make_pair(sm,aResVec));
+ SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+ 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<int> aVec(SMDSEntity_Last);
+ for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
+ if(IsQuadratic) {
+ aVec[SMDSEntity_Quad_Triangle] = nbFaces;
+ aVec[SMDSEntity_Node] = (int)( nbNodes + nbFaces*3 - (NbSeg-1) );
+ }
+ else {
+ aVec[SMDSEntity_Node] = nbNodes;
+ aVec[SMDSEntity_Triangle] = nbFaces;
+ }
+ SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+ aResMap.insert(std::make_pair(sm,aVec));
+
+ return true;