1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 // File : StdMeshers_Quadrangle_2D.cxx
24 // Author : Paul RASCLE, EDF
27 #include "StdMeshers_Quadrangle_2D.hxx"
29 #include "SMDS_EdgePosition.hxx"
30 #include "SMDS_FacePosition.hxx"
31 #include "SMDS_MeshElement.hxx"
32 #include "SMDS_MeshNode.hxx"
33 #include "SMESH_Block.hxx"
34 #include "SMESH_Comment.hxx"
35 #include "SMESH_Gen.hxx"
36 #include "SMESH_Mesh.hxx"
37 #include "SMESH_MesherHelper.hxx"
38 #include "SMESH_subMesh.hxx"
39 #include "StdMeshers_FaceSide.hxx"
40 #include "StdMeshers_QuadrangleParams.hxx"
41 #include "StdMeshers_ViscousLayers2D.hxx"
43 #include <BRep_Tool.hxx>
44 #include <Geom_Surface.hxx>
45 #include <NCollection_DefineArray2.hxx>
46 #include <Precision.hxx>
47 #include <TColStd_SequenceOfReal.hxx>
48 #include <TColStd_SequenceOfInteger.hxx>
49 #include <TColgp_SequenceOfXY.hxx>
51 #include <TopExp_Explorer.hxx>
52 #include <TopTools_ListIteratorOfListOfShape.hxx>
53 #include <TopTools_MapOfShape.hxx>
56 #include "utilities.h"
57 #include "Utils_ExceptHandlers.hxx"
59 #ifndef StdMeshers_Array2OfNode_HeaderFile
60 #define StdMeshers_Array2OfNode_HeaderFile
61 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
62 DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
63 DEFINE_ARRAY2(StdMeshers_Array2OfNode,
64 StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
70 typedef SMESH_Comment TComm;
72 //=============================================================================
76 //=============================================================================
78 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
80 : SMESH_2D_Algo(hypId, studyId, gen),
81 myQuadranglePreference(false),
82 myTrianglePreference(false),
85 myQuadType(QUAD_STANDARD),
88 MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
89 _name = "Quadrangle_2D";
90 _shapeType = (1 << TopAbs_FACE);
91 _compatibleHypothesis.push_back("QuadrangleParams");
92 _compatibleHypothesis.push_back("QuadranglePreference");
93 _compatibleHypothesis.push_back("TrianglePreference");
94 _compatibleHypothesis.push_back("ViscousLayers2D");
97 //=============================================================================
101 //=============================================================================
103 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
105 MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
108 //=============================================================================
112 //=============================================================================
114 bool StdMeshers_Quadrangle_2D::CheckHypothesis
116 const TopoDS_Shape& aShape,
117 SMESH_Hypothesis::Hypothesis_Status& aStatus)
120 aStatus = SMESH_Hypothesis::HYP_OK;
122 const list <const SMESHDS_Hypothesis * >& hyps =
123 GetUsedHypothesis(aMesh, aShape, false);
124 const SMESHDS_Hypothesis * aHyp = 0;
127 myQuadType = QUAD_STANDARD;
128 myQuadranglePreference = false;
129 myTrianglePreference = false;
130 myQuadStruct.reset();
132 bool isFirstParams = true;
134 // First assigned hypothesis (if any) is processed now
135 if (hyps.size() > 0) {
137 if (strcmp("QuadrangleParams", aHyp->GetName()) == 0) {
138 const StdMeshers_QuadrangleParams* aHyp1 =
139 (const StdMeshers_QuadrangleParams*)aHyp;
140 myTriaVertexID = aHyp1->GetTriaVertex();
141 myQuadType = aHyp1->GetQuadType();
142 if (myQuadType == QUAD_QUADRANGLE_PREF ||
143 myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
144 myQuadranglePreference = true;
145 else if (myQuadType == QUAD_TRIANGLE_PREF)
146 myTrianglePreference = true;
148 else if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) {
149 isFirstParams = false;
150 myQuadranglePreference = true;
152 else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){
153 isFirstParams = false;
154 myTrianglePreference = true;
157 isFirstParams = false;
161 // Second(last) assigned hypothesis (if any) is processed now
162 if (hyps.size() > 1) {
165 if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) {
166 myQuadranglePreference = true;
167 myTrianglePreference = false;
168 myQuadType = QUAD_STANDARD;
170 else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){
171 myQuadranglePreference = false;
172 myTrianglePreference = true;
173 myQuadType = QUAD_STANDARD;
177 const StdMeshers_QuadrangleParams* aHyp2 =
178 (const StdMeshers_QuadrangleParams*)aHyp;
179 myTriaVertexID = aHyp2->GetTriaVertex();
181 if (!myQuadranglePreference && !myTrianglePreference) { // priority of hypos
182 myQuadType = aHyp2->GetQuadType();
183 if (myQuadType == QUAD_QUADRANGLE_PREF ||
184 myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
185 myQuadranglePreference = true;
186 else if (myQuadType == QUAD_TRIANGLE_PREF)
187 myTrianglePreference = true;
195 //=============================================================================
199 //=============================================================================
201 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
202 const TopoDS_Shape& aShape)
204 const TopoDS_Face& F = TopoDS::Face(aShape);
205 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
207 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
208 aMesh.GetSubMesh(aShape);
210 SMESH_MesherHelper helper (aMesh);
213 myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
217 _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
218 myNeedSmooth = false;
220 FaceQuadStruct::Ptr quad = CheckNbEdges(aMesh, aShape);
225 if (myQuadranglePreference) {
226 int n1 = quad->side[0]->NbPoints();
227 int n2 = quad->side[1]->NbPoints();
228 int n3 = quad->side[2]->NbPoints();
229 int n4 = quad->side[3]->NbPoints();
230 int nfull = n1+n2+n3+n4;
233 if (nfull == ntmp && ((n1 != n3) || (n2 != n4))) {
234 // special path for using only quandrangle faces
235 bool ok = ComputeQuadPref(aMesh, aShape, quad);
236 if ( ok && myNeedSmooth )
241 else if (myQuadType == QUAD_REDUCED) {
242 int n1 = quad->side[0]->NbPoints();
243 int n2 = quad->side[1]->NbPoints();
244 int n3 = quad->side[2]->NbPoints();
245 int n4 = quad->side[3]->NbPoints();
248 int n13tmp = n13/2; n13tmp = n13tmp*2;
249 int n24tmp = n24/2; n24tmp = n24tmp*2;
250 if ((n1 == n3 && n2 != n4 && n24tmp == n24) ||
251 (n2 == n4 && n1 != n3 && n13tmp == n13)) {
252 bool ok = ComputeReduced(aMesh, aShape, quad);
253 if ( ok && myNeedSmooth )
257 if ( n1 != n3 && n2 != n4 )
258 error( COMPERR_WARNING,
259 "To use 'Reduced' transition, "
260 "two opposite sides should have same number of segments, "
261 "but actual number of segments is different on all sides. "
262 "'Standard' transion has been used.");
264 error( COMPERR_WARNING,
265 "To use 'Reduced' transition, "
266 "two opposite sides should have an even difference in number of segments. "
267 "'Standard' transion has been used.");
270 // set normalized grid on unit square in parametric domain
272 if (!SetNormalizedGrid(aMesh, aShape, quad))
275 // --- compute 3D values on points, store points & quadrangles
277 int nbdown = quad->side[0]->NbPoints();
278 int nbup = quad->side[2]->NbPoints();
280 int nbright = quad->side[1]->NbPoints();
281 int nbleft = quad->side[3]->NbPoints();
283 int nbhoriz = Min(nbdown, nbup);
284 int nbvertic = Min(nbright, nbleft);
286 // internal mesh nodes
287 int i, j, geomFaceID = meshDS->ShapeToIndex(F);
288 for (i = 1; i < nbhoriz - 1; i++) {
289 for (j = 1; j < nbvertic - 1; j++) {
290 int ij = j * nbhoriz + i;
291 double u = quad->uv_grid[ij].u;
292 double v = quad->uv_grid[ij].v;
293 gp_Pnt P = S->Value(u, v);
294 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
295 meshDS->SetNodeOnFace(node, geomFaceID, u, v);
296 quad->uv_grid[ij].node = node;
303 // --.--.--.--.--.-- nbvertic
309 // ---.----.----.--- 0
310 // 0 > > > > > > > > nbhoriz
316 int iup = nbhoriz - 1;
317 if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
320 int jup = nbvertic - 1;
321 if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
323 // regular quadrangles
324 for (i = ilow; i < iup; i++) {
325 for (j = jlow; j < jup; j++) {
326 const SMDS_MeshNode *a, *b, *c, *d;
327 a = quad->uv_grid[j * nbhoriz + i ].node;
328 b = quad->uv_grid[j * nbhoriz + i + 1].node;
329 c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
330 d = quad->uv_grid[(j + 1) * nbhoriz + i ].node;
331 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
333 meshDS->SetMeshElementOnShape(face, geomFaceID);
338 const vector<UVPtStruct>& uv_e0 = quad->side[0]->GetUVPtStruct(true,0);
339 const vector<UVPtStruct>& uv_e1 = quad->side[1]->GetUVPtStruct(false,1);
340 const vector<UVPtStruct>& uv_e2 = quad->side[2]->GetUVPtStruct(true,1);
341 const vector<UVPtStruct>& uv_e3 = quad->side[3]->GetUVPtStruct(false,0);
343 if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
344 return error(COMPERR_BAD_INPUT_MESH);
346 double eps = Precision::Confusion();
348 // Boundary quadrangles
350 if (quad->isEdgeOut[0]) {
353 // |___|___|___|___|___|___|
355 // |___|___|___|___|___|___|
357 // |___|___|___|___|___|___| __ first row of the regular grid
358 // . . . . . . . . . __ down edge nodes
360 // >->->->->->->->->->->->-> -- direction of processing
362 int g = 0; // number of last processed node in the regular grid
364 // number of last node of the down edge to be processed
365 int stop = nbdown - 1;
366 // if right edge is out, we will stop at a node, previous to the last one
367 if (quad->isEdgeOut[1]) stop--;
369 // for each node of the down edge find nearest node
370 // in the first row of the regular grid and link them
371 for (i = 0; i < stop; i++) {
372 const SMDS_MeshNode *a, *b, *c, *d;
374 b = uv_e0[i + 1].node;
375 gp_Pnt pb (b->X(), b->Y(), b->Z());
377 // find node c in the regular grid, which will be linked with node b
380 // right bound reached, link with the rightmost node
382 c = quad->uv_grid[nbhoriz + iup].node;
385 // find in the grid node c, nearest to the b
386 double mind = RealLast();
387 for (int k = g; k <= iup; k++) {
389 const SMDS_MeshNode *nk;
390 if (k < ilow) // this can be, if left edge is out
391 nk = uv_e3[1].node; // get node from the left edge
393 nk = quad->uv_grid[nbhoriz + k].node; // get one of middle nodes
395 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
396 double dist = pb.Distance(pnk);
397 if (dist < mind - eps) {
407 if (near == g) { // make triangle
408 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
409 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
411 else { // make quadrangle
415 d = quad->uv_grid[nbhoriz + near - 1].node;
416 //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
418 if (!myTrianglePreference){
419 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
420 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
423 SplitQuad(meshDS, geomFaceID, a, b, c, d);
426 // if node d is not at position g - make additional triangles
428 for (int k = near - 1; k > g; k--) {
429 c = quad->uv_grid[nbhoriz + k].node;
433 d = quad->uv_grid[nbhoriz + k - 1].node;
434 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
435 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
442 if (quad->isEdgeOut[2]) {
445 // <-<-<-<-<-<-<-<-<-<-<-<-< -- direction of processing
447 // . . . . . . . . . __ up edge nodes
448 // ___ ___ ___ ___ ___ ___ __ first row of the regular grid
450 // |___|___|___|___|___|___|
452 // |___|___|___|___|___|___|
455 int g = nbhoriz - 1; // last processed node in the regular grid
458 // if left edge is out, we will stop at a second node
459 if (quad->isEdgeOut[3]) stop++;
461 // for each node of the up edge find nearest node
462 // in the first row of the regular grid and link them
463 for (i = nbup - 1; i > stop; i--) {
464 const SMDS_MeshNode *a, *b, *c, *d;
466 b = uv_e2[i - 1].node;
467 gp_Pnt pb (b->X(), b->Y(), b->Z());
469 // find node c in the grid, which will be linked with node b
471 if (i == stop + 1) { // left bound reached, link with the leftmost node
472 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
475 // find node c in the grid, nearest to the b
476 double mind = RealLast();
477 for (int k = g; k >= ilow; k--) {
478 const SMDS_MeshNode *nk;
480 nk = uv_e1[nbright - 2].node;
482 nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
483 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
484 double dist = pb.Distance(pnk);
485 if (dist < mind - eps) {
495 if (near == g) { // make triangle
496 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
497 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
499 else { // make quadrangle
501 d = uv_e1[nbright - 2].node;
503 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
504 //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
505 if (!myTrianglePreference){
506 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
507 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
510 SplitQuad(meshDS, geomFaceID, a, b, c, d);
513 if (near + 1 < g) { // if d not is at g - make additional triangles
514 for (int k = near + 1; k < g; k++) {
515 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
517 d = uv_e1[nbright - 2].node;
519 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
520 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
521 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
530 // right or left boundary quadrangles
531 if (quad->isEdgeOut[1]) {
532 // MESSAGE("right edge is out");
533 int g = 0; // last processed node in the grid
534 int stop = nbright - 1;
535 if (quad->isEdgeOut[2]) stop--;
536 for (i = 0; i < stop; i++) {
537 const SMDS_MeshNode *a, *b, *c, *d;
539 b = uv_e1[i + 1].node;
540 gp_Pnt pb (b->X(), b->Y(), b->Z());
542 // find node c in the grid, nearest to the b
544 if (i == stop - 1) { // up bondary reached
545 c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
548 double mind = RealLast();
549 for (int k = g; k <= jup; k++) {
550 const SMDS_MeshNode *nk;
552 nk = uv_e0[nbdown - 2].node;
554 nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
555 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
556 double dist = pb.Distance(pnk);
557 if (dist < mind - eps) {
567 if (near == g) { // make triangle
568 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
569 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
571 else { // make quadrangle
573 d = uv_e0[nbdown - 2].node;
575 d = quad->uv_grid[nbhoriz*near - 2].node;
576 //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
578 if (!myTrianglePreference){
579 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
580 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
583 SplitQuad(meshDS, geomFaceID, a, b, c, d);
586 if (near - 1 > g) { // if d not is at g - make additional triangles
587 for (int k = near - 1; k > g; k--) {
588 c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
590 d = uv_e0[nbdown - 2].node;
592 d = quad->uv_grid[nbhoriz*k - 2].node;
593 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
594 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
601 if (quad->isEdgeOut[3]) {
602 // MESSAGE("left edge is out");
603 int g = nbvertic - 1; // last processed node in the grid
605 if (quad->isEdgeOut[0]) stop++;
606 for (i = nbleft - 1; i > stop; i--) {
607 const SMDS_MeshNode *a, *b, *c, *d;
609 b = uv_e3[i - 1].node;
610 gp_Pnt pb (b->X(), b->Y(), b->Z());
612 // find node c in the grid, nearest to the b
614 if (i == stop + 1) { // down bondary reached
615 c = quad->uv_grid[nbhoriz*jlow + 1].node;
618 double mind = RealLast();
619 for (int k = g; k >= jlow; k--) {
620 const SMDS_MeshNode *nk;
624 nk = quad->uv_grid[nbhoriz*k + 1].node;
625 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
626 double dist = pb.Distance(pnk);
627 if (dist < mind - eps) {
637 if (near == g) { // make triangle
638 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
639 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
641 else { // make quadrangle
645 d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
646 //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
647 if (!myTrianglePreference){
648 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
649 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
652 SplitQuad(meshDS, geomFaceID, a, b, c, d);
655 if (near + 1 < g) { // if d not is at g - make additional triangles
656 for (int k = near + 1; k < g; k++) {
657 c = quad->uv_grid[nbhoriz*k + 1].node;
661 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
662 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
663 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
680 //=============================================================================
684 //=============================================================================
686 bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
687 const TopoDS_Shape& aShape,
688 MapShapeNbElems& aResMap)
691 aMesh.GetSubMesh(aShape);
693 std::vector<int> aNbNodes(4);
694 bool IsQuadratic = false;
695 if (!CheckNbEdgesForEvaluate(aMesh, aShape, aResMap, aNbNodes, IsQuadratic)) {
696 std::vector<int> aResVec(SMDSEntity_Last);
697 for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
698 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
699 aResMap.insert(std::make_pair(sm,aResVec));
700 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
701 smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
705 if (myQuadranglePreference) {
706 int n1 = aNbNodes[0];
707 int n2 = aNbNodes[1];
708 int n3 = aNbNodes[2];
709 int n4 = aNbNodes[3];
710 int nfull = n1+n2+n3+n4;
713 if (nfull==ntmp && ((n1!=n3) || (n2!=n4))) {
714 // special path for using only quandrangle faces
715 return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
720 int nbdown = aNbNodes[0];
721 int nbup = aNbNodes[2];
723 int nbright = aNbNodes[1];
724 int nbleft = aNbNodes[3];
726 int nbhoriz = Min(nbdown, nbup);
727 int nbvertic = Min(nbright, nbleft);
729 int dh = Max(nbdown, nbup) - nbhoriz;
730 int dv = Max(nbright, nbleft) - nbvertic;
737 int nbNodes = (nbhoriz-2)*(nbvertic-2);
738 //int nbFaces3 = dh + dv + kdh*(nbvertic-1)*2 + kdv*(nbhoriz-1)*2;
739 int nbFaces3 = dh + dv;
740 //if (kdh==1 && kdv==1) nbFaces3 -= 2;
741 //if (dh>0 && dv>0) nbFaces3 -= 2;
742 //int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv);
743 int nbFaces4 = (nbhoriz-1)*(nbvertic-1);
745 std::vector<int> aVec(SMDSEntity_Last);
746 for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
748 aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
749 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
750 int nbbndedges = nbdown + nbup + nbright + nbleft -4;
751 int nbintedges = (nbFaces4*4 + nbFaces3*3 - nbbndedges) / 2;
752 aVec[SMDSEntity_Node] = nbNodes + nbintedges;
753 if (aNbNodes.size()==5) {
754 aVec[SMDSEntity_Quad_Triangle] = nbFaces3 + aNbNodes[3] -1;
755 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4 - aNbNodes[3] +1;
759 aVec[SMDSEntity_Node] = nbNodes;
760 aVec[SMDSEntity_Triangle] = nbFaces3;
761 aVec[SMDSEntity_Quadrangle] = nbFaces4;
762 if (aNbNodes.size()==5) {
763 aVec[SMDSEntity_Triangle] = nbFaces3 + aNbNodes[3] - 1;
764 aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1;
767 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
768 aResMap.insert(std::make_pair(sm,aVec));
774 //================================================================================
776 * \brief Return true if only two given edges meat at their common vertex
778 //================================================================================
780 static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
781 const TopoDS_Edge& e2,
785 if (!TopExp::CommonVertex(e1, e2, v))
787 TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v));
788 for (; ancestIt.More() ; ancestIt.Next())
789 if (ancestIt.Value().ShapeType() == TopAbs_EDGE)
790 if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value()))
795 //=============================================================================
799 //=============================================================================
801 FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh,
802 const TopoDS_Shape & aShape)
804 if ( myQuadStruct && myQuadStruct->face.IsSame( aShape ))
807 TopoDS_Face F = TopoDS::Face(aShape);
808 if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
809 const bool ignoreMediumNodes = _quadraticMesh;
811 // verify 1 wire only, with 4 edges
812 list< TopoDS_Edge > edges;
813 list< int > nbEdgesInWire;
814 int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
816 error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
817 return FaceQuadStruct::Ptr();
819 FaceQuadStruct::Ptr quad( new FaceQuadStruct );
821 quad->side.reserve(nbEdgesInWire.front());
825 list< TopoDS_Edge >::iterator edgeIt = edges.begin();
826 if (nbEdgesInWire.front() == 3) // exactly 3 edges
828 SMESH_Comment comment;
829 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
830 if (myTriaVertexID < 1)
832 comment << "No Base vertex parameter provided for a trilateral geometrical face";
836 TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
838 TopoDS_Edge E1,E2,E3;
839 for (; edgeIt != edges.end(); ++edgeIt) {
840 TopoDS_Edge E = *edgeIt;
841 TopoDS_Vertex VF, VL;
842 TopExp::Vertices(E, VF, VL, true);
845 else if (VL.IsSame(V))
850 if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull())
852 quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true,
853 ignoreMediumNodes, myProxyMesh));
854 quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true,
855 ignoreMediumNodes, myProxyMesh));
856 quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,
857 ignoreMediumNodes, myProxyMesh));
858 const vector<UVPtStruct>& UVPSleft = quad->side[0]->GetUVPtStruct(true,0);
859 /* vector<UVPtStruct>& UVPStop = */quad->side[1]->GetUVPtStruct(false,1);
860 /* vector<UVPtStruct>& UVPSright = */quad->side[2]->GetUVPtStruct(true,1);
861 const SMDS_MeshNode* aNode = UVPSleft[0].node;
862 gp_Pnt2d aPnt2d(UVPSleft[0].u, UVPSleft[0].v);
863 quad->side.push_back(new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]));
867 comment << "Invalid Base vertex parameter: " << myTriaVertexID << " is not among [";
868 TopTools_MapOfShape vMap;
869 for (TopExp_Explorer v(aShape, TopAbs_VERTEX); v.More(); v.Next())
870 if (vMap.Add(v.Current()))
871 comment << meshDS->ShapeToIndex(v.Current()) << (vMap.Extent()==3 ? "]" : ", ");
877 else if (nbEdgesInWire.front() == 4) // exactly 4 edges
879 for (; edgeIt != edges.end(); ++edgeIt, nbSides++)
880 quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh, nbSides < QUAD_TOP_SIDE,
881 ignoreMediumNodes, myProxyMesh));
883 else if (nbEdgesInWire.front() > 4) // more than 4 edges - try to unite some
885 list< TopoDS_Edge > sideEdges;
886 vector< int > degenSides;
887 while (!edges.empty()) {
889 sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
890 bool sameSide = true;
891 while (!edges.empty() && sameSide) {
892 sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
894 sideEdges.splice(sideEdges.end(), edges, edges.begin());
896 if (nbSides == 0) { // go backward from the first edge
898 while (!edges.empty() && sameSide) {
899 sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
901 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
904 if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() ))
905 degenSides.push_back( nbSides );
907 quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh, nbSides < QUAD_TOP_SIDE,
908 ignoreMediumNodes, myProxyMesh));
911 if ( !degenSides.empty() && nbSides - degenSides.size() == 4 )
914 for ( unsigned i = QUAD_TOP_SIDE; i < quad->side.size(); ++i )
915 quad->side[i]->Reverse();
917 for ( int i = degenSides.size()-1; i > -1; --i )
919 StdMeshers_FaceSide* degenSide = quad->side[ degenSides[ i ]];
921 quad->side.erase( quad->side.begin() + degenSides[ i ] );
923 for ( unsigned i = QUAD_TOP_SIDE; i < quad->side.size(); ++i )
924 quad->side[i]->Reverse();
926 nbSides -= degenSides.size();
928 // issue 20222. Try to unite only edges shared by two same faces
931 quad.reset( new FaceQuadStruct );
932 quad->side.reserve(nbEdgesInWire.front());
935 SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
936 while (!edges.empty()) {
938 sideEdges.splice(sideEdges.end(), edges, edges.begin());
939 bool sameSide = true;
940 while (!edges.empty() && sameSide) {
942 SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
943 twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
945 sideEdges.splice(sideEdges.end(), edges, edges.begin());
947 if (nbSides == 0) { // go backward from the first edge
949 while (!edges.empty() && sameSide) {
951 SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
952 twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
954 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
957 quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
958 nbSides < QUAD_TOP_SIDE,
959 ignoreMediumNodes, myProxyMesh));
966 MESSAGE ("StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n");
967 for (int i = 0; i < nbSides; ++i) {
969 for (int e = 0; e < quad->side[i]->NbEdges(); ++e)
970 MESSAGE (myHelper->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " ");
975 nbSides = nbEdgesInWire.front();
976 error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
984 //=============================================================================
988 //=============================================================================
990 bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
991 const TopoDS_Shape & aShape,
992 MapShapeNbElems& aResMap,
993 std::vector<int>& aNbNodes,
997 const TopoDS_Face & F = TopoDS::Face(aShape);
999 // verify 1 wire only, with 4 edges
1000 list< TopoDS_Edge > edges;
1001 list< int > nbEdgesInWire;
1002 int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
1010 list< TopoDS_Edge >::iterator edgeIt = edges.begin();
1011 SMESH_subMesh * sm = aMesh.GetSubMesh(*edgeIt);
1012 MapShapeNbElemsItr anIt = aResMap.find(sm);
1013 if (anIt==aResMap.end()) {
1016 std::vector<int> aVec = (*anIt).second;
1017 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1018 if (nbEdgesInWire.front() == 3) { // exactly 3 edges
1019 if (myTriaVertexID>0) {
1020 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
1021 TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
1023 TopoDS_Edge E1,E2,E3;
1024 for (; edgeIt != edges.end(); ++edgeIt) {
1025 TopoDS_Edge E = TopoDS::Edge(*edgeIt);
1026 TopoDS_Vertex VF, VL;
1027 TopExp::Vertices(E, VF, VL, true);
1030 else if (VL.IsSame(V))
1035 SMESH_subMesh * sm = aMesh.GetSubMesh(E1);
1036 MapShapeNbElemsItr anIt = aResMap.find(sm);
1037 if (anIt==aResMap.end()) return false;
1038 std::vector<int> aVec = (*anIt).second;
1040 aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1042 aNbNodes[0] = aVec[SMDSEntity_Node] + 2;
1043 sm = aMesh.GetSubMesh(E2);
1044 anIt = aResMap.find(sm);
1045 if (anIt==aResMap.end()) return false;
1046 aVec = (*anIt).second;
1048 aNbNodes[1] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1050 aNbNodes[1] = aVec[SMDSEntity_Node] + 2;
1051 sm = aMesh.GetSubMesh(E3);
1052 anIt = aResMap.find(sm);
1053 if (anIt==aResMap.end()) return false;
1054 aVec = (*anIt).second;
1056 aNbNodes[2] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1058 aNbNodes[2] = aVec[SMDSEntity_Node] + 2;
1059 aNbNodes[3] = aNbNodes[1];
1065 if (nbEdgesInWire.front() == 4) { // exactly 4 edges
1066 for (; edgeIt != edges.end(); edgeIt++) {
1067 SMESH_subMesh * sm = aMesh.GetSubMesh(*edgeIt);
1068 MapShapeNbElemsItr anIt = aResMap.find(sm);
1069 if (anIt==aResMap.end()) {
1072 std::vector<int> aVec = (*anIt).second;
1074 aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2;
1076 aNbNodes[nbSides] = aVec[SMDSEntity_Node] + 2;
1080 else if (nbEdgesInWire.front() > 4) { // more than 4 edges - try to unite some
1081 list< TopoDS_Edge > sideEdges;
1082 while (!edges.empty()) {
1084 sideEdges.splice(sideEdges.end(), edges, edges.begin()); // edges.front() -> sideEdges.end()
1085 bool sameSide = true;
1086 while (!edges.empty() && sameSide) {
1087 sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
1089 sideEdges.splice(sideEdges.end(), edges, edges.begin());
1091 if (nbSides == 0) { // go backward from the first edge
1093 while (!edges.empty() && sameSide) {
1094 sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
1096 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
1099 list<TopoDS_Edge>::iterator ite = sideEdges.begin();
1100 aNbNodes[nbSides] = 1;
1101 for (; ite!=sideEdges.end(); ite++) {
1102 SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
1103 MapShapeNbElemsItr anIt = aResMap.find(sm);
1104 if (anIt==aResMap.end()) {
1107 std::vector<int> aVec = (*anIt).second;
1109 aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
1111 aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
1115 // issue 20222. Try to unite only edges shared by two same faces
1118 SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
1119 while (!edges.empty()) {
1121 sideEdges.splice(sideEdges.end(), edges, edges.begin());
1122 bool sameSide = true;
1123 while (!edges.empty() && sameSide) {
1125 SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
1126 twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
1128 sideEdges.splice(sideEdges.end(), edges, edges.begin());
1130 if (nbSides == 0) { // go backward from the first edge
1132 while (!edges.empty() && sameSide) {
1134 SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
1135 twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
1137 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
1140 list<TopoDS_Edge>::iterator ite = sideEdges.begin();
1141 aNbNodes[nbSides] = 1;
1142 for (; ite!=sideEdges.end(); ite++) {
1143 SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
1144 MapShapeNbElemsItr anIt = aResMap.find(sm);
1145 if (anIt==aResMap.end()) {
1148 std::vector<int> aVec = (*anIt).second;
1150 aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
1152 aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
1160 nbSides = nbEdgesInWire.front();
1161 error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
1169 //=============================================================================
1173 //=============================================================================
1176 StdMeshers_Quadrangle_2D::CheckAnd2Dcompute (SMESH_Mesh & aMesh,
1177 const TopoDS_Shape & aShape,
1178 const bool CreateQuadratic)
1180 _quadraticMesh = CreateQuadratic;
1182 FaceQuadStruct::Ptr quad = CheckNbEdges(aMesh, aShape);
1185 // set normalized grid on unit square in parametric domain
1186 if (!SetNormalizedGrid(aMesh, aShape, quad))
1192 //=============================================================================
1196 //=============================================================================
1198 faceQuadStruct::~faceQuadStruct()
1200 for (size_t i = 0; i < side.size(); i++) {
1203 for (size_t j = i+1; j < side.size(); j++)
1204 if ( side[i] == side[j] )
1218 inline const vector<UVPtStruct>& getUVPtStructIn(FaceQuadStruct::Ptr& quad, int i, int nbSeg)
1220 bool isXConst = (i == QUAD_BOTTOM_SIDE || i == QUAD_TOP_SIDE);
1221 double constValue = (i == QUAD_BOTTOM_SIDE || i == QUAD_LEFT_SIDE) ? 0 : 1;
1223 quad->isEdgeOut[i] ?
1224 quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
1225 quad->side[i]->GetUVPtStruct(isXConst,constValue);
1227 inline gp_UV calcUV(double x, double y,
1228 const gp_UV& a0,const gp_UV& a1,const gp_UV& a2,const gp_UV& a3,
1229 const gp_UV& p0,const gp_UV& p1,const gp_UV& p2,const gp_UV& p3)
1232 ((1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3 ) -
1233 ((1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3);
1237 //=============================================================================
1241 //=============================================================================
1243 bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
1244 const TopoDS_Shape& aShape,
1245 FaceQuadStruct::Ptr & quad)
1247 // Algorithme décrit dans "Génération automatique de maillages"
1248 // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
1249 // traitement dans le domaine paramétrique 2d u,v
1250 // transport - projection sur le carré unité
1252 // MESSAGE("StdMeshers_Quadrangle_2D::SetNormalizedGrid");
1253 // const TopoDS_Face& F = TopoDS::Face(aShape);
1255 // 1 --- find orientation of the 4 edges, by test on extrema
1258 // |<----north-2-------^ a3 -------------> a2
1260 // west-3 east-1 =right | |
1264 // v----south-0--------> a0 -------------> a1
1269 // 3 --- 2D normalized values on unit square [0..1][0..1]
1271 int nbhoriz = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
1272 int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
1274 quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints());
1275 quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints());
1276 quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints());
1277 quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints());
1279 UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
1281 const vector<UVPtStruct>& uv_e0 = getUVPtStructIn(quad, 0, nbhoriz - 1);
1282 const vector<UVPtStruct>& uv_e1 = getUVPtStructIn(quad, 1, nbvertic - 1);
1283 const vector<UVPtStruct>& uv_e2 = getUVPtStructIn(quad, 2, nbhoriz - 1);
1284 const vector<UVPtStruct>& uv_e3 = getUVPtStructIn(quad, 3, nbvertic - 1);
1286 if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
1287 //return error("Can't find nodes on sides");
1288 return error(COMPERR_BAD_INPUT_MESH);
1291 UpdateDegenUV( quad );
1293 // copy data of face boundary
1294 /*if (! quad->isEdgeOut[0])*/ {
1296 for (int i = 0; i < nbhoriz; i++) // down
1297 uv_grid[ j * nbhoriz + i ] = uv_e0[i];
1299 /*if (! quad->isEdgeOut[1])*/ {
1300 const int i = nbhoriz - 1;
1301 for (int j = 0; j < nbvertic; j++) // right
1302 uv_grid[ j * nbhoriz + i ] = uv_e1[j];
1304 /*if (! quad->isEdgeOut[2])*/ {
1305 const int j = nbvertic - 1;
1306 for (int i = 0; i < nbhoriz; i++) // up
1307 uv_grid[ j * nbhoriz + i ] = uv_e2[i];
1309 /*if (! quad->isEdgeOut[3])*/ {
1311 for (int j = 0; j < nbvertic; j++) // left
1312 uv_grid[ j * nbhoriz + i ] = uv_e3[j];
1315 // normalized 2d parameters on grid
1316 for (int i = 0; i < nbhoriz; i++) {
1317 for (int j = 0; j < nbvertic; j++) {
1318 int ij = j * nbhoriz + i;
1319 // --- droite i cste : x = x0 + y(x1-x0)
1320 double x0 = uv_e0[i].normParam; // bas - sud
1321 double x1 = uv_e2[i].normParam; // haut - nord
1322 // --- droite j cste : y = y0 + x(y1-y0)
1323 double y0 = uv_e3[j].normParam; // gauche-ouest
1324 double y1 = uv_e1[j].normParam; // droite - est
1325 // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
1326 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1327 double y = y0 + x * (y1 - y0);
1330 //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
1331 //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
1335 // 4 --- projection on 2d domain (u,v)
1336 gp_UV a0 (uv_e0.front().u, uv_e0.front().v);
1337 gp_UV a1 (uv_e0.back().u, uv_e0.back().v );
1338 gp_UV a2 (uv_e2.back().u, uv_e2.back().v );
1339 gp_UV a3 (uv_e2.front().u, uv_e2.front().v);
1341 for (int i = 0; i < nbhoriz; i++)
1343 gp_UV p0( uv_e0[i].u, uv_e0[i].v );
1344 gp_UV p2( uv_e2[i].u, uv_e2[i].v );
1345 for (int j = 0; j < nbvertic; j++)
1347 gp_UV p1( uv_e1[j].u, uv_e1[j].v );
1348 gp_UV p3( uv_e3[j].u, uv_e3[j].v );
1350 int ij = j * nbhoriz + i;
1351 double x = uv_grid[ij].x;
1352 double y = uv_grid[ij].y;
1354 gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
1356 uv_grid[ij].u = uv.X();
1357 uv_grid[ij].v = uv.Y();
1363 //=======================================================================
1364 //function : ShiftQuad
1365 //purpose : auxilary function for ComputeQuadPref
1366 //=======================================================================
1368 static void shiftQuad(FaceQuadStruct::Ptr& quad, const int num, bool)
1370 quad->shift( num, /*ori=*/true );
1373 //================================================================================
1375 * \brief Rotate sides of a quad by nb
1376 * \param nb - number of rotation quartes
1377 * \param ori - to keep orientation of sides as in an unit quad or not
1379 //================================================================================
1381 void FaceQuadStruct::shift( size_t nb, bool ori )
1383 if ( nb == 0 ) return;
1384 StdMeshers_FaceSide* sideArr[4] = { side[0], side[1], side[2], side[3] };
1385 for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) {
1386 int id = (i + nb) % NB_QUAD_SIDES;
1387 bool wasForward = (i < QUAD_TOP_SIDE);
1388 bool newForward = (id < QUAD_TOP_SIDE);
1389 if (ori && wasForward != newForward)
1390 sideArr[ i ]->Reverse();
1391 side[ id ] = sideArr[ i ];
1395 //=======================================================================
1397 //purpose : auxilary function for ComputeQuadPref
1398 //=======================================================================
1400 static gp_UV calcUV(double x0, double x1, double y0, double y1,
1401 FaceQuadStruct::Ptr& quad,
1402 const gp_UV& a0, const gp_UV& a1,
1403 const gp_UV& a2, const gp_UV& a3)
1405 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
1406 double y = y0 + x * (y1 - y0);
1408 gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE]->Value2d(x).XY();
1409 gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ]->Value2d(y).XY();
1410 gp_UV p2 = quad->side[QUAD_TOP_SIDE ]->Value2d(x).XY();
1411 gp_UV p3 = quad->side[QUAD_LEFT_SIDE ]->Value2d(y).XY();
1413 gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
1418 //=======================================================================
1419 //function : calcUV2
1420 //purpose : auxilary function for ComputeQuadPref
1421 //=======================================================================
1423 static gp_UV calcUV2(double x, double y,
1424 FaceQuadStruct::Ptr& quad,
1425 const gp_UV& a0, const gp_UV& a1,
1426 const gp_UV& a2, const gp_UV& a3)
1428 gp_UV p0 = quad->side[QUAD_BOTTOM_SIDE]->Value2d(x).XY();
1429 gp_UV p1 = quad->side[QUAD_RIGHT_SIDE ]->Value2d(y).XY();
1430 gp_UV p2 = quad->side[QUAD_TOP_SIDE ]->Value2d(x).XY();
1431 gp_UV p3 = quad->side[QUAD_LEFT_SIDE ]->Value2d(y).XY();
1433 gp_UV uv = calcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
1439 //=======================================================================
1441 * Create only quandrangle faces
1443 //=======================================================================
1445 bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh,
1446 const TopoDS_Shape& aShape,
1447 FaceQuadStruct::Ptr quad)
1449 // Auxilary key in order to keep old variant
1450 // of meshing after implementation new variant
1451 // for bug 0016220 from Mantis.
1452 bool OldVersion = false;
1453 if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
1456 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
1457 const TopoDS_Face& F = TopoDS::Face(aShape);
1458 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1460 int i,j,geomFaceID = meshDS->ShapeToIndex(F);
1462 int nb = quad->side[0]->NbPoints();
1463 int nr = quad->side[1]->NbPoints();
1464 int nt = quad->side[2]->NbPoints();
1465 int nl = quad->side[3]->NbPoints();
1466 int dh = abs(nb-nt);
1467 int dv = abs(nr-nl);
1471 // it is a base case => not shift quad but me be replacement is need
1472 shiftQuad(quad,0,WisF);
1475 // we have to shift quad on 2
1476 shiftQuad(quad,2,WisF);
1481 // we have to shift quad on 1
1482 shiftQuad(quad,1,WisF);
1485 // we have to shift quad on 3
1486 shiftQuad(quad,3,WisF);
1490 nb = quad->side[0]->NbPoints();
1491 nr = quad->side[1]->NbPoints();
1492 nt = quad->side[2]->NbPoints();
1493 nl = quad->side[3]->NbPoints();
1496 int nbh = Max(nb,nt);
1497 int nbv = Max(nr,nl);
1501 // ----------- Old version ---------------
1502 // orientation of face and 3 main domain for future faces
1508 // left | | | | rigth
1515 // ----------- New version ---------------
1516 // orientation of face and 3 main domain for future faces
1522 // left |/________\| rigth
1538 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
1539 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
1540 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
1541 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
1543 if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
1544 return error(COMPERR_BAD_INPUT_MESH);
1547 UpdateDegenUV( quad );
1549 // arrays for normalized params
1550 TColStd_SequenceOfReal npb, npr, npt, npl;
1551 for (i=0; i<nb; i++) {
1552 npb.Append(uv_eb[i].normParam);
1554 for (i=0; i<nr; i++) {
1555 npr.Append(uv_er[i].normParam);
1557 for (i=0; i<nt; i++) {
1558 npt.Append(uv_et[i].normParam);
1560 for (i=0; i<nl; i++) {
1561 npl.Append(uv_el[i].normParam);
1566 // add some params to right and left after the first param
1569 double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
1570 for (i=1; i<=dr; i++) {
1571 npr.InsertAfter(1,npr.Value(2)-dpr);
1575 dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
1576 for (i=1; i<=dl; i++) {
1577 npl.InsertAfter(1,npl.Value(2)-dpr);
1581 gp_XY a0(uv_eb.front().u, uv_eb.front().v);
1582 gp_XY a1(uv_eb.back().u, uv_eb.back().v);
1583 gp_XY a2(uv_et.back().u, uv_et.back().v);
1584 gp_XY a3(uv_et.front().u, uv_et.front().v);
1586 int nnn = Min(nr,nl);
1587 // auxilary sequence of XY for creation nodes
1588 // in the bottom part of central domain
1589 // Length of UVL and UVR must be == nbv-nnn
1590 TColgp_SequenceOfXY UVL, UVR, UVT;
1593 // step1: create faces for left domain
1594 StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
1596 for (j=1; j<=nl; j++)
1597 NodesL.SetValue(1,j,uv_el[j-1].node);
1600 for (i=1; i<=dl; i++)
1601 NodesL.SetValue(i+1,nl,uv_et[i].node);
1602 // create and add needed nodes
1603 TColgp_SequenceOfXY UVtmp;
1604 for (i=1; i<=dl; i++) {
1605 double x0 = npt.Value(i+1);
1608 double y0 = npl.Value(i+1);
1609 double y1 = npr.Value(i+1);
1610 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1611 gp_Pnt P = S->Value(UV.X(),UV.Y());
1612 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1613 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1614 NodesL.SetValue(i+1,1,N);
1615 if (UVL.Length()<nbv-nnn) UVL.Append(UV);
1617 for (j=2; j<nl; j++) {
1618 double y0 = npl.Value(dl+j);
1619 double y1 = npr.Value(dl+j);
1620 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1621 gp_Pnt P = S->Value(UV.X(),UV.Y());
1622 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1623 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1624 NodesL.SetValue(i+1,j,N);
1625 if (i==dl) UVtmp.Append(UV);
1628 for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn; i++) {
1629 UVL.Append(UVtmp.Value(i));
1632 for (i=1; i<=dl; i++) {
1633 for (j=1; j<nl; j++) {
1636 myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
1637 NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
1638 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1642 myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
1643 NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
1644 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1650 // fill UVL using c2d
1651 for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn; i++) {
1652 UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
1656 // step2: create faces for right domain
1657 StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
1659 for (j=1; j<=nr; j++)
1660 NodesR.SetValue(1,j,uv_er[nr-j].node);
1663 for (i=1; i<=dr; i++)
1664 NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
1665 // create and add needed nodes
1666 TColgp_SequenceOfXY UVtmp;
1667 for (i=1; i<=dr; i++) {
1668 double x0 = npt.Value(nt-i);
1671 double y0 = npl.Value(i+1);
1672 double y1 = npr.Value(i+1);
1673 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1674 gp_Pnt P = S->Value(UV.X(),UV.Y());
1675 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1676 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1677 NodesR.SetValue(i+1,nr,N);
1678 if (UVR.Length()<nbv-nnn) UVR.Append(UV);
1680 for (j=2; j<nr; j++) {
1681 double y0 = npl.Value(nbv-j+1);
1682 double y1 = npr.Value(nbv-j+1);
1683 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1684 gp_Pnt P = S->Value(UV.X(),UV.Y());
1685 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1686 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1687 NodesR.SetValue(i+1,j,N);
1688 if (i==dr) UVtmp.Prepend(UV);
1691 for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn; i++) {
1692 UVR.Append(UVtmp.Value(i));
1695 for (i=1; i<=dr; i++) {
1696 for (j=1; j<nr; j++) {
1699 myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
1700 NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
1701 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1705 myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
1706 NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
1707 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1713 // fill UVR using c2d
1714 for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn; i++) {
1715 UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
1719 // step3: create faces for central domain
1720 StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
1721 // add first line using NodesL
1722 for (i=1; i<=dl+1; i++)
1723 NodesC.SetValue(1,i,NodesL(i,1));
1724 for (i=2; i<=nl; i++)
1725 NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
1726 // add last line using NodesR
1727 for (i=1; i<=dr+1; i++)
1728 NodesC.SetValue(nb,i,NodesR(i,nr));
1729 for (i=1; i<nr; i++)
1730 NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
1731 // add top nodes (last columns)
1732 for (i=dl+2; i<nbh-dr; i++)
1733 NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
1734 // add bottom nodes (first columns)
1735 for (i=2; i<nb; i++)
1736 NodesC.SetValue(i,1,uv_eb[i-1].node);
1738 // create and add needed nodes
1739 // add linear layers
1740 for (i=2; i<nb; i++) {
1741 double x0 = npt.Value(dl+i);
1743 for (j=1; j<nnn; j++) {
1744 double y0 = npl.Value(nbv-nnn+j);
1745 double y1 = npr.Value(nbv-nnn+j);
1746 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
1747 gp_Pnt P = S->Value(UV.X(),UV.Y());
1748 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1749 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1750 NodesC.SetValue(i,nbv-nnn+j,N);
1755 // add diagonal layers
1756 gp_UV A2 = UVR.Value(nbv-nnn);
1757 gp_UV A3 = UVL.Value(nbv-nnn);
1758 for (i=1; i<nbv-nnn; i++) {
1759 gp_UV p1 = UVR.Value(i);
1760 gp_UV p3 = UVL.Value(i);
1761 double y = i / double(nbv-nnn);
1762 for (j=2; j<nb; j++) {
1763 double x = npb.Value(j);
1764 gp_UV p0( uv_eb[j-1].u, uv_eb[j-1].v );
1765 gp_UV p2 = UVT.Value( j-1 );
1766 gp_UV UV = calcUV(x, y, a0, a1, A2, A3, p0,p1,p2,p3 );
1767 gp_Pnt P = S->Value(UV.X(),UV.Y());
1768 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1769 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
1770 NodesC.SetValue(j,i+1,N);
1774 for (i=1; i<nb; i++) {
1775 for (j=1; j<nbv; j++) {
1778 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1779 NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1780 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1784 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1785 NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1786 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1792 else { // New version (!OldVersion)
1793 // step1: create faces for bottom rectangle domain
1794 StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
1795 // fill UVL and UVR using c2d
1796 for (j=0; j<nb; j++) {
1797 NodesBRD.SetValue(j+1,1,uv_eb[j].node);
1799 for (i=1; i<nnn-1; i++) {
1800 NodesBRD.SetValue(1,i+1,uv_el[i].node);
1801 NodesBRD.SetValue(nb,i+1,uv_er[i].node);
1802 for (j=2; j<nb; j++) {
1803 double x = npb.Value(j);
1804 double y = (1-x) * npl.Value(i+1) + x * npr.Value(i+1);
1805 gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
1806 gp_Pnt P = S->Value(UV.X(),UV.Y());
1807 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1808 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
1809 NodesBRD.SetValue(j,i+1,N);
1812 for (j=1; j<nnn-1; j++) {
1813 for (i=1; i<nb; i++) {
1816 myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
1817 NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
1818 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1822 myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
1823 NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
1824 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1828 int drl = abs(nr-nl);
1829 // create faces for region C
1830 StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
1831 // add nodes from previous region
1832 for (j=1; j<=nb; j++) {
1833 NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
1835 if ((drl+addv) > 0) {
1840 TColgp_SequenceOfXY UVtmp;
1841 double drparam = npr.Value(nr) - npr.Value(nnn-1);
1842 double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
1844 for (i=1; i<=drl; i++) {
1845 // add existed nodes from right edge
1846 NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
1847 //double dtparam = npt.Value(i+1);
1848 y1 = npr.Value(nnn+i-1); // param on right edge
1849 double dpar = (y1 - npr.Value(nnn-1))/drparam;
1850 y0 = npl.Value(nnn-1) + dpar*dlparam; // param on left edge
1851 double dy = y1 - y0;
1852 for (j=1; j<nb; j++) {
1853 double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
1854 double y = y0 + dy*x;
1855 gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
1856 gp_Pnt P = S->Value(UV.X(),UV.Y());
1857 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1858 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1859 NodesC.SetValue(j,i+1,N);
1862 double dy0 = (1-y0)/(addv+1);
1863 double dy1 = (1-y1)/(addv+1);
1864 for (i=1; i<=addv; i++) {
1865 double yy0 = y0 + dy0*i;
1866 double yy1 = y1 + dy1*i;
1867 double dyy = yy1 - yy0;
1868 for (j=1; j<=nb; j++) {
1869 double x = npt.Value(i+1+drl) +
1870 npb.Value(j) * (npt.Value(nt-i) - npt.Value(i+1+drl));
1871 double y = yy0 + dyy*x;
1872 gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
1873 gp_Pnt P = S->Value(UV.X(),UV.Y());
1874 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1875 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1876 NodesC.SetValue(j,i+drl+1,N);
1883 TColgp_SequenceOfXY UVtmp;
1884 double dlparam = npl.Value(nl) - npl.Value(nnn-1);
1885 double drparam = npr.Value(nnn) - npr.Value(nnn-1);
1886 double y0 = npl.Value(nnn-1);
1887 double y1 = npr.Value(nnn-1);
1888 for (i=1; i<=drl; i++) {
1889 // add existed nodes from right edge
1890 NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
1891 y0 = npl.Value(nnn+i-1); // param on left edge
1892 double dpar = (y0 - npl.Value(nnn-1))/dlparam;
1893 y1 = npr.Value(nnn-1) + dpar*drparam; // param on right edge
1894 double dy = y1 - y0;
1895 for (j=2; j<=nb; j++) {
1896 double x = npb.Value(j)*npt.Value(nt-i);
1897 double y = y0 + dy*x;
1898 gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
1899 gp_Pnt P = S->Value(UV.X(),UV.Y());
1900 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1901 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1902 NodesC.SetValue(j,i+1,N);
1905 double dy0 = (1-y0)/(addv+1);
1906 double dy1 = (1-y1)/(addv+1);
1907 for (i=1; i<=addv; i++) {
1908 double yy0 = y0 + dy0*i;
1909 double yy1 = y1 + dy1*i;
1910 double dyy = yy1 - yy0;
1911 for (j=1; j<=nb; j++) {
1912 double x = npt.Value(i+1) +
1913 npb.Value(j) * (npt.Value(nt-i-drl) - npt.Value(i+1));
1914 double y = yy0 + dyy*x;
1915 gp_UV UV = calcUV2(x, y, quad, a0, a1, a2, a3);
1916 gp_Pnt P = S->Value(UV.X(),UV.Y());
1917 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
1918 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
1919 NodesC.SetValue(j,i+drl+1,N);
1924 for (j=1; j<=drl+addv; j++) {
1925 for (i=1; i<nb; i++) {
1928 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
1929 NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
1930 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1934 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
1935 NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
1936 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1941 StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
1942 for (i=1; i<=nt; i++) {
1943 NodesLast.SetValue(i,2,uv_et[i-1].node);
1946 for (i=n1; i<drl+addv+1; i++) {
1948 NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
1950 for (i=1; i<=nb; i++) {
1952 NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
1954 for (i=drl+addv; i>=n2; i--) {
1956 NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
1958 for (i=1; i<nt; i++) {
1961 myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
1962 NodesLast.Value(i+1,2), NodesLast.Value(i,2));
1963 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1967 myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
1968 NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
1969 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
1972 } // if ((drl+addv) > 0)
1974 } // end new version implementation
1981 //=======================================================================
1983 * Evaluate only quandrangle faces
1985 //=======================================================================
1987 bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh,
1988 const TopoDS_Shape& aShape,
1989 std::vector<int>& aNbNodes,
1990 MapShapeNbElems& aResMap,
1993 // Auxilary key in order to keep old variant
1994 // of meshing after implementation new variant
1995 // for bug 0016220 from Mantis.
1996 bool OldVersion = false;
1997 if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
2000 const TopoDS_Face& F = TopoDS::Face(aShape);
2001 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
2003 int nb = aNbNodes[0];
2004 int nr = aNbNodes[1];
2005 int nt = aNbNodes[2];
2006 int nl = aNbNodes[3];
2007 int dh = abs(nb-nt);
2008 int dv = abs(nr-nl);
2012 // it is a base case => not shift
2015 // we have to shift on 2
2024 // we have to shift quad on 1
2031 // we have to shift quad on 3
2041 int nbh = Max(nb,nt);
2042 int nbv = Max(nr,nl);
2057 // add some params to right and left after the first param
2064 int nnn = Min(nr,nl);
2069 // step1: create faces for left domain
2071 nbNodes += dl*(nl-1);
2072 nbFaces += dl*(nl-1);
2074 // step2: create faces for right domain
2076 nbNodes += dr*(nr-1);
2077 nbFaces += dr*(nr-1);
2079 // step3: create faces for central domain
2080 nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2);
2081 nbFaces += (nb-1)*(nbv-1);
2083 else { // New version (!OldVersion)
2084 nbNodes += (nnn-2)*(nb-2);
2085 nbFaces += (nnn-2)*(nb-1);
2086 int drl = abs(nr-nl);
2087 nbNodes += drl*(nb-1) + addv*nb;
2088 nbFaces += (drl+addv)*(nb-1) + (nt-1);
2089 } // end new version implementation
2091 std::vector<int> aVec(SMDSEntity_Last);
2092 for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
2094 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces;
2095 aVec[SMDSEntity_Node] = nbNodes + nbFaces*4;
2096 if (aNbNodes.size()==5) {
2097 aVec[SMDSEntity_Quad_Triangle] = aNbNodes[3] - 1;
2098 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces - aNbNodes[3] + 1;
2102 aVec[SMDSEntity_Node] = nbNodes;
2103 aVec[SMDSEntity_Quadrangle] = nbFaces;
2104 if (aNbNodes.size()==5) {
2105 aVec[SMDSEntity_Triangle] = aNbNodes[3] - 1;
2106 aVec[SMDSEntity_Quadrangle] = nbFaces - aNbNodes[3] + 1;
2109 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
2110 aResMap.insert(std::make_pair(sm,aVec));
2116 //=============================================================================
2117 /*! Split quadrangle in to 2 triangles by smallest diagonal
2120 //=============================================================================
2121 void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
2123 const SMDS_MeshNode* theNode1,
2124 const SMDS_MeshNode* theNode2,
2125 const SMDS_MeshNode* theNode3,
2126 const SMDS_MeshNode* theNode4)
2128 gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
2129 gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
2130 gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
2131 gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
2132 SMDS_MeshFace* face;
2133 if (a.Distance(c) > b.Distance(d)){
2134 face = myHelper->AddFace(theNode2, theNode4 , theNode1);
2135 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2136 face = myHelper->AddFace(theNode2, theNode3, theNode4);
2137 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2141 face = myHelper->AddFace(theNode1, theNode2 ,theNode3);
2142 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2143 face = myHelper->AddFace(theNode1, theNode3, theNode4);
2144 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
2150 enum uvPos { UV_A0, UV_A1, UV_A2, UV_A3, UV_B, UV_R, UV_T, UV_L, UV_SIZE };
2152 inline SMDS_MeshNode* makeNode( UVPtStruct & uvPt,
2154 FaceQuadStruct::Ptr& quad,
2156 SMESH_MesherHelper* helper,
2157 Handle(Geom_Surface) S)
2159 const vector<UVPtStruct>& uv_eb = quad->side[QUAD_BOTTOM_SIDE]->GetUVPtStruct();
2160 const vector<UVPtStruct>& uv_et = quad->side[QUAD_TOP_SIDE ]->GetUVPtStruct();
2161 double rBot = ( uv_eb.size() - 1 ) * uvPt.normParam;
2162 double rTop = ( uv_et.size() - 1 ) * uvPt.normParam;
2163 int iBot = int( rBot );
2164 int iTop = int( rTop );
2165 double xBot = uv_eb[ iBot ].normParam + ( rBot - iBot ) * ( uv_eb[ iBot+1 ].normParam - uv_eb[ iBot ].normParam );
2166 double xTop = uv_et[ iTop ].normParam + ( rTop - iTop ) * ( uv_et[ iTop+1 ].normParam - uv_et[ iTop ].normParam );
2167 double x = xBot + y * ( xTop - xBot );
2169 gp_UV uv = calcUV(/*x,y=*/x, y,
2170 /*a0,...=*/UVs[UV_A0], UVs[UV_A1], UVs[UV_A2], UVs[UV_A3],
2171 /*p0=*/quad->side[QUAD_BOTTOM_SIDE]->Value2d( x ).XY(),
2173 /*p2=*/quad->side[QUAD_TOP_SIDE ]->Value2d( x ).XY(),
2174 /*p3=*/UVs[ UV_L ]);
2175 gp_Pnt P = S->Value( uv.X(), uv.Y() );
2178 return helper->AddNode(P.X(), P.Y(), P.Z(), 0, uv.X(), uv.Y() );
2181 void reduce42( const vector<UVPtStruct>& curr_base,
2182 vector<UVPtStruct>& next_base,
2184 int & next_base_len,
2185 FaceQuadStruct::Ptr& quad,
2188 SMESH_MesherHelper* helper,
2189 Handle(Geom_Surface)& S)
2191 // add one "HH": nodes a,b,c,d,e and faces 1,2,3,4,5,6
2193 // .-----a-----b i + 1
2204 const SMDS_MeshNode*& Na = next_base[ ++next_base_len ].node;
2206 Na = makeNode( next_base[ next_base_len ], y, quad, UVs, helper, S );
2209 const SMDS_MeshNode*& Nb = next_base[ ++next_base_len ].node;
2211 Nb = makeNode( next_base[ next_base_len ], y, quad, UVs, helper, S );
2214 double u = (curr_base[j + 2].u + next_base[next_base_len - 2].u) / 2.0;
2215 double v = (curr_base[j + 2].v + next_base[next_base_len - 2].v) / 2.0;
2216 gp_Pnt P = S->Value(u,v);
2217 SMDS_MeshNode* Nc = helper->AddNode(P.X(), P.Y(), P.Z(), 0, u, v);
2220 u = (curr_base[j + 2].u + next_base[next_base_len - 1].u) / 2.0;
2221 v = (curr_base[j + 2].v + next_base[next_base_len - 1].v) / 2.0;
2223 SMDS_MeshNode* Nd = helper->AddNode(P.X(), P.Y(), P.Z(), 0, u, v);
2226 u = (curr_base[j + 2].u + next_base[next_base_len].u) / 2.0;
2227 v = (curr_base[j + 2].v + next_base[next_base_len].v) / 2.0;
2229 SMDS_MeshNode* Ne = helper->AddNode(P.X(), P.Y(), P.Z(), 0, u, v);
2232 helper->AddFace(curr_base[j + 0].node,
2233 curr_base[j + 1].node, Nc,
2234 next_base[next_base_len - 2].node);
2236 helper->AddFace(curr_base[j + 1].node,
2237 curr_base[j + 2].node, Nd, Nc);
2239 helper->AddFace(curr_base[j + 2].node,
2240 curr_base[j + 3].node, Ne, Nd);
2242 helper->AddFace(curr_base[j + 3].node,
2243 curr_base[j + 4].node, Nb, Ne);
2245 helper->AddFace(Nc, Nd, Na, next_base[next_base_len - 2].node);
2247 helper->AddFace(Nd, Ne, Nb, Na);
2250 void reduce31( const vector<UVPtStruct>& curr_base,
2251 vector<UVPtStruct>& next_base,
2253 int & next_base_len,
2254 FaceQuadStruct::Ptr& quad,
2257 SMESH_MesherHelper* helper,
2258 Handle(Geom_Surface)& S)
2260 // add one "H": nodes b,c,e and faces 1,2,4,5
2262 // .---------b i + 1
2273 const SMDS_MeshNode*& Nb = next_base[ ++next_base_len ].node;
2275 Nb = makeNode( next_base[ next_base_len ], y, quad, UVs, helper, S );
2278 double u1 = (curr_base[ j ].u + next_base[ next_base_len-1 ].u ) / 2.0;
2279 double u2 = (curr_base[ j+3 ].u + next_base[ next_base_len ].u ) / 2.0;
2280 double u3 = (u2 - u1) / 3.0;
2282 double v1 = (curr_base[ j ].v + next_base[ next_base_len-1 ].v ) / 2.0;
2283 double v2 = (curr_base[ j+3 ].v + next_base[ next_base_len ].v ) / 2.0;
2284 double v3 = (v2 - v1) / 3.0;
2288 gp_Pnt P = S->Value(u,v);
2289 SMDS_MeshNode* Nc = helper->AddNode( P.X(), P.Y(), P.Z(), 0, u, v );
2294 SMDS_MeshNode* Ne = helper->AddNode( P.X(), P.Y(), P.Z(), 0, u, v );
2298 helper->AddFace( curr_base[ j + 0 ].node,
2299 curr_base[ j + 1 ].node,
2301 next_base[ next_base_len - 1 ].node);
2303 helper->AddFace( curr_base[ j + 1 ].node,
2304 curr_base[ j + 2 ].node, Ne, Nc);
2306 helper->AddFace( curr_base[ j + 2 ].node,
2307 curr_base[ j + 3 ].node, Nb, Ne);
2309 helper->AddFace(Nc, Ne, Nb,
2310 next_base[ next_base_len - 1 ].node);
2313 typedef void (* PReduceFunction) ( const vector<UVPtStruct>& curr_base,
2314 vector<UVPtStruct>& next_base,
2316 int & next_base_len,
2317 FaceQuadStruct::Ptr & quad,
2320 SMESH_MesherHelper* helper,
2321 Handle(Geom_Surface)& S);
2325 //=======================================================================
2327 * Implementation of Reduced algorithm (meshing with quadrangles only)
2329 //=======================================================================
2331 bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh,
2332 const TopoDS_Shape& aShape,
2333 FaceQuadStruct::Ptr quad)
2335 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
2336 const TopoDS_Face& F = TopoDS::Face(aShape);
2337 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
2338 int i,j,geomFaceID = meshDS->ShapeToIndex(F);
2340 int nb = quad->side[0]->NbPoints(); // bottom
2341 int nr = quad->side[1]->NbPoints(); // right
2342 int nt = quad->side[2]->NbPoints(); // top
2343 int nl = quad->side[3]->NbPoints(); // left
2345 // Simple Reduce 10->8->6->4 (3 steps) Multiple Reduce 10->4 (1 step)
2347 // .-----.-----.-----.-----. .-----.-----.-----.-----.
2348 // | / \ | / \ | | / \ | / \ |
2349 // | / .--.--. \ | | / \ | / \ |
2350 // | / / | \ \ | | / .----.----. \ |
2351 // .---.---.---.---.---.---. | / / \ | / \ \ |
2352 // | / / \ | / \ \ | | / / \ | / \ \ |
2353 // | / / .-.-. \ \ | | / / .---.---. \ \ |
2354 // | / / / | \ \ \ | | / / / \ | / \ \ \ |
2355 // .--.--.--.--.--.--.--.--. | / / / \ | / \ \ \ |
2356 // | / / / \ | / \ \ \ | | / / / .-.-. \ \ \ |
2357 // | / / / .-.-. \ \ \ | | / / / / | \ \ \ \ |
2358 // | / / / / | \ \ \ \ | | / / / / | \ \ \ \ |
2359 // .-.-.-.--.--.--.--.-.-.-. .-.-.-.--.--.--.--.-.-.-.
2361 bool MultipleReduce = false;
2373 else if (nb == nt) {
2374 nr1 = nb; // and == nt
2388 // number of rows and columns
2389 int nrows = nr1 - 1;
2390 int ncol_top = nt1 - 1;
2391 int ncol_bot = nb1 - 1;
2392 // number of rows needed to reduce ncol_bot to ncol_top using simple 3->1 "tree" (see below)
2394 int( ceil( log( double(ncol_bot) / ncol_top) / log( 3.))); // = log x base 3
2395 if ( nrows < nrows_tree31 )
2397 MultipleReduce = true;
2398 error( COMPERR_WARNING,
2399 SMESH_Comment("To use 'Reduced' transition, "
2400 "number of face rows should be at least ")
2401 << nrows_tree31 << ". Actual number of face rows is " << nrows << ". "
2402 "'Quadrangle preference (reversed)' transion has been used.");
2406 if (MultipleReduce) { // == ComputeQuadPref QUAD_QUADRANGLE_PREF_REVERSED
2407 //==================================================
2408 int dh = abs(nb-nt);
2409 int dv = abs(nr-nl);
2413 // it is a base case => not shift quad but may be replacement is need
2414 shiftQuad(quad,0,true);
2417 // we have to shift quad on 2
2418 shiftQuad(quad,2,true);
2423 // we have to shift quad on 1
2424 shiftQuad(quad,1,true);
2427 // we have to shift quad on 3
2428 shiftQuad(quad,3,true);
2432 nb = quad->side[0]->NbPoints();
2433 nr = quad->side[1]->NbPoints();
2434 nt = quad->side[2]->NbPoints();
2435 nl = quad->side[3]->NbPoints();
2438 int nbh = Max(nb,nt);
2439 int nbv = Max(nr,nl);
2452 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
2453 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
2454 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
2455 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
2457 if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
2458 return error(COMPERR_BAD_INPUT_MESH);
2461 UpdateDegenUV( quad );
2463 // arrays for normalized params
2464 TColStd_SequenceOfReal npb, npr, npt, npl;
2465 for (j = 0; j < nb; j++) {
2466 npb.Append(uv_eb[j].normParam);
2468 for (i = 0; i < nr; i++) {
2469 npr.Append(uv_er[i].normParam);
2471 for (j = 0; j < nt; j++) {
2472 npt.Append(uv_et[j].normParam);
2474 for (i = 0; i < nl; i++) {
2475 npl.Append(uv_el[i].normParam);
2479 // orientation of face and 3 main domain for future faces
2485 // left | | | | rigth
2492 // add some params to right and left after the first param
2495 double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
2496 for (i=1; i<=dr; i++) {
2497 npr.InsertAfter(1,npr.Value(2)-dpr);
2501 dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
2502 for (i=1; i<=dl; i++) {
2503 npl.InsertAfter(1,npl.Value(2)-dpr);
2506 gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
2507 gp_XY a1 (uv_eb.back().u, uv_eb.back().v);
2508 gp_XY a2 (uv_et.back().u, uv_et.back().v);
2509 gp_XY a3 (uv_et.front().u, uv_et.front().v);
2511 int nnn = Min(nr,nl);
2512 // auxilary sequence of XY for creation of nodes
2513 // in the bottom part of central domain
2514 // it's length must be == nbv-nnn-1
2515 TColgp_SequenceOfXY UVL;
2516 TColgp_SequenceOfXY UVR;
2517 //==================================================
2519 // step1: create faces for left domain
2520 StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
2522 for (j=1; j<=nl; j++)
2523 NodesL.SetValue(1,j,uv_el[j-1].node);
2526 for (i=1; i<=dl; i++)
2527 NodesL.SetValue(i+1,nl,uv_et[i].node);
2528 // create and add needed nodes
2529 TColgp_SequenceOfXY UVtmp;
2530 for (i=1; i<=dl; i++) {
2531 double x0 = npt.Value(i+1);
2534 double y0 = npl.Value(i+1);
2535 double y1 = npr.Value(i+1);
2536 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2537 gp_Pnt P = S->Value(UV.X(),UV.Y());
2538 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2539 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2540 NodesL.SetValue(i+1,1,N);
2541 if (UVL.Length()<nbv-nnn-1) UVL.Append(UV);
2543 for (j=2; j<nl; j++) {
2544 double y0 = npl.Value(dl+j);
2545 double y1 = npr.Value(dl+j);
2546 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2547 gp_Pnt P = S->Value(UV.X(),UV.Y());
2548 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2549 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2550 NodesL.SetValue(i+1,j,N);
2551 if (i==dl) UVtmp.Append(UV);
2554 for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
2555 UVL.Append(UVtmp.Value(i));
2558 for (i=1; i<=dl; i++) {
2559 for (j=1; j<nl; j++) {
2561 myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
2562 NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
2563 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2568 // fill UVL using c2d
2569 for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
2570 UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
2574 // step2: create faces for right domain
2575 StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
2577 for (j=1; j<=nr; j++)
2578 NodesR.SetValue(1,j,uv_er[nr-j].node);
2581 for (i=1; i<=dr; i++)
2582 NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
2583 // create and add needed nodes
2584 TColgp_SequenceOfXY UVtmp;
2585 for (i=1; i<=dr; i++) {
2586 double x0 = npt.Value(nt-i);
2589 double y0 = npl.Value(i+1);
2590 double y1 = npr.Value(i+1);
2591 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2592 gp_Pnt P = S->Value(UV.X(),UV.Y());
2593 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2594 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2595 NodesR.SetValue(i+1,nr,N);
2596 if (UVR.Length()<nbv-nnn-1) UVR.Append(UV);
2598 for (j=2; j<nr; j++) {
2599 double y0 = npl.Value(nbv-j+1);
2600 double y1 = npr.Value(nbv-j+1);
2601 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2602 gp_Pnt P = S->Value(UV.X(),UV.Y());
2603 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2604 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2605 NodesR.SetValue(i+1,j,N);
2606 if (i==dr) UVtmp.Prepend(UV);
2609 for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
2610 UVR.Append(UVtmp.Value(i));
2613 for (i=1; i<=dr; i++) {
2614 for (j=1; j<nr; j++) {
2616 myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
2617 NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
2618 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2623 // fill UVR using c2d
2624 for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
2625 UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
2629 // step3: create faces for central domain
2630 StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
2631 // add first line using NodesL
2632 for (i=1; i<=dl+1; i++)
2633 NodesC.SetValue(1,i,NodesL(i,1));
2634 for (i=2; i<=nl; i++)
2635 NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
2636 // add last line using NodesR
2637 for (i=1; i<=dr+1; i++)
2638 NodesC.SetValue(nb,i,NodesR(i,nr));
2639 for (i=1; i<nr; i++)
2640 NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
2641 // add top nodes (last columns)
2642 for (i=dl+2; i<nbh-dr; i++)
2643 NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
2644 // add bottom nodes (first columns)
2645 for (i=2; i<nb; i++)
2646 NodesC.SetValue(i,1,uv_eb[i-1].node);
2648 // create and add needed nodes
2649 // add linear layers
2650 for (i=2; i<nb; i++) {
2651 double x0 = npt.Value(dl+i);
2653 for (j=1; j<nnn; j++) {
2654 double y0 = npl.Value(nbv-nnn+j);
2655 double y1 = npr.Value(nbv-nnn+j);
2656 gp_UV UV = calcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
2657 gp_Pnt P = S->Value(UV.X(),UV.Y());
2658 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2659 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
2660 NodesC.SetValue(i,nbv-nnn+j,N);
2663 // add diagonal layers
2664 for (i=1; i<nbv-nnn; i++) {
2665 double du = UVR.Value(i).X() - UVL.Value(i).X();
2666 double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
2667 for (j=2; j<nb; j++) {
2668 double u = UVL.Value(i).X() + du*npb.Value(j);
2669 double v = UVL.Value(i).Y() + dv*npb.Value(j);
2670 gp_Pnt P = S->Value(u,v);
2671 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
2672 meshDS->SetNodeOnFace(N, geomFaceID, u, v);
2673 NodesC.SetValue(j,i+1,N);
2677 for (i=1; i<nb; i++) {
2678 for (j=1; j<nbv; j++) {
2680 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
2681 NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
2682 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
2685 } // end Multiple Reduce implementation
2686 else { // Simple Reduce (!MultipleReduce)
2687 //=========================================================
2690 // it is a base case => not shift quad
2691 //shiftQuad(quad,0,true);
2694 // we have to shift quad on 2
2695 shiftQuad(quad,2,true);
2700 // we have to shift quad on 1
2701 shiftQuad(quad,1,true);
2704 // we have to shift quad on 3
2705 shiftQuad(quad,3,true);
2709 nb = quad->side[0]->NbPoints();
2710 nr = quad->side[1]->NbPoints();
2711 nt = quad->side[2]->NbPoints();
2712 nl = quad->side[3]->NbPoints();
2714 // number of rows and columns
2715 int nrows = nr - 1; // and also == nl - 1
2716 int ncol_top = nt - 1;
2717 int ncol_bot = nb - 1;
2718 int npair_top = ncol_top / 2;
2719 // maximum number of bottom elements for "linear" simple reduce 4->2
2720 int max_lin42 = ncol_top + npair_top * 2 * nrows;
2721 // maximum number of bottom elements for "linear" simple reduce 3->1
2722 int max_lin31 = ncol_top + ncol_top * 2 * nrows;
2723 // maximum number of bottom elements for "tree" simple reduce 4->2
2725 // number of rows needed to reduce ncol_bot to ncol_top using simple 4->2 "tree"
2726 int nrows_tree42 = int( log( (double)(ncol_bot / ncol_top) )/log((double)2) ); // needed to avoid overflow at pow(2) while computing max_tree42
2727 if (nrows_tree42 < nrows) {
2728 max_tree42 = npair_top * pow(2.0, nrows + 1);
2729 if ( ncol_top > npair_top * 2 ) {
2730 int delta = ncol_bot - max_tree42;
2731 for (int irow = 1; irow < nrows; irow++) {
2732 int nfour = delta / 4;
2735 if (delta <= (ncol_top - npair_top * 2))
2736 max_tree42 = ncol_bot;
2739 // maximum number of bottom elements for "tree" simple reduce 3->1
2740 //int max_tree31 = ncol_top * pow(3.0, nrows);
2741 bool is_lin_31 = false;
2742 bool is_lin_42 = false;
2743 bool is_tree_31 = false;
2744 bool is_tree_42 = false;
2745 int max_lin = max_lin42;
2746 if (ncol_bot > max_lin42) {
2747 if (ncol_bot <= max_lin31) {
2749 max_lin = max_lin31;
2753 // if ncol_bot is a 3*n or not 2*n
2754 if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
2756 max_lin = max_lin31;
2762 if (ncol_bot > max_lin) { // not "linear"
2763 is_tree_31 = (ncol_bot > max_tree42);
2764 if (ncol_bot <= max_tree42) {
2765 if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
2774 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
2775 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
2776 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
2777 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
2779 if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
2780 return error(COMPERR_BAD_INPUT_MESH);
2782 myHelper->SetElementsOnShape( true );
2784 gp_UV uv[ UV_SIZE ];
2785 uv[ UV_A0 ].SetCoord( uv_eb.front().u, uv_eb.front().v);
2786 uv[ UV_A1 ].SetCoord( uv_eb.back().u, uv_eb.back().v );
2787 uv[ UV_A2 ].SetCoord( uv_et.back().u, uv_et.back().v );
2788 uv[ UV_A3 ].SetCoord( uv_et.front().u, uv_et.front().v);
2790 vector<UVPtStruct> curr_base = uv_eb, next_base;
2792 UVPtStruct nullUVPtStruct; nullUVPtStruct.node = 0;
2794 int curr_base_len = nb;
2795 int next_base_len = 0;
2798 { // ------------------------------------------------------------------
2799 // New algorithm implemented by request of IPAL22856
2800 // "2D quadrangle mesher of reduced type works wrong"
2801 // http://bugtracker.opencascade.com/show_bug.cgi?id=22856
2803 // the algorithm is following: all reduces are centred in horizontal
2804 // direction and are distributed among all rows
2806 if (ncol_bot > max_tree42) {
2810 if ((ncol_top/3)*3 == ncol_top ) {
2818 const int col_top_size = is_lin_42 ? 2 : 1;
2819 const int col_base_size = is_lin_42 ? 4 : 3;
2821 // Compute nb of "columns" (like in "linear" simple reducing) in all rows
2823 vector<int> nb_col_by_row;
2825 int delta_all = nb - nt;
2826 int delta_one_col = nrows * 2;
2827 int nb_col = delta_all / delta_one_col;
2828 int remainder = delta_all - nb_col * delta_one_col;
2829 if (remainder > 0) {
2832 if ( nb_col * col_top_size >= nt ) // == "tree" reducing situation
2834 // top row is full (all elements reduced), add "columns" one by one
2835 // in rows below until all bottom elements are reduced
2836 nb_col = ( nt - 1 ) / col_top_size;
2837 nb_col_by_row.resize( nrows, nb_col );
2838 int nbrows_not_full = nrows - 1;
2839 int cur_top_size = nt - 1;
2840 remainder = delta_all - nb_col * delta_one_col;
2841 while ( remainder > 0 )
2843 delta_one_col = nbrows_not_full * 2;
2844 int nb_col_add = remainder / delta_one_col;
2845 cur_top_size += 2 * nb_col_by_row[ nbrows_not_full ];
2846 int nb_col_free = cur_top_size / col_top_size - nb_col_by_row[ nbrows_not_full-1 ];
2847 if ( nb_col_add > nb_col_free )
2848 nb_col_add = nb_col_free;
2849 for ( int irow = 0; irow < nbrows_not_full; ++irow )
2850 nb_col_by_row[ irow ] += nb_col_add;
2852 remainder -= nb_col_add * delta_one_col;
2855 else // == "linear" reducing situation
2857 nb_col_by_row.resize( nrows, nb_col );
2859 for ( int irow = remainder / 2; irow < nrows; ++irow )
2860 nb_col_by_row[ irow ]--;
2865 PReduceFunction reduceFunction = & ( is_lin_42 ? reduce42 : reduce31 );
2867 const int reduce_grp_size = is_lin_42 ? 4 : 3;
2869 for (i = 1; i < nr; i++) // layer by layer
2871 nb_col = nb_col_by_row[ i-1 ];
2872 int nb_next = curr_base_len - nb_col * 2;
2873 if (nb_next < nt) nb_next = nt;
2875 const double y = uv_el[ i ].normParam;
2877 if ( i + 1 == nr ) // top
2884 next_base.resize( nb_next, nullUVPtStruct );
2885 next_base.front() = uv_el[i];
2886 next_base.back() = uv_er[i];
2888 // compute normalized param u
2889 double du = 1. / ( nb_next - 1 );
2890 next_base[0].normParam = 0.;
2891 for ( j = 1; j < nb_next; ++j )
2892 next_base[j].normParam = next_base[j-1].normParam + du;
2894 uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v );
2895 uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v );
2897 int free_left = ( curr_base_len - 1 - nb_col * col_base_size ) / 2;
2898 int free_middle = curr_base_len - 1 - nb_col * col_base_size - 2 * free_left;
2900 // not reduced left elements
2901 for (j = 0; j < free_left; j++)
2904 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
2906 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
2908 myHelper->AddFace(curr_base[ j ].node,
2909 curr_base[ j+1 ].node,
2911 next_base[ next_base_len-1 ].node);
2914 for (int icol = 1; icol <= nb_col; icol++)
2917 reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S );
2919 j += reduce_grp_size;
2921 // elements in the middle of "columns" added for symmetry
2922 if ( free_middle > 0 && ( nb_col % 2 == 0 ) && icol == nb_col / 2 )
2924 for (int imiddle = 1; imiddle <= free_middle; imiddle++) {
2925 // f (i + 1, j + imiddle)
2926 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
2928 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
2930 myHelper->AddFace(curr_base[ j-1+imiddle ].node,
2931 curr_base[ j +imiddle ].node,
2933 next_base[ next_base_len-1 ].node);
2939 // not reduced right elements
2940 for (; j < curr_base_len-1; j++) {
2942 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
2944 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
2946 myHelper->AddFace(curr_base[ j ].node,
2947 curr_base[ j+1 ].node,
2949 next_base[ next_base_len-1 ].node);
2952 curr_base_len = next_base_len + 1;
2954 curr_base.swap( next_base );
2958 else if ( is_tree_42 || is_tree_31 )
2960 // "tree" simple reduce "42": 2->4->8->16->32->...
2962 // .-------------------------------.-------------------------------. nr
2964 // | \ .---------------.---------------. / |
2966 // .---------------.---------------.---------------.---------------.
2967 // | \ | / | \ | / |
2968 // | \ .-------.-------. / | \ .-------.-------. / |
2969 // | | | | | | | | |
2970 // .-------.-------.-------.-------.-------.-------.-------.-------. i
2971 // |\ | /|\ | /|\ | /|\ | /|
2972 // | \.---.---./ | \.---.---./ | \.---.---./ | \.---.---./ |
2973 // | | | | | | | | | | | | | | | | |
2974 // .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.
2975 // |\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|\ | /|
2976 // | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. | .-.-. |
2977 // | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | |
2978 // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
2981 // "tree" simple reduce "31": 1->3->9->27->...
2983 // .-----------------------------------------------------. nr
2985 // | .-----------------. |
2987 // .-----------------.-----------------.-----------------.
2988 // | \ / | \ / | \ / |
2989 // | .-----. | .-----. | .-----. | i
2990 // | | | | | | | | | |
2991 // .-----.-----.-----.-----.-----.-----.-----.-----.-----.
2992 // |\ /|\ /|\ /|\ /|\ /|\ /|\ /|\ /|\ /|
2993 // | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. | .-. |
2994 // | | | | | | | | | | | | | | | | | | | | | | | | | | | |
2995 // .-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. 1
2998 PReduceFunction reduceFunction = & ( is_tree_42 ? reduce42 : reduce31 );
3000 const int reduce_grp_size = is_tree_42 ? 4 : 3;
3002 for (i = 1; i < nr; i++) // layer by layer
3004 // to stop reducing, if number of nodes reaches nt
3005 int delta = curr_base_len - nt;
3007 // to calculate normalized parameter, we must know number of points in next layer
3008 int nb_reduce_groups = (curr_base_len - 1) / reduce_grp_size;
3009 int nb_next = nb_reduce_groups * (reduce_grp_size-2) + (curr_base_len - nb_reduce_groups*reduce_grp_size);
3010 if (nb_next < nt) nb_next = nt;
3012 const double y = uv_el[ i ].normParam;
3014 if ( i + 1 == nr ) // top
3021 next_base.resize( nb_next, nullUVPtStruct );
3022 next_base.front() = uv_el[i];
3023 next_base.back() = uv_er[i];
3025 // compute normalized param u
3026 double du = 1. / ( nb_next - 1 );
3027 next_base[0].normParam = 0.;
3028 for ( j = 1; j < nb_next; ++j )
3029 next_base[j].normParam = next_base[j-1].normParam + du;
3031 uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v );
3032 uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v );
3034 for (j = 0; j+reduce_grp_size < curr_base_len && delta > 0; j+=reduce_grp_size, delta-=2)
3036 reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S );
3039 // not reduced side elements (if any)
3040 for (; j < curr_base_len-1; j++)
3043 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
3045 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
3047 myHelper->AddFace(curr_base[ j ].node,
3048 curr_base[ j+1 ].node,
3050 next_base[ next_base_len-1 ].node);
3052 curr_base_len = next_base_len + 1;
3054 curr_base.swap( next_base );
3056 } // end "tree" simple reduce
3058 else if ( is_lin_42 || is_lin_31 ) {
3059 // "linear" simple reduce "31": 2->6->10->14
3061 // .-----------------------------.-----------------------------. nr
3063 // | .---------. | .---------. |
3065 // .---------.---------.---------.---------.---------.---------.
3066 // | / \ / \ | / \ / \ |
3067 // | / .-----. \ | / .-----. \ | i
3068 // | / | | \ | / | | \ |
3069 // .-----.-----.-----.-----.-----.-----.-----.-----.-----.-----.
3070 // | / / \ / \ \ | / / \ / \ \ |
3071 // | / / .-. \ \ | / / .-. \ \ |
3072 // | / / / \ \ \ | / / / \ \ \ |
3073 // .--.----.---.-----.---.-----.-.--.----.---.-----.---.-----.-. 1
3076 // "linear" simple reduce "42": 4->8->12->16
3078 // .---------------.---------------.---------------.---------------. nr
3079 // | \ | / | \ | / |
3080 // | \ .-------.-------. / | \ .-------.-------. / |
3081 // | | | | | | | | |
3082 // .-------.-------.-------.-------.-------.-------.-------.-------.
3083 // | / \ | / \ | / \ | / \ |
3084 // | / \.----.----./ \ | / \.----.----./ \ | i
3085 // | / | | | \ | / | | | \ |
3086 // .-----.----.----.----.----.-----.-----.----.----.----.----.-----.
3087 // | / / \ | / \ \ | / / \ | / \ \ |
3088 // | / / .-.-. \ \ | / / .-.-. \ \ |
3089 // | / / / | \ \ \ | / / / | \ \ \ |
3090 // .---.---.---.---.---.---.---.---.---.---.---.---.---.---.---.---. 1
3093 // nt = 5, nb = 7, nr = 4
3094 //int delta_all = 2;
3095 //int delta_one_col = 6;
3097 //int remainder = 2;
3098 //if (remainder > 0) nb_col++;
3100 //int free_left = 1;
3102 //int free_middle = 4;
3104 int delta_all = nb - nt;
3105 int delta_one_col = (nr - 1) * 2;
3106 int nb_col = delta_all / delta_one_col;
3107 int remainder = delta_all - nb_col * delta_one_col;
3108 if (remainder > 0) {
3111 const int col_top_size = is_lin_42 ? 2 : 1;
3112 int free_left = ((nt - 1) - nb_col * col_top_size) / 2;
3113 free_left += nr - 2;
3114 int free_middle = (nr - 2) * 2;
3115 if (remainder > 0 && nb_col == 1) {
3116 int nb_rows_short_col = remainder / 2;
3117 int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
3118 free_left -= nb_rows_thrown;
3121 // nt = 5, nb = 17, nr = 4
3122 //int delta_all = 12;
3123 //int delta_one_col = 6;
3125 //int remainder = 0;
3126 //int free_left = 2;
3127 //int free_middle = 4;
3129 PReduceFunction reduceFunction = & ( is_lin_42 ? reduce42 : reduce31 );
3131 const int reduce_grp_size = is_lin_42 ? 4 : 3;
3133 for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) // layer by layer
3135 // to calculate normalized parameter, we must know number of points in next layer
3136 int nb_next = curr_base_len - nb_col * 2;
3137 if (remainder > 0 && i > remainder / 2)
3138 // take into account short "column"
3140 if (nb_next < nt) nb_next = nt;
3142 const double y = uv_el[ i ].normParam;
3144 if ( i + 1 == nr ) // top
3151 next_base.resize( nb_next, nullUVPtStruct );
3152 next_base.front() = uv_el[i];
3153 next_base.back() = uv_er[i];
3155 // compute normalized param u
3156 double du = 1. / ( nb_next - 1 );
3157 next_base[0].normParam = 0.;
3158 for ( j = 1; j < nb_next; ++j )
3159 next_base[j].normParam = next_base[j-1].normParam + du;
3161 uv[ UV_L ].SetCoord( next_base.front().u, next_base.front().v );
3162 uv[ UV_R ].SetCoord( next_base.back().u, next_base.back().v );
3164 // not reduced left elements
3165 for (j = 0; j < free_left; j++)
3168 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
3170 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
3172 myHelper->AddFace(curr_base[ j ].node,
3173 curr_base[ j+1 ].node,
3175 next_base[ next_base_len-1 ].node);
3178 for (int icol = 1; icol <= nb_col; icol++) {
3180 if (remainder > 0 && icol == nb_col && i > remainder / 2)
3181 // stop short "column"
3185 reduceFunction( curr_base, next_base, j, next_base_len, quad, uv, y, myHelper, S );
3187 j += reduce_grp_size;
3189 // not reduced middle elements
3190 if (icol < nb_col) {
3191 if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
3192 // pass middle elements before stopped short "column"
3195 int free_add = free_middle;
3196 if (remainder > 0 && icol == nb_col - 1)
3197 // next "column" is short
3198 free_add -= (nr - 1) - (remainder / 2);
3200 for (int imiddle = 1; imiddle <= free_add; imiddle++) {
3201 // f (i + 1, j + imiddle)
3202 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
3204 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
3206 myHelper->AddFace(curr_base[ j-1+imiddle ].node,
3207 curr_base[ j +imiddle ].node,
3209 next_base[ next_base_len-1 ].node);
3215 // not reduced right elements
3216 for (; j < curr_base_len-1; j++) {
3218 const SMDS_MeshNode*& Nf = next_base[++next_base_len].node;
3220 Nf = makeNode( next_base[ next_base_len ], y, quad, uv, myHelper, S );
3222 myHelper->AddFace(curr_base[ j ].node,
3223 curr_base[ j+1 ].node,
3225 next_base[ next_base_len-1 ].node);
3228 curr_base_len = next_base_len + 1;
3230 curr_base.swap( next_base );
3233 } // end "linear" simple reduce
3238 } // end Simple Reduce implementation
3244 //================================================================================
3245 namespace // data for smoothing
3248 // --------------------------------------------------------------------------------
3250 * \brief Structure used to check validity of node position after smoothing.
3251 * It holds two nodes connected to a smoothed node and belonging to
3258 TTriangle( TSmoothNode* n1=0, TSmoothNode* n2=0 ): _n1(n1), _n2(n2) {}
3260 inline bool IsForward( gp_UV uv ) const;
3262 // --------------------------------------------------------------------------------
3264 * \brief Data of a smoothed node
3269 vector< TTriangle > _triangles; // if empty, then node is not movable
3271 // --------------------------------------------------------------------------------
3272 inline bool TTriangle::IsForward( gp_UV uv ) const
3274 gp_Vec2d v1( uv, _n1->_uv ), v2( uv, _n2->_uv );
3280 //================================================================================
3282 * \brief Set UV of nodes on degenerated VERTEXes in the middle of degenerated EDGE
3284 * WARNING: this method must be called AFTER retrieving UVPtStruct's from quad
3286 //================================================================================
3288 void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct::Ptr quad)
3290 for ( unsigned i = 0; i < quad->side.size(); ++i )
3292 StdMeshers_FaceSide* side = quad->side[i];
3293 const vector<UVPtStruct>& uvVec = side->GetUVPtStruct();
3295 // find which end of the side is on degenerated shape
3297 if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() ))
3299 else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() ))
3300 degenInd = uvVec.size() - 1;
3304 // find another side sharing the degenerated shape
3305 bool isPrev = ( degenInd == 0 );
3306 if ( i >= QUAD_TOP_SIDE )
3308 int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4;
3309 StdMeshers_FaceSide* side2 = quad->side[ i2 ];
3310 const vector<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
3312 if ( uvVec[ degenInd ].node == uvVec2[0].node )
3314 else if ( uvVec[ degenInd ].node == uvVec2.back().node )
3315 degenInd2 = uvVec2.size() - 1;
3317 throw SALOME_Exception( LOCALIZED( "Logical error" ));
3319 // move UV in the middle
3320 uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd ]);
3321 uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
3322 uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
3323 uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
3327 //================================================================================
3329 * \brief Perform smoothing of 2D elements on a FACE with ignored degenerated EDGE
3331 //================================================================================
3333 void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct::Ptr quad)
3335 if ( !myNeedSmooth ) return;
3337 // Get nodes to smooth
3339 typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
3340 TNo2SmooNoMap smooNoMap;
3342 const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
3343 SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
3344 SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
3345 SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
3346 while ( nIt->more() ) // loop on nodes bound to a FACE
3348 const SMDS_MeshNode* node = nIt->next();
3349 TSmoothNode & sNode = smooNoMap[ node ];
3350 sNode._uv = myHelper->GetNodeUV( geomFace, node );
3352 // set sNode._triangles
3353 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face );
3354 while ( fIt->more() )
3356 const SMDS_MeshElement* face = fIt->next();
3357 const int nbN = face->NbCornerNodes();
3358 const int nInd = face->GetNodeIndex( node );
3359 const int prevInd = myHelper->WrapIndex( nInd - 1, nbN );
3360 const int nextInd = myHelper->WrapIndex( nInd + 1, nbN );
3361 const SMDS_MeshNode* prevNode = face->GetNode( prevInd );
3362 const SMDS_MeshNode* nextNode = face->GetNode( nextInd );
3363 sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ],
3364 & smooNoMap[ nextNode ]));
3367 // set _uv of smooth nodes on FACE boundary
3368 for ( unsigned i = 0; i < quad->side.size(); ++i )
3370 const vector<UVPtStruct>& uvVec = quad->side[i]->GetUVPtStruct();
3371 for ( unsigned j = 0; j < uvVec.size(); ++j )
3373 TSmoothNode & sNode = smooNoMap[ uvVec[j].node ];
3374 sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v );
3378 // define refernce orientation in 2D
3379 TNo2SmooNoMap::iterator n2sn = smooNoMap.begin();
3380 for ( ; n2sn != smooNoMap.end(); ++n2sn )
3381 if ( !n2sn->second._triangles.empty() )
3383 if ( n2sn == smooNoMap.end() ) return;
3384 const TSmoothNode & sampleNode = n2sn->second;
3385 const bool refForward = ( sampleNode._triangles[0].IsForward( sampleNode._uv ));
3389 for ( int iLoop = 0; iLoop < 5; ++iLoop )
3391 for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
3393 TSmoothNode& sNode = n2sn->second;
3394 if ( sNode._triangles.empty() )
3395 continue; // not movable node
3399 for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
3400 newUV += sNode._triangles[i]._n1->_uv;
3401 newUV /= sNode._triangles.size();
3403 // check validity of the newUV
3404 bool isValid = true;
3405 for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
3406 isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
3413 // Set new XYZ to the smoothed nodes
3415 Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
3417 for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
3419 TSmoothNode& sNode = n2sn->second;
3420 if ( sNode._triangles.empty() )
3421 continue; // not movable node
3423 SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first );
3424 gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() );
3425 meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
3428 node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() )));
3431 // Move medium nodes in quadratic mesh
3432 if ( _quadraticMesh )
3434 const TLinkNodeMap& links = myHelper->GetTLinkNodeMap();
3435 TLinkNodeMap::const_iterator linkIt = links.begin();
3436 for ( ; linkIt != links.end(); ++linkIt )
3438 const SMESH_TLink& link = linkIt->first;
3439 SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( linkIt->second );
3441 if ( node->getshapeId() != myHelper->GetSubShapeID() )
3442 continue; // medium node is on EDGE or VERTEX
3444 gp_XY uv1 = myHelper->GetNodeUV( geomFace, link.node1(), node );
3445 gp_XY uv2 = myHelper->GetNodeUV( geomFace, link.node2(), node );
3447 gp_XY uv = myHelper->GetMiddleUV( surface, uv1, uv2 );
3448 node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
3450 gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
3451 meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );