00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "StdMeshers_Quadrangle_2D.hxx"
00028
00029 #include "StdMeshers_FaceSide.hxx"
00030
00031 #include "StdMeshers_QuadrangleParams.hxx"
00032
00033 #include "SMESH_Gen.hxx"
00034 #include "SMESH_Mesh.hxx"
00035 #include "SMESH_subMesh.hxx"
00036 #include "SMESH_MesherHelper.hxx"
00037 #include "SMESH_Block.hxx"
00038 #include "SMESH_Comment.hxx"
00039
00040 #include "SMDS_MeshElement.hxx"
00041 #include "SMDS_MeshNode.hxx"
00042 #include "SMDS_EdgePosition.hxx"
00043 #include "SMDS_FacePosition.hxx"
00044
00045 #include <BRep_Tool.hxx>
00046 #include <Geom_Surface.hxx>
00047 #include <NCollection_DefineArray2.hxx>
00048 #include <Precision.hxx>
00049 #include <TColStd_SequenceOfReal.hxx>
00050 #include <TColStd_SequenceOfInteger.hxx>
00051 #include <TColgp_SequenceOfXY.hxx>
00052 #include <TopExp.hxx>
00053 #include <TopExp_Explorer.hxx>
00054 #include <TopTools_ListIteratorOfListOfShape.hxx>
00055 #include <TopTools_MapOfShape.hxx>
00056 #include <TopoDS.hxx>
00057
00058 #include "utilities.h"
00059 #include "Utils_ExceptHandlers.hxx"
00060
00061 #ifndef StdMeshers_Array2OfNode_HeaderFile
00062 #define StdMeshers_Array2OfNode_HeaderFile
00063 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
00064 DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
00065 DEFINE_ARRAY2(StdMeshers_Array2OfNode,
00066 StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
00067 #endif
00068
00069 using namespace std;
00070
00071 typedef gp_XY gp_UV;
00072 typedef SMESH_Comment TComm;
00073
00074
00078
00079
00080 StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
00081 SMESH_Gen* gen)
00082 : SMESH_2D_Algo(hypId, studyId, gen)
00083 {
00084 MESSAGE("StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D");
00085 _name = "Quadrangle_2D";
00086 _shapeType = (1 << TopAbs_FACE);
00087 _compatibleHypothesis.push_back("QuadrangleParams");
00088 _compatibleHypothesis.push_back("QuadranglePreference");
00089 _compatibleHypothesis.push_back("TrianglePreference");
00090 myHelper = 0;
00091 }
00092
00093
00097
00098
00099 StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D()
00100 {
00101 MESSAGE("StdMeshers_Quadrangle_2D::~StdMeshers_Quadrangle_2D");
00102 }
00103
00104
00108
00109
00110 bool StdMeshers_Quadrangle_2D::CheckHypothesis
00111 (SMESH_Mesh& aMesh,
00112 const TopoDS_Shape& aShape,
00113 SMESH_Hypothesis::Hypothesis_Status& aStatus)
00114 {
00115 bool isOk = true;
00116 aStatus = SMESH_Hypothesis::HYP_OK;
00117
00118 const list <const SMESHDS_Hypothesis * >& hyps =
00119 GetUsedHypothesis(aMesh, aShape, false);
00120 const SMESHDS_Hypothesis * aHyp = 0;
00121
00122 myTriaVertexID = -1;
00123 myQuadType = QUAD_STANDARD;
00124 myQuadranglePreference = false;
00125 myTrianglePreference = false;
00126
00127 bool isFirstParams = true;
00128
00129
00130 if (hyps.size() > 0) {
00131 aHyp = hyps.front();
00132 if (strcmp("QuadrangleParams", aHyp->GetName()) == 0) {
00133 const StdMeshers_QuadrangleParams* aHyp1 =
00134 (const StdMeshers_QuadrangleParams*)aHyp;
00135 myTriaVertexID = aHyp1->GetTriaVertex();
00136 myQuadType = aHyp1->GetQuadType();
00137 if (myQuadType == QUAD_QUADRANGLE_PREF ||
00138 myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
00139 myQuadranglePreference = true;
00140 else if (myQuadType == QUAD_TRIANGLE_PREF)
00141 myTrianglePreference = true;
00142 }
00143 else if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) {
00144 isFirstParams = false;
00145 myQuadranglePreference = true;
00146 }
00147 else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){
00148 isFirstParams = false;
00149 myTrianglePreference = true;
00150 }
00151 else {
00152 isFirstParams = false;
00153 }
00154 }
00155
00156
00157 if (hyps.size() > 1) {
00158 aHyp = hyps.back();
00159 if (isFirstParams) {
00160 if (strcmp("QuadranglePreference", aHyp->GetName()) == 0) {
00161 myQuadranglePreference = true;
00162 myTrianglePreference = false;
00163 myQuadType = QUAD_STANDARD;
00164 }
00165 else if (strcmp("TrianglePreference", aHyp->GetName()) == 0){
00166 myQuadranglePreference = false;
00167 myTrianglePreference = true;
00168 myQuadType = QUAD_STANDARD;
00169 }
00170 }
00171 else {
00172 const StdMeshers_QuadrangleParams* aHyp2 =
00173 (const StdMeshers_QuadrangleParams*)aHyp;
00174 myTriaVertexID = aHyp2->GetTriaVertex();
00175
00176 if (!myQuadranglePreference && !myTrianglePreference) {
00177 myQuadType = aHyp2->GetQuadType();
00178 if (myQuadType == QUAD_QUADRANGLE_PREF ||
00179 myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
00180 myQuadranglePreference = true;
00181 else if (myQuadType == QUAD_TRIANGLE_PREF)
00182 myTrianglePreference = true;
00183 }
00184 }
00185 }
00186
00187 return isOk;
00188 }
00189
00190
00194
00195
00196 bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh& aMesh,
00197 const TopoDS_Shape& aShape)
00198 {
00199
00200
00201
00202 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00203 aMesh.GetSubMesh(aShape);
00204
00205 SMESH_MesherHelper helper (aMesh);
00206 myHelper = &helper;
00207
00208 _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
00209 myNeedSmooth = false;
00210
00211 FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
00212 std::auto_ptr<FaceQuadStruct> quadDeleter (quad);
00213 if (!quad)
00214 return false;
00215
00216 if (myQuadranglePreference) {
00217 int n1 = quad->side[0]->NbPoints();
00218 int n2 = quad->side[1]->NbPoints();
00219 int n3 = quad->side[2]->NbPoints();
00220 int n4 = quad->side[3]->NbPoints();
00221 int nfull = n1+n2+n3+n4;
00222 int ntmp = nfull/2;
00223 ntmp = ntmp*2;
00224 if (nfull == ntmp && ((n1 != n3) || (n2 != n4))) {
00225
00226 bool ok = ComputeQuadPref(aMesh, aShape, quad);
00227 if ( ok && myNeedSmooth )
00228 Smooth( quad );
00229 return ok;
00230 }
00231 }
00232 else if (myQuadType == QUAD_REDUCED) {
00233 int n1 = quad->side[0]->NbPoints();
00234 int n2 = quad->side[1]->NbPoints();
00235 int n3 = quad->side[2]->NbPoints();
00236 int n4 = quad->side[3]->NbPoints();
00237 int n13 = n1 - n3;
00238 int n24 = n2 - n4;
00239 int n13tmp = n13/2; n13tmp = n13tmp*2;
00240 int n24tmp = n24/2; n24tmp = n24tmp*2;
00241 if ((n1 == n3 && n2 != n4 && n24tmp == n24) ||
00242 (n2 == n4 && n1 != n3 && n13tmp == n13)) {
00243 bool ok = ComputeReduced(aMesh, aShape, quad);
00244 if ( ok && myNeedSmooth )
00245 Smooth( quad );
00246 return ok;
00247 }
00248 }
00249
00250
00251
00252 if (!SetNormalizedGrid(aMesh, aShape, quad))
00253 return false;
00254
00255
00256
00257 int nbdown = quad->side[0]->NbPoints();
00258 int nbup = quad->side[2]->NbPoints();
00259
00260 int nbright = quad->side[1]->NbPoints();
00261 int nbleft = quad->side[3]->NbPoints();
00262
00263 int nbhoriz = Min(nbdown, nbup);
00264 int nbvertic = Min(nbright, nbleft);
00265
00266 const TopoDS_Face& F = TopoDS::Face(aShape);
00267 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
00268
00269
00270 int i, j, geomFaceID = meshDS->ShapeToIndex(F);
00271 for (i = 1; i < nbhoriz - 1; i++) {
00272 for (j = 1; j < nbvertic - 1; j++) {
00273 int ij = j * nbhoriz + i;
00274 double u = quad->uv_grid[ij].u;
00275 double v = quad->uv_grid[ij].v;
00276 gp_Pnt P = S->Value(u, v);
00277 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
00278 meshDS->SetNodeOnFace(node, geomFaceID, u, v);
00279 quad->uv_grid[ij].node = node;
00280 }
00281 }
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 i = 0;
00298 int ilow = 0;
00299 int iup = nbhoriz - 1;
00300 if (quad->isEdgeOut[3]) { ilow++; } else { if (quad->isEdgeOut[1]) iup--; }
00301
00302 int jlow = 0;
00303 int jup = nbvertic - 1;
00304 if (quad->isEdgeOut[0]) { jlow++; } else { if (quad->isEdgeOut[2]) jup--; }
00305
00306
00307 for (i = ilow; i < iup; i++) {
00308 for (j = jlow; j < jup; j++) {
00309 const SMDS_MeshNode *a, *b, *c, *d;
00310 a = quad->uv_grid[j * nbhoriz + i].node;
00311 b = quad->uv_grid[j * nbhoriz + i + 1].node;
00312 c = quad->uv_grid[(j + 1) * nbhoriz + i + 1].node;
00313 d = quad->uv_grid[(j + 1) * nbhoriz + i].node;
00314 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
00315 if (face) {
00316 meshDS->SetMeshElementOnShape(face, geomFaceID);
00317 }
00318 }
00319 }
00320
00321 const vector<UVPtStruct>& uv_e0 = quad->side[0]->GetUVPtStruct(true,0);
00322 const vector<UVPtStruct>& uv_e1 = quad->side[1]->GetUVPtStruct(false,1);
00323 const vector<UVPtStruct>& uv_e2 = quad->side[2]->GetUVPtStruct(true,1);
00324 const vector<UVPtStruct>& uv_e3 = quad->side[3]->GetUVPtStruct(false,0);
00325
00326 if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
00327 return error(COMPERR_BAD_INPUT_MESH);
00328
00329 double eps = Precision::Confusion();
00330
00331
00332
00333 if (quad->isEdgeOut[0]) {
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345 int g = 0;
00346
00347
00348 int stop = nbdown - 1;
00349
00350 if (quad->isEdgeOut[1]) stop--;
00351
00352
00353
00354 for (i = 0; i < stop; i++) {
00355 const SMDS_MeshNode *a, *b, *c, *d;
00356 a = uv_e0[i].node;
00357 b = uv_e0[i + 1].node;
00358 gp_Pnt pb (b->X(), b->Y(), b->Z());
00359
00360
00361 int near = g;
00362 if (i == stop - 1) {
00363
00364 near = iup;
00365 c = quad->uv_grid[nbhoriz + iup].node;
00366 }
00367 else {
00368
00369 double mind = RealLast();
00370 for (int k = g; k <= iup; k++) {
00371
00372 const SMDS_MeshNode *nk;
00373 if (k < ilow)
00374 nk = uv_e3[1].node;
00375 else
00376 nk = quad->uv_grid[nbhoriz + k].node;
00377
00378 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
00379 double dist = pb.Distance(pnk);
00380 if (dist < mind - eps) {
00381 c = nk;
00382 near = k;
00383 mind = dist;
00384 } else {
00385 break;
00386 }
00387 }
00388 }
00389
00390 if (near == g) {
00391 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
00392 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00393 }
00394 else {
00395 if (near - 1 < ilow)
00396 d = uv_e3[1].node;
00397 else
00398 d = quad->uv_grid[nbhoriz + near - 1].node;
00399
00400
00401 if (!myTrianglePreference){
00402 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
00403 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00404 }
00405 else {
00406 SplitQuad(meshDS, geomFaceID, a, b, c, d);
00407 }
00408
00409
00410 if (near - 1 > g) {
00411 for (int k = near - 1; k > g; k--) {
00412 c = quad->uv_grid[nbhoriz + k].node;
00413 if (k - 1 < ilow)
00414 d = uv_e3[1].node;
00415 else
00416 d = quad->uv_grid[nbhoriz + k - 1].node;
00417 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
00418 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00419 }
00420 }
00421 g = near;
00422 }
00423 }
00424 } else {
00425 if (quad->isEdgeOut[2]) {
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438 int g = nbhoriz - 1;
00439
00440 int stop = 0;
00441
00442 if (quad->isEdgeOut[3]) stop++;
00443
00444
00445
00446 for (i = nbup - 1; i > stop; i--) {
00447 const SMDS_MeshNode *a, *b, *c, *d;
00448 a = uv_e2[i].node;
00449 b = uv_e2[i - 1].node;
00450 gp_Pnt pb (b->X(), b->Y(), b->Z());
00451
00452
00453 int near = g;
00454 if (i == stop + 1) {
00455 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + ilow].node;
00456 near = ilow;
00457 } else {
00458
00459 double mind = RealLast();
00460 for (int k = g; k >= ilow; k--) {
00461 const SMDS_MeshNode *nk;
00462 if (k > iup)
00463 nk = uv_e1[nbright - 2].node;
00464 else
00465 nk = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
00466 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
00467 double dist = pb.Distance(pnk);
00468 if (dist < mind - eps) {
00469 c = nk;
00470 near = k;
00471 mind = dist;
00472 } else {
00473 break;
00474 }
00475 }
00476 }
00477
00478 if (near == g) {
00479 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
00480 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00481 }
00482 else {
00483 if (near + 1 > iup)
00484 d = uv_e1[nbright - 2].node;
00485 else
00486 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + near + 1].node;
00487
00488 if (!myTrianglePreference){
00489 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
00490 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00491 }
00492 else {
00493 SplitQuad(meshDS, geomFaceID, a, b, c, d);
00494 }
00495
00496 if (near + 1 < g) {
00497 for (int k = near + 1; k < g; k++) {
00498 c = quad->uv_grid[nbhoriz*(nbvertic - 2) + k].node;
00499 if (k + 1 > iup)
00500 d = uv_e1[nbright - 2].node;
00501 else
00502 d = quad->uv_grid[nbhoriz*(nbvertic - 2) + k + 1].node;
00503 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
00504 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00505 }
00506 }
00507 g = near;
00508 }
00509 }
00510 }
00511 }
00512
00513
00514 if (quad->isEdgeOut[1]) {
00515
00516 int g = 0;
00517 int stop = nbright - 1;
00518 if (quad->isEdgeOut[2]) stop--;
00519 for (i = 0; i < stop; i++) {
00520 const SMDS_MeshNode *a, *b, *c, *d;
00521 a = uv_e1[i].node;
00522 b = uv_e1[i + 1].node;
00523 gp_Pnt pb (b->X(), b->Y(), b->Z());
00524
00525
00526 int near = g;
00527 if (i == stop - 1) {
00528 c = quad->uv_grid[nbhoriz*(jup + 1) - 2].node;
00529 near = jup;
00530 } else {
00531 double mind = RealLast();
00532 for (int k = g; k <= jup; k++) {
00533 const SMDS_MeshNode *nk;
00534 if (k < jlow)
00535 nk = uv_e0[nbdown - 2].node;
00536 else
00537 nk = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
00538 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
00539 double dist = pb.Distance(pnk);
00540 if (dist < mind - eps) {
00541 c = nk;
00542 near = k;
00543 mind = dist;
00544 } else {
00545 break;
00546 }
00547 }
00548 }
00549
00550 if (near == g) {
00551 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
00552 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00553 }
00554 else {
00555 if (near - 1 < jlow)
00556 d = uv_e0[nbdown - 2].node;
00557 else
00558 d = quad->uv_grid[nbhoriz*near - 2].node;
00559
00560
00561 if (!myTrianglePreference){
00562 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
00563 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00564 }
00565 else {
00566 SplitQuad(meshDS, geomFaceID, a, b, c, d);
00567 }
00568
00569 if (near - 1 > g) {
00570 for (int k = near - 1; k > g; k--) {
00571 c = quad->uv_grid[nbhoriz*(k + 1) - 2].node;
00572 if (k - 1 < jlow)
00573 d = uv_e0[nbdown - 2].node;
00574 else
00575 d = quad->uv_grid[nbhoriz*k - 2].node;
00576 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
00577 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00578 }
00579 }
00580 g = near;
00581 }
00582 }
00583 } else {
00584 if (quad->isEdgeOut[3]) {
00585
00586 int g = nbvertic - 1;
00587 int stop = 0;
00588 if (quad->isEdgeOut[0]) stop++;
00589 for (i = nbleft - 1; i > stop; i--) {
00590 const SMDS_MeshNode *a, *b, *c, *d;
00591 a = uv_e3[i].node;
00592 b = uv_e3[i - 1].node;
00593 gp_Pnt pb (b->X(), b->Y(), b->Z());
00594
00595
00596 int near = g;
00597 if (i == stop + 1) {
00598 c = quad->uv_grid[nbhoriz*jlow + 1].node;
00599 near = jlow;
00600 } else {
00601 double mind = RealLast();
00602 for (int k = g; k >= jlow; k--) {
00603 const SMDS_MeshNode *nk;
00604 if (k > jup)
00605 nk = uv_e2[1].node;
00606 else
00607 nk = quad->uv_grid[nbhoriz*k + 1].node;
00608 gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
00609 double dist = pb.Distance(pnk);
00610 if (dist < mind - eps) {
00611 c = nk;
00612 near = k;
00613 mind = dist;
00614 } else {
00615 break;
00616 }
00617 }
00618 }
00619
00620 if (near == g) {
00621 SMDS_MeshFace* face = myHelper->AddFace(a, b, c);
00622 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00623 }
00624 else {
00625 if (near + 1 > jup)
00626 d = uv_e2[1].node;
00627 else
00628 d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
00629
00630 if (!myTrianglePreference){
00631 SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
00632 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00633 }
00634 else {
00635 SplitQuad(meshDS, geomFaceID, a, b, c, d);
00636 }
00637
00638 if (near + 1 < g) {
00639 for (int k = near + 1; k < g; k++) {
00640 c = quad->uv_grid[nbhoriz*k + 1].node;
00641 if (k + 1 > jup)
00642 d = uv_e2[1].node;
00643 else
00644 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
00645 SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
00646 if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
00647 }
00648 }
00649 g = near;
00650 }
00651 }
00652 }
00653 }
00654
00655 if ( myNeedSmooth )
00656 Smooth( quad );
00657
00658 bool isOk = true;
00659 return isOk;
00660 }
00661
00662
00663
00667
00668
00669 bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh& aMesh,
00670 const TopoDS_Shape& aShape,
00671 MapShapeNbElems& aResMap)
00672
00673 {
00674 aMesh.GetSubMesh(aShape);
00675
00676 std::vector<int> aNbNodes(4);
00677 bool IsQuadratic = false;
00678 if (!CheckNbEdgesForEvaluate(aMesh, aShape, aResMap, aNbNodes, IsQuadratic)) {
00679 std::vector<int> aResVec(SMDSEntity_Last);
00680 for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00681 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00682 aResMap.insert(std::make_pair(sm,aResVec));
00683 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00684 smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
00685 return false;
00686 }
00687
00688 if (myQuadranglePreference) {
00689 int n1 = aNbNodes[0];
00690 int n2 = aNbNodes[1];
00691 int n3 = aNbNodes[2];
00692 int n4 = aNbNodes[3];
00693 int nfull = n1+n2+n3+n4;
00694 int ntmp = nfull/2;
00695 ntmp = ntmp*2;
00696 if (nfull==ntmp && ((n1!=n3) || (n2!=n4))) {
00697
00698 return EvaluateQuadPref(aMesh, aShape, aNbNodes, aResMap, IsQuadratic);
00699
00700 }
00701 }
00702
00703 int nbdown = aNbNodes[0];
00704 int nbup = aNbNodes[2];
00705
00706 int nbright = aNbNodes[1];
00707 int nbleft = aNbNodes[3];
00708
00709 int nbhoriz = Min(nbdown, nbup);
00710 int nbvertic = Min(nbright, nbleft);
00711
00712 int dh = Max(nbdown, nbup) - nbhoriz;
00713 int dv = Max(nbright, nbleft) - nbvertic;
00714
00715
00716
00717
00718
00719
00720 int nbNodes = (nbhoriz-2)*(nbvertic-2);
00721
00722 int nbFaces3 = dh + dv;
00723
00724
00725
00726 int nbFaces4 = (nbhoriz-1)*(nbvertic-1);
00727
00728 std::vector<int> aVec(SMDSEntity_Last);
00729 for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
00730 if (IsQuadratic) {
00731 aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
00732 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
00733 int nbbndedges = nbdown + nbup + nbright + nbleft -4;
00734 int nbintedges = (nbFaces4*4 + nbFaces3*3 - nbbndedges) / 2;
00735 aVec[SMDSEntity_Node] = nbNodes + nbintedges;
00736 if (aNbNodes.size()==5) {
00737 aVec[SMDSEntity_Quad_Triangle] = nbFaces3 + aNbNodes[3] -1;
00738 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4 - aNbNodes[3] +1;
00739 }
00740 }
00741 else {
00742 aVec[SMDSEntity_Node] = nbNodes;
00743 aVec[SMDSEntity_Triangle] = nbFaces3;
00744 aVec[SMDSEntity_Quadrangle] = nbFaces4;
00745 if (aNbNodes.size()==5) {
00746 aVec[SMDSEntity_Triangle] = nbFaces3 + aNbNodes[3] - 1;
00747 aVec[SMDSEntity_Quadrangle] = nbFaces4 - aNbNodes[3] + 1;
00748 }
00749 }
00750 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00751 aResMap.insert(std::make_pair(sm,aVec));
00752
00753 return true;
00754 }
00755
00756
00757
00761
00762
00763 static bool twoEdgesMeatAtVertex(const TopoDS_Edge& e1,
00764 const TopoDS_Edge& e2,
00765 SMESH_Mesh & mesh)
00766 {
00767 TopoDS_Vertex v;
00768 if (!TopExp::CommonVertex(e1, e2, v))
00769 return false;
00770 TopTools_ListIteratorOfListOfShape ancestIt(mesh.GetAncestors(v));
00771 for (; ancestIt.More() ; ancestIt.Next())
00772 if (ancestIt.Value().ShapeType() == TopAbs_EDGE)
00773 if (!e1.IsSame(ancestIt.Value()) && !e2.IsSame(ancestIt.Value()))
00774 return false;
00775 return true;
00776 }
00777
00778
00782
00783
00784 FaceQuadStruct* StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh & aMesh,
00785 const TopoDS_Shape & aShape)
00786
00787 {
00788 TopoDS_Face F = TopoDS::Face(aShape);
00789 if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
00790 const bool ignoreMediumNodes = _quadraticMesh;
00791
00792
00793 TopoDS_Vertex V;
00794 list< TopoDS_Edge > edges;
00795 list< int > nbEdgesInWire;
00796 int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
00797 if (nbWire != 1) {
00798 error(COMPERR_BAD_SHAPE, TComm("Wrong number of wires: ") << nbWire);
00799 return 0;
00800 }
00801 FaceQuadStruct* quad = new FaceQuadStruct;
00802 quad->uv_grid = 0;
00803 quad->side.reserve(nbEdgesInWire.front());
00804 quad->face = F;
00805
00806 int nbSides = 0;
00807 list< TopoDS_Edge >::iterator edgeIt = edges.begin();
00808 if (nbEdgesInWire.front() == 3)
00809 {
00810 SMESH_Comment comment;
00811 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
00812 if (myTriaVertexID == -1)
00813 {
00814 comment << "No Base vertex parameter provided for a trilateral geometrical face";
00815 }
00816 else
00817 {
00818 TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
00819 if (!V.IsNull()) {
00820 TopoDS_Edge E1,E2,E3;
00821 for (; edgeIt != edges.end(); ++edgeIt) {
00822 TopoDS_Edge E = *edgeIt;
00823 TopoDS_Vertex VF, VL;
00824 TopExp::Vertices(E, VF, VL, true);
00825 if (VF.IsSame(V))
00826 E1 = E;
00827 else if (VL.IsSame(V))
00828 E3 = E;
00829 else
00830 E2 = E;
00831 }
00832 if (!E1.IsNull() && !E2.IsNull() && !E3.IsNull())
00833 {
00834 quad->side.push_back(new StdMeshers_FaceSide(F, E1, &aMesh, true, ignoreMediumNodes));
00835 quad->side.push_back(new StdMeshers_FaceSide(F, E2, &aMesh, true, ignoreMediumNodes));
00836 quad->side.push_back(new StdMeshers_FaceSide(F, E3, &aMesh, false,ignoreMediumNodes));
00837 const vector<UVPtStruct>& UVPSleft = quad->side[0]->GetUVPtStruct(true,0);
00838 quad->side[1]->GetUVPtStruct(false,1);
00839 quad->side[2]->GetUVPtStruct(true,1);
00840 const SMDS_MeshNode* aNode = UVPSleft[0].node;
00841 gp_Pnt2d aPnt2d(UVPSleft[0].u, UVPSleft[0].v);
00842 quad->side.push_back(new StdMeshers_FaceSide(aNode, aPnt2d, quad->side[1]));
00843 return quad;
00844 }
00845 }
00846 comment << "Invalid Base vertex parameter: " << myTriaVertexID << " is not among [";
00847 TopTools_MapOfShape vMap;
00848 for (TopExp_Explorer v(aShape, TopAbs_VERTEX); v.More(); v.Next())
00849 if (vMap.Add(v.Current()))
00850 comment << meshDS->ShapeToIndex(v.Current()) << (vMap.Extent()==3 ? "]" : ", ");
00851 }
00852 error(comment);
00853 delete quad;
00854 return quad = 0;
00855 }
00856 else if (nbEdgesInWire.front() == 4)
00857 {
00858 for (; edgeIt != edges.end(); ++edgeIt, nbSides++)
00859 quad->side.push_back(new StdMeshers_FaceSide(F, *edgeIt, &aMesh,
00860 nbSides<TOP_SIDE, ignoreMediumNodes));
00861 }
00862 else if (nbEdgesInWire.front() > 4)
00863 {
00864 list< TopoDS_Edge > sideEdges;
00865 vector< int > degenSides;
00866 while (!edges.empty()) {
00867 sideEdges.clear();
00868 sideEdges.splice(sideEdges.end(), edges, edges.begin());
00869 bool sameSide = true;
00870 while (!edges.empty() && sameSide) {
00871 sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
00872 if (sameSide)
00873 sideEdges.splice(sideEdges.end(), edges, edges.begin());
00874 }
00875 if (nbSides == 0) {
00876 sameSide = true;
00877 while (!edges.empty() && sameSide) {
00878 sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
00879 if (sameSide)
00880 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
00881 }
00882 }
00883 if ( sideEdges.size() == 1 && BRep_Tool::Degenerated( sideEdges.front() ))
00884 degenSides.push_back( nbSides );
00885
00886 quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
00887 nbSides<TOP_SIDE, ignoreMediumNodes));
00888 ++nbSides;
00889 }
00890 if ( !degenSides.empty() && nbSides - degenSides.size() == 4 )
00891 {
00892 myNeedSmooth = true;
00893 for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
00894 quad->side[i]->Reverse();
00895
00896 for ( int i = degenSides.size()-1; i > -1; --i )
00897 {
00898 StdMeshers_FaceSide * & degenSide = quad->side[ degenSides[ i ]];
00899 delete degenSide;
00900 quad->side.erase( vector<StdMeshers_FaceSide*>::iterator( & degenSide ));
00901 }
00902 for ( unsigned i = TOP_SIDE; i < quad->side.size(); ++i )
00903 quad->side[i]->Reverse();
00904
00905 nbSides -= degenSides.size();
00906 }
00907
00908 if (nbSides < 4) {
00909
00910 { FaceQuadStruct cleaner(*quad); }
00911 quad->side.clear();
00912 quad->side.reserve(nbEdgesInWire.front());
00913 nbSides = 0;
00914
00915 SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
00916 while (!edges.empty()) {
00917 sideEdges.clear();
00918 sideEdges.splice(sideEdges.end(), edges, edges.begin());
00919 bool sameSide = true;
00920 while (!edges.empty() && sameSide) {
00921 sameSide =
00922 SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
00923 twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
00924 if (sameSide)
00925 sideEdges.splice(sideEdges.end(), edges, edges.begin());
00926 }
00927 if (nbSides == 0) {
00928 sameSide = true;
00929 while (!edges.empty() && sameSide) {
00930 sameSide =
00931 SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
00932 twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
00933 if (sameSide)
00934 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
00935 }
00936 }
00937 quad->side.push_back(new StdMeshers_FaceSide(F, sideEdges, &aMesh,
00938 nbSides<TOP_SIDE, ignoreMediumNodes));
00939 ++nbSides;
00940 }
00941 }
00942 }
00943 if (nbSides != 4) {
00944 #ifdef _DEBUG_
00945 MESSAGE ("StdMeshers_Quadrangle_2D. Edge IDs of " << nbSides << " sides:\n");
00946 for (int i = 0; i < nbSides; ++i) {
00947 MESSAGE (" (");
00948 for (int e = 0; e < quad->side[i]->NbEdges(); ++e)
00949 MESSAGE (myHelper->GetMeshDS()->ShapeToIndex(quad->side[i]->Edge(e)) << " ");
00950 MESSAGE (")\n");
00951 }
00952
00953 #endif
00954 if (!nbSides)
00955 nbSides = nbEdgesInWire.front();
00956 error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
00957 delete quad;
00958 quad = 0;
00959 }
00960
00961 return quad;
00962 }
00963
00964
00965
00969
00970
00971 bool StdMeshers_Quadrangle_2D::CheckNbEdgesForEvaluate(SMESH_Mesh& aMesh,
00972 const TopoDS_Shape & aShape,
00973 MapShapeNbElems& aResMap,
00974 std::vector<int>& aNbNodes,
00975 bool& IsQuadratic)
00976
00977 {
00978 const TopoDS_Face & F = TopoDS::Face(aShape);
00979
00980
00981 TopoDS_Vertex V;
00982 list< TopoDS_Edge > edges;
00983 list< int > nbEdgesInWire;
00984 int nbWire = SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
00985 if (nbWire != 1) {
00986 return false;
00987 }
00988
00989 aNbNodes.resize(4);
00990
00991 int nbSides = 0;
00992 list< TopoDS_Edge >::iterator edgeIt = edges.begin();
00993 SMESH_subMesh * sm = aMesh.GetSubMesh(*edgeIt);
00994 MapShapeNbElemsItr anIt = aResMap.find(sm);
00995 if (anIt==aResMap.end()) {
00996 return false;
00997 }
00998 std::vector<int> aVec = (*anIt).second;
00999 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
01000 if (nbEdgesInWire.front() == 3) {
01001 if (myTriaVertexID>0) {
01002 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
01003 TopoDS_Vertex V = TopoDS::Vertex(meshDS->IndexToShape(myTriaVertexID));
01004 if (!V.IsNull()) {
01005 TopoDS_Edge E1,E2,E3;
01006 for (; edgeIt != edges.end(); ++edgeIt) {
01007 TopoDS_Edge E = TopoDS::Edge(*edgeIt);
01008 TopoDS_Vertex VF, VL;
01009 TopExp::Vertices(E, VF, VL, true);
01010 if (VF.IsSame(V))
01011 E1 = E;
01012 else if (VL.IsSame(V))
01013 E3 = E;
01014 else
01015 E2 = E;
01016 }
01017 SMESH_subMesh * sm = aMesh.GetSubMesh(E1);
01018 MapShapeNbElemsItr anIt = aResMap.find(sm);
01019 if (anIt==aResMap.end()) return false;
01020 std::vector<int> aVec = (*anIt).second;
01021 if (IsQuadratic)
01022 aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2;
01023 else
01024 aNbNodes[0] = aVec[SMDSEntity_Node] + 2;
01025 sm = aMesh.GetSubMesh(E2);
01026 anIt = aResMap.find(sm);
01027 if (anIt==aResMap.end()) return false;
01028 aVec = (*anIt).second;
01029 if (IsQuadratic)
01030 aNbNodes[1] = (aVec[SMDSEntity_Node]-1)/2 + 2;
01031 else
01032 aNbNodes[1] = aVec[SMDSEntity_Node] + 2;
01033 sm = aMesh.GetSubMesh(E3);
01034 anIt = aResMap.find(sm);
01035 if (anIt==aResMap.end()) return false;
01036 aVec = (*anIt).second;
01037 if (IsQuadratic)
01038 aNbNodes[2] = (aVec[SMDSEntity_Node]-1)/2 + 2;
01039 else
01040 aNbNodes[2] = aVec[SMDSEntity_Node] + 2;
01041 aNbNodes[3] = aNbNodes[1];
01042 aNbNodes.resize(5);
01043 nbSides = 4;
01044 }
01045 }
01046 }
01047 if (nbEdgesInWire.front() == 4) {
01048 for (; edgeIt != edges.end(); edgeIt++) {
01049 SMESH_subMesh * sm = aMesh.GetSubMesh(*edgeIt);
01050 MapShapeNbElemsItr anIt = aResMap.find(sm);
01051 if (anIt==aResMap.end()) {
01052 return false;
01053 }
01054 std::vector<int> aVec = (*anIt).second;
01055 if (IsQuadratic)
01056 aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2;
01057 else
01058 aNbNodes[nbSides] = aVec[SMDSEntity_Node] + 2;
01059 nbSides++;
01060 }
01061 }
01062 else if (nbEdgesInWire.front() > 4) {
01063 list< TopoDS_Edge > sideEdges;
01064 while (!edges.empty()) {
01065 sideEdges.clear();
01066 sideEdges.splice(sideEdges.end(), edges, edges.begin());
01067 bool sameSide = true;
01068 while (!edges.empty() && sameSide) {
01069 sameSide = SMESH_Algo::IsContinuous(sideEdges.back(), edges.front());
01070 if (sameSide)
01071 sideEdges.splice(sideEdges.end(), edges, edges.begin());
01072 }
01073 if (nbSides == 0) {
01074 sameSide = true;
01075 while (!edges.empty() && sameSide) {
01076 sameSide = SMESH_Algo::IsContinuous(sideEdges.front(), edges.back());
01077 if (sameSide)
01078 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
01079 }
01080 }
01081 list<TopoDS_Edge>::iterator ite = sideEdges.begin();
01082 aNbNodes[nbSides] = 1;
01083 for (; ite!=sideEdges.end(); ite++) {
01084 SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
01085 MapShapeNbElemsItr anIt = aResMap.find(sm);
01086 if (anIt==aResMap.end()) {
01087 return false;
01088 }
01089 std::vector<int> aVec = (*anIt).second;
01090 if (IsQuadratic)
01091 aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
01092 else
01093 aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
01094 }
01095 ++nbSides;
01096 }
01097
01098 if (nbSides < 4) {
01099 nbSides = 0;
01100 SMESH_Block::GetOrderedEdges (F, V, edges, nbEdgesInWire);
01101 while (!edges.empty()) {
01102 sideEdges.clear();
01103 sideEdges.splice(sideEdges.end(), edges, edges.begin());
01104 bool sameSide = true;
01105 while (!edges.empty() && sameSide) {
01106 sameSide =
01107 SMESH_Algo::IsContinuous(sideEdges.back(), edges.front()) &&
01108 twoEdgesMeatAtVertex(sideEdges.back(), edges.front(), aMesh);
01109 if (sameSide)
01110 sideEdges.splice(sideEdges.end(), edges, edges.begin());
01111 }
01112 if (nbSides == 0) {
01113 sameSide = true;
01114 while (!edges.empty() && sameSide) {
01115 sameSide =
01116 SMESH_Algo::IsContinuous(sideEdges.front(), edges.back()) &&
01117 twoEdgesMeatAtVertex(sideEdges.front(), edges.back(), aMesh);
01118 if (sameSide)
01119 sideEdges.splice(sideEdges.begin(), edges, --edges.end());
01120 }
01121 }
01122 list<TopoDS_Edge>::iterator ite = sideEdges.begin();
01123 aNbNodes[nbSides] = 1;
01124 for (; ite!=sideEdges.end(); ite++) {
01125 SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
01126 MapShapeNbElemsItr anIt = aResMap.find(sm);
01127 if (anIt==aResMap.end()) {
01128 return false;
01129 }
01130 std::vector<int> aVec = (*anIt).second;
01131 if (IsQuadratic)
01132 aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
01133 else
01134 aNbNodes[nbSides] += aVec[SMDSEntity_Node] + 1;
01135 }
01136 ++nbSides;
01137 }
01138 }
01139 }
01140 if (nbSides != 4) {
01141 if (!nbSides)
01142 nbSides = nbEdgesInWire.front();
01143 error(COMPERR_BAD_SHAPE, TComm("Face must have 4 sides but not ") << nbSides);
01144 return false;
01145 }
01146
01147 return true;
01148 }
01149
01150
01151
01155
01156
01157 FaceQuadStruct *StdMeshers_Quadrangle_2D::CheckAnd2Dcompute
01158 (SMESH_Mesh & aMesh,
01159 const TopoDS_Shape & aShape,
01160 const bool CreateQuadratic)
01161 {
01162 _quadraticMesh = CreateQuadratic;
01163
01164 FaceQuadStruct *quad = CheckNbEdges(aMesh, aShape);
01165
01166 if (!quad) return 0;
01167
01168
01169 bool stat = SetNormalizedGrid(aMesh, aShape, quad);
01170 if (!stat) {
01171 if (quad) delete quad;
01172 quad = 0;
01173 }
01174
01175 return quad;
01176 }
01177
01178
01182
01183
01184 faceQuadStruct::~faceQuadStruct()
01185 {
01186 for (int i = 0; i < side.size(); i++) {
01187 if (side[i]) delete side[i];
01188 }
01189 if (uv_grid) delete [] uv_grid;
01190 }
01191
01192 namespace {
01193 inline const vector<UVPtStruct>& GetUVPtStructIn(FaceQuadStruct* quad, int i, int nbSeg)
01194 {
01195 bool isXConst = (i == BOTTOM_SIDE || i == TOP_SIDE);
01196 double constValue = (i == BOTTOM_SIDE || i == LEFT_SIDE) ? 0 : 1;
01197 return
01198 quad->isEdgeOut[i] ?
01199 quad->side[i]->SimulateUVPtStruct(nbSeg,isXConst,constValue) :
01200 quad->side[i]->GetUVPtStruct(isXConst,constValue);
01201 }
01202 inline gp_UV CalcUV(double x, double y,
01203 const gp_UV& a0,const gp_UV& a1,const gp_UV& a2,const gp_UV& a3,
01204 const gp_UV& p0,const gp_UV& p1,const gp_UV& p2,const gp_UV& p3)
01205 {
01206 return
01207 ((1 - y) * p0 + x * p1 + y * p2 + (1 - x) * p3 ) -
01208 ((1 - x) * (1 - y) * a0 + x * (1 - y) * a1 + x * y * a2 + (1 - x) * y * a3);
01209 }
01210 }
01211
01212
01216
01217
01218 bool StdMeshers_Quadrangle_2D::SetNormalizedGrid (SMESH_Mesh & aMesh,
01219 const TopoDS_Shape& aShape,
01220 FaceQuadStruct* & quad)
01221 {
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 int nbhoriz = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
01247 int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
01248
01249 quad->isEdgeOut[0] = (quad->side[0]->NbPoints() > quad->side[2]->NbPoints());
01250 quad->isEdgeOut[1] = (quad->side[1]->NbPoints() > quad->side[3]->NbPoints());
01251 quad->isEdgeOut[2] = (quad->side[2]->NbPoints() > quad->side[0]->NbPoints());
01252 quad->isEdgeOut[3] = (quad->side[3]->NbPoints() > quad->side[1]->NbPoints());
01253
01254 UVPtStruct *uv_grid = quad->uv_grid = new UVPtStruct[nbvertic * nbhoriz];
01255
01256 const vector<UVPtStruct>& uv_e0 = GetUVPtStructIn(quad, 0, nbhoriz - 1);
01257 const vector<UVPtStruct>& uv_e1 = GetUVPtStructIn(quad, 1, nbvertic - 1);
01258 const vector<UVPtStruct>& uv_e2 = GetUVPtStructIn(quad, 2, nbhoriz - 1);
01259 const vector<UVPtStruct>& uv_e3 = GetUVPtStructIn(quad, 3, nbvertic - 1);
01260
01261 if (uv_e0.empty() || uv_e1.empty() || uv_e2.empty() || uv_e3.empty())
01262
01263 return error(COMPERR_BAD_INPUT_MESH);
01264
01265 if ( myNeedSmooth )
01266 UpdateDegenUV( quad );
01267
01268
01269 if (! quad->isEdgeOut[0]) {
01270 int j = 0;
01271 for (int i = 0; i < nbhoriz; i++) {
01272 int ij = j * nbhoriz + i;
01273 uv_grid[ij].node = uv_e0[i].node;
01274 }
01275 }
01276 if (! quad->isEdgeOut[1]) {
01277 int i = nbhoriz - 1;
01278 for (int j = 0; j < nbvertic; j++) {
01279 int ij = j * nbhoriz + i;
01280 uv_grid[ij].node = uv_e1[j].node;
01281 }
01282 }
01283 if (! quad->isEdgeOut[2]) {
01284 int j = nbvertic - 1;
01285 for (int i = 0; i < nbhoriz; i++) {
01286 int ij = j * nbhoriz + i;
01287 uv_grid[ij].node = uv_e2[i].node;
01288 }
01289 }
01290 if (! quad->isEdgeOut[3]) {
01291 int i = 0;
01292 for (int j = 0; j < nbvertic; j++) {
01293 int ij = j * nbhoriz + i;
01294 uv_grid[ij].node = uv_e3[j].node;
01295 }
01296 }
01297
01298
01299 for (int i = 0; i < nbhoriz; i++) {
01300 for (int j = 0; j < nbvertic; j++) {
01301 int ij = j * nbhoriz + i;
01302
01303 double x0 = uv_e0[i].normParam;
01304 double x1 = uv_e2[i].normParam;
01305
01306 double y0 = uv_e3[j].normParam;
01307 double y1 = uv_e1[j].normParam;
01308
01309 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
01310 double y = y0 + x * (y1 - y0);
01311 uv_grid[ij].x = x;
01312 uv_grid[ij].y = y;
01313
01314
01315 }
01316 }
01317
01318
01319 gp_UV a0(uv_e0.front().u, uv_e0.front().v);
01320 gp_UV a1(uv_e0.back().u, uv_e0.back().v);
01321 gp_UV a2(uv_e2.back().u, uv_e2.back().v);
01322 gp_UV a3(uv_e2.front().u, uv_e2.front().v);
01323
01324 for (int i = 0; i < nbhoriz; i++) {
01325 for (int j = 0; j < nbvertic; j++) {
01326 int ij = j * nbhoriz + i;
01327 double x = uv_grid[ij].x;
01328 double y = uv_grid[ij].y;
01329 double param_0 = uv_e0[0].normParam + x * (uv_e0.back().normParam - uv_e0[0].normParam);
01330 double param_2 = uv_e2[0].normParam + x * (uv_e2.back().normParam - uv_e2[0].normParam);
01331 double param_1 = uv_e1[0].normParam + y * (uv_e1.back().normParam - uv_e1[0].normParam);
01332 double param_3 = uv_e3[0].normParam + y * (uv_e3.back().normParam - uv_e3[0].normParam);
01333
01334
01335 gp_UV p0 = quad->side[0]->Value2d(param_0).XY();
01336 gp_UV p1 = quad->side[1]->Value2d(param_1).XY();
01337 gp_UV p2 = quad->side[2]->Value2d(param_2).XY();
01338 gp_UV p3 = quad->side[3]->Value2d(param_3).XY();
01339
01340 gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
01341
01342 uv_grid[ij].u = uv.X();
01343 uv_grid[ij].v = uv.Y();
01344 }
01345 }
01346 return true;
01347 }
01348
01349
01350
01351
01352
01353
01354 static void ShiftQuad(FaceQuadStruct* quad, const int num, bool)
01355 {
01356 StdMeshers_FaceSide* side[4] = { quad->side[0], quad->side[1], quad->side[2], quad->side[3] };
01357 for (int i = BOTTOM_SIDE; i < NB_SIDES; ++i) {
01358 int id = (i + num) % NB_SIDES;
01359 bool wasForward = (i < TOP_SIDE);
01360 bool newForward = (id < TOP_SIDE);
01361 if (wasForward != newForward)
01362 side[ i ]->Reverse();
01363 quad->side[ id ] = side[ i ];
01364 }
01365 }
01366
01367
01368
01369
01370
01371
01372 static gp_UV CalcUV(double x0, double x1, double y0, double y1,
01373 FaceQuadStruct* quad,
01374 const gp_UV& a0, const gp_UV& a1,
01375 const gp_UV& a2, const gp_UV& a3)
01376 {
01377 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
01378 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
01379 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
01380 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
01381
01382 double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0));
01383 double y = y0 + x * (y1 - y0);
01384
01385 double param_b = uv_eb[0].normParam + x * (uv_eb.back().normParam - uv_eb[0].normParam);
01386 double param_t = uv_et[0].normParam + x * (uv_et.back().normParam - uv_et[0].normParam);
01387 double param_r = uv_er[0].normParam + y * (uv_er.back().normParam - uv_er[0].normParam);
01388 double param_l = uv_el[0].normParam + y * (uv_el.back().normParam - uv_el[0].normParam);
01389
01390 gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(param_b).XY();
01391 gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(param_r).XY();
01392 gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(param_t).XY();
01393 gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(param_l).XY();
01394
01395 gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
01396
01397 return uv;
01398 }
01399
01400
01401
01402
01403
01404
01405 static gp_UV CalcUV2(double x, double y,
01406 FaceQuadStruct* quad,
01407 const gp_UV& a0, const gp_UV& a1,
01408 const gp_UV& a2, const gp_UV& a3)
01409 {
01410 gp_UV p0 = quad->side[BOTTOM_SIDE]->Value2d(x).XY();
01411 gp_UV p1 = quad->side[RIGHT_SIDE ]->Value2d(y).XY();
01412 gp_UV p2 = quad->side[TOP_SIDE ]->Value2d(x).XY();
01413 gp_UV p3 = quad->side[LEFT_SIDE ]->Value2d(y).XY();
01414
01415 gp_UV uv = CalcUV(x,y, a0,a1,a2,a3, p0,p1,p2,p3);
01416
01417 return uv;
01418 }
01419
01420
01421
01425
01426
01427 bool StdMeshers_Quadrangle_2D::ComputeQuadPref (SMESH_Mesh & aMesh,
01428 const TopoDS_Shape& aShape,
01429 FaceQuadStruct* quad)
01430 {
01431
01432
01433
01434 bool OldVersion = false;
01435 if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
01436 OldVersion = true;
01437
01438 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
01439 const TopoDS_Face& F = TopoDS::Face(aShape);
01440 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
01441 bool WisF = true;
01442 int i,j,geomFaceID = meshDS->ShapeToIndex(F);
01443
01444 int nb = quad->side[0]->NbPoints();
01445 int nr = quad->side[1]->NbPoints();
01446 int nt = quad->side[2]->NbPoints();
01447 int nl = quad->side[3]->NbPoints();
01448 int dh = abs(nb-nt);
01449 int dv = abs(nr-nl);
01450
01451 if (dh>=dv) {
01452 if (nt>nb) {
01453
01454 ShiftQuad(quad,0,WisF);
01455 }
01456 else {
01457
01458 ShiftQuad(quad,2,WisF);
01459 }
01460 }
01461 else {
01462 if (nr>nl) {
01463
01464 ShiftQuad(quad,1,WisF);
01465 }
01466 else {
01467
01468 ShiftQuad(quad,3,WisF);
01469 }
01470 }
01471
01472 nb = quad->side[0]->NbPoints();
01473 nr = quad->side[1]->NbPoints();
01474 nt = quad->side[2]->NbPoints();
01475 nl = quad->side[3]->NbPoints();
01476 dh = abs(nb-nt);
01477 dv = abs(nr-nl);
01478 int nbh = Max(nb,nt);
01479 int nbv = Max(nr,nl);
01480 int addh = 0;
01481 int addv = 0;
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511 if (dh>dv) {
01512 addv = (dh-dv)/2;
01513 nbv = nbv + addv;
01514 }
01515 else {
01516 addh = (dv-dh)/2;
01517 nbh = nbh + addh;
01518 }
01519
01520 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
01521 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
01522 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
01523 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
01524
01525 if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
01526 return error(COMPERR_BAD_INPUT_MESH);
01527
01528 if ( myNeedSmooth )
01529 UpdateDegenUV( quad );
01530
01531
01532
01533 TColStd_SequenceOfReal npb, npr, npt, npl;
01534 for (i=0; i<nb; i++) {
01535 npb.Append(uv_eb[i].normParam);
01536
01537
01538
01539 }
01540 for (i=0; i<nr; i++) {
01541 npr.Append(uv_er[i].normParam);
01542 }
01543 for (i=0; i<nt; i++) {
01544 npt.Append(uv_et[i].normParam);
01545 }
01546 for (i=0; i<nl; i++) {
01547 npl.Append(uv_el[i].normParam);
01548 }
01549
01550 int dl,dr;
01551 if (OldVersion) {
01552
01553
01554 dr = nbv - nr;
01555 double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
01556 for (i=1; i<=dr; i++) {
01557 npr.InsertAfter(1,npr.Value(2)-dpr);
01558 }
01559
01560 dl = nbv - nl;
01561 dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
01562 for (i=1; i<=dl; i++) {
01563 npl.InsertAfter(1,npl.Value(2)-dpr);
01564 }
01565 }
01566
01567
01568
01569
01570
01571
01572 gp_XY a0(uv_eb.front().u, uv_eb.front().v);
01573 gp_XY a1(uv_eb.back().u, uv_eb.back().v);
01574 gp_XY a2(uv_et.back().u, uv_et.back().v);
01575 gp_XY a3(uv_et.front().u, uv_et.front().v);
01576
01577
01578
01579 int nnn = Min(nr,nl);
01580
01581
01582
01583 TColgp_SequenceOfXY UVL, UVR, UVT;
01584
01585 if (OldVersion) {
01586
01587 StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
01588
01589 for (j=1; j<=nl; j++)
01590 NodesL.SetValue(1,j,uv_el[j-1].node);
01591 if (dl>0) {
01592
01593 for (i=1; i<=dl; i++)
01594 NodesL.SetValue(i+1,nl,uv_et[i].node);
01595
01596 TColgp_SequenceOfXY UVtmp;
01597 for (i=1; i<=dl; i++) {
01598 double x0 = npt.Value(i+1);
01599 double x1 = x0;
01600
01601 double y0 = npl.Value(i+1);
01602 double y1 = npr.Value(i+1);
01603 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
01604 gp_Pnt P = S->Value(UV.X(),UV.Y());
01605 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01606 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01607 NodesL.SetValue(i+1,1,N);
01608 if (UVL.Length()<nbv-nnn) UVL.Append(UV);
01609
01610 for (j=2; j<nl; j++) {
01611 double y0 = npl.Value(dl+j);
01612 double y1 = npr.Value(dl+j);
01613 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
01614 gp_Pnt P = S->Value(UV.X(),UV.Y());
01615 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01616 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01617 NodesL.SetValue(i+1,j,N);
01618 if (i==dl) UVtmp.Append(UV);
01619 }
01620 }
01621 for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn; i++) {
01622 UVL.Append(UVtmp.Value(i));
01623 }
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633 for (i=1; i<=dl; i++) {
01634 for (j=1; j<nl; j++) {
01635 if (WisF) {
01636 SMDS_MeshFace* F =
01637 myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
01638 NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
01639 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01640 }
01641 else {
01642 SMDS_MeshFace* F =
01643 myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i,j+1),
01644 NodesL.Value(i+1,j+1), NodesL.Value(i+1,j));
01645 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01646 }
01647 }
01648 }
01649 }
01650 else {
01651
01652 for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn; i++) {
01653 UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
01654 }
01655 }
01656
01657
01658 StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
01659
01660 for (j=1; j<=nr; j++)
01661 NodesR.SetValue(1,j,uv_er[nr-j].node);
01662 if (dr>0) {
01663
01664 for (i=1; i<=dr; i++)
01665 NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
01666
01667 TColgp_SequenceOfXY UVtmp;
01668 for (i=1; i<=dr; i++) {
01669 double x0 = npt.Value(nt-i);
01670 double x1 = x0;
01671
01672 double y0 = npl.Value(i+1);
01673 double y1 = npr.Value(i+1);
01674 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
01675 gp_Pnt P = S->Value(UV.X(),UV.Y());
01676 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01677 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01678 NodesR.SetValue(i+1,nr,N);
01679 if (UVR.Length()<nbv-nnn) UVR.Append(UV);
01680
01681 for (j=2; j<nr; j++) {
01682 double y0 = npl.Value(nbv-j+1);
01683 double y1 = npr.Value(nbv-j+1);
01684 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
01685 gp_Pnt P = S->Value(UV.X(),UV.Y());
01686 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01687 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01688 NodesR.SetValue(i+1,j,N);
01689 if (i==dr) UVtmp.Prepend(UV);
01690 }
01691 }
01692 for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn; i++) {
01693 UVR.Append(UVtmp.Value(i));
01694 }
01695
01696 for (i=1; i<=dr; i++) {
01697 for (j=1; j<nr; j++) {
01698 if (WisF) {
01699 SMDS_MeshFace* F =
01700 myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
01701 NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
01702 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01703 }
01704 else {
01705 SMDS_MeshFace* F =
01706 myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i,j+1),
01707 NodesR.Value(i+1,j+1), NodesR.Value(i+1,j));
01708 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01709 }
01710 }
01711 }
01712 }
01713 else {
01714
01715 for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn; i++) {
01716 UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
01717 }
01718 }
01719
01720
01721 StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
01722
01723 for (i=1; i<=dl+1; i++)
01724 NodesC.SetValue(1,i,NodesL(i,1));
01725 for (i=2; i<=nl; i++)
01726 NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
01727
01728 for (i=1; i<=dr+1; i++)
01729 NodesC.SetValue(nb,i,NodesR(i,nr));
01730 for (i=1; i<nr; i++)
01731 NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
01732
01733 for (i=dl+2; i<nbh-dr; i++)
01734 NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
01735
01736 for (i=2; i<nb; i++)
01737 NodesC.SetValue(i,1,uv_eb[i-1].node);
01738
01739
01740
01741 for (i=2; i<nb; i++) {
01742 double x0 = npt.Value(dl+i);
01743 double x1 = x0;
01744 for (j=1; j<nnn; j++) {
01745 double y0 = npl.Value(nbv-nnn+j);
01746 double y1 = npr.Value(nbv-nnn+j);
01747 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
01748 gp_Pnt P = S->Value(UV.X(),UV.Y());
01749 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01750 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01751 NodesC.SetValue(i,nbv-nnn+j,N);
01752 if ( j==1 )
01753 UVT.Append( UV );
01754 }
01755 }
01756
01757
01758
01759
01760
01761
01762
01763 gp_UV A2 = UVR.Value(nbv-nnn);
01764 gp_UV A3 = UVL.Value(nbv-nnn);
01765 for (i=1; i<nbv-nnn; i++) {
01766 gp_UV p1 = UVR.Value(i);
01767 gp_UV p3 = UVL.Value(i);
01768 double y = i / double(nbv-nnn);
01769 for (j=2; j<nb; j++) {
01770 double x = npb.Value(j);
01771 gp_UV p0( uv_eb[j-1].u, uv_eb[j-1].v );
01772 gp_UV p2 = UVT.Value( j-1 );
01773 gp_UV UV = CalcUV(x, y, a0, a1, A2, A3, p0,p1,p2,p3 );
01774 gp_Pnt P = S->Value(UV.X(),UV.Y());
01775 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01776 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
01777 NodesC.SetValue(j,i+1,N);
01778 }
01779 }
01780
01781 for (i=1; i<nb; i++) {
01782 for (j=1; j<nbv; j++) {
01783 if (WisF) {
01784 SMDS_MeshFace* F =
01785 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
01786 NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
01787 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01788 }
01789 else {
01790 SMDS_MeshFace* F =
01791 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
01792 NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
01793 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01794 }
01795 }
01796 }
01797 }
01798
01799 else {
01800
01801 StdMeshers_Array2OfNode NodesBRD(1,nb,1,nnn-1);
01802
01803 for (j=0; j<nb; j++) {
01804 NodesBRD.SetValue(j+1,1,uv_eb[j].node);
01805 }
01806 for (i=1; i<nnn-1; i++) {
01807 NodesBRD.SetValue(1,i+1,uv_el[i].node);
01808 NodesBRD.SetValue(nb,i+1,uv_er[i].node);
01809 for (j=2; j<nb; j++) {
01810 double x = npb.Value(j);
01811 double y = (1-x) * npl.Value(i+1) + x * npr.Value(i+1);
01812 gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
01813 gp_Pnt P = S->Value(UV.X(),UV.Y());
01814 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01815 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(),UV.Y());
01816 NodesBRD.SetValue(j,i+1,N);
01817 }
01818 }
01819 for (j=1; j<nnn-1; j++) {
01820 for (i=1; i<nb; i++) {
01821 if (WisF) {
01822 SMDS_MeshFace* F =
01823 myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i+1,j),
01824 NodesBRD.Value(i+1,j+1), NodesBRD.Value(i,j+1));
01825 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01826 }
01827 else {
01828 SMDS_MeshFace* F =
01829 myHelper->AddFace(NodesBRD.Value(i,j), NodesBRD.Value(i,j+1),
01830 NodesBRD.Value(i+1,j+1), NodesBRD.Value(i+1,j));
01831 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01832 }
01833 }
01834 }
01835 int drl = abs(nr-nl);
01836
01837 StdMeshers_Array2OfNode NodesC(1,nb,1,drl+1+addv);
01838
01839 for (j=1; j<=nb; j++) {
01840 NodesC.SetValue(j,1,NodesBRD.Value(j,nnn-1));
01841 }
01842 if ((drl+addv) > 0) {
01843 int n1,n2;
01844 if (nr>nl) {
01845 n1 = 1;
01846 n2 = drl + 1;
01847 TColgp_SequenceOfXY UVtmp;
01848 double drparam = npr.Value(nr) - npr.Value(nnn-1);
01849 double dlparam = npl.Value(nnn) - npl.Value(nnn-1);
01850 double y0,y1;
01851 for (i=1; i<=drl; i++) {
01852
01853 NodesC.SetValue(nb,i+1,uv_er[nnn+i-2].node);
01854
01855 y1 = npr.Value(nnn+i-1);
01856 double dpar = (y1 - npr.Value(nnn-1))/drparam;
01857 y0 = npl.Value(nnn-1) + dpar*dlparam;
01858 double dy = y1 - y0;
01859 for (j=1; j<nb; j++) {
01860 double x = npt.Value(i+1) + npb.Value(j)*(1-npt.Value(i+1));
01861 double y = y0 + dy*x;
01862 gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
01863 gp_Pnt P = S->Value(UV.X(),UV.Y());
01864 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01865 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01866 NodesC.SetValue(j,i+1,N);
01867 }
01868 }
01869 double dy0 = (1-y0)/(addv+1);
01870 double dy1 = (1-y1)/(addv+1);
01871 for (i=1; i<=addv; i++) {
01872 double yy0 = y0 + dy0*i;
01873 double yy1 = y1 + dy1*i;
01874 double dyy = yy1 - yy0;
01875 for (j=1; j<=nb; j++) {
01876 double x = npt.Value(i+1+drl) +
01877 npb.Value(j) * (npt.Value(nt-i) - npt.Value(i+1+drl));
01878 double y = yy0 + dyy*x;
01879 gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
01880 gp_Pnt P = S->Value(UV.X(),UV.Y());
01881 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01882 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01883 NodesC.SetValue(j,i+drl+1,N);
01884 }
01885 }
01886 }
01887 else {
01888 n2 = 1;
01889 n1 = drl + 1;
01890 TColgp_SequenceOfXY UVtmp;
01891 double dlparam = npl.Value(nl) - npl.Value(nnn-1);
01892 double drparam = npr.Value(nnn) - npr.Value(nnn-1);
01893 double y0 = npl.Value(nnn-1);
01894 double y1 = npr.Value(nnn-1);
01895 for (i=1; i<=drl; i++) {
01896
01897 NodesC.SetValue(1,i+1,uv_el[nnn+i-2].node);
01898 y0 = npl.Value(nnn+i-1);
01899 double dpar = (y0 - npl.Value(nnn-1))/dlparam;
01900 y1 = npr.Value(nnn-1) + dpar*drparam;
01901 double dy = y1 - y0;
01902 for (j=2; j<=nb; j++) {
01903 double x = npb.Value(j)*npt.Value(nt-i);
01904 double y = y0 + dy*x;
01905 gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
01906 gp_Pnt P = S->Value(UV.X(),UV.Y());
01907 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01908 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01909 NodesC.SetValue(j,i+1,N);
01910 }
01911 }
01912 double dy0 = (1-y0)/(addv+1);
01913 double dy1 = (1-y1)/(addv+1);
01914 for (i=1; i<=addv; i++) {
01915 double yy0 = y0 + dy0*i;
01916 double yy1 = y1 + dy1*i;
01917 double dyy = yy1 - yy0;
01918 for (j=1; j<=nb; j++) {
01919 double x = npt.Value(i+1) +
01920 npb.Value(j) * (npt.Value(nt-i-drl) - npt.Value(i+1));
01921 double y = yy0 + dyy*x;
01922 gp_UV UV = CalcUV2(x, y, quad, a0, a1, a2, a3);
01923 gp_Pnt P = S->Value(UV.X(),UV.Y());
01924 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
01925 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
01926 NodesC.SetValue(j,i+drl+1,N);
01927 }
01928 }
01929 }
01930
01931 for (j=1; j<=drl+addv; j++) {
01932 for (i=1; i<nb; i++) {
01933 if (WisF) {
01934 SMDS_MeshFace* F =
01935 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
01936 NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
01937 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01938 }
01939 else {
01940 SMDS_MeshFace* F =
01941 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i,j+1),
01942 NodesC.Value(i+1,j+1), NodesC.Value(i+1,j));
01943 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01944 }
01945 }
01946 }
01947
01948 StdMeshers_Array2OfNode NodesLast(1,nt,1,2);
01949 for (i=1; i<=nt; i++) {
01950 NodesLast.SetValue(i,2,uv_et[i-1].node);
01951 }
01952 int nnn=0;
01953 for (i=n1; i<drl+addv+1; i++) {
01954 nnn++;
01955 NodesLast.SetValue(nnn,1,NodesC.Value(1,i));
01956 }
01957 for (i=1; i<=nb; i++) {
01958 nnn++;
01959 NodesLast.SetValue(nnn,1,NodesC.Value(i,drl+addv+1));
01960 }
01961 for (i=drl+addv; i>=n2; i--) {
01962 nnn++;
01963 NodesLast.SetValue(nnn,1,NodesC.Value(nb,i));
01964 }
01965 for (i=1; i<nt; i++) {
01966 if (WisF) {
01967 SMDS_MeshFace* F =
01968 myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i+1,1),
01969 NodesLast.Value(i+1,2), NodesLast.Value(i,2));
01970 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01971 }
01972 else {
01973 SMDS_MeshFace* F =
01974 myHelper->AddFace(NodesLast.Value(i,1), NodesLast.Value(i,2),
01975 NodesLast.Value(i+1,2), NodesLast.Value(i+1,2));
01976 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
01977 }
01978 }
01979 }
01980
01981 }
01982
01983 bool isOk = true;
01984 return isOk;
01985 }
01986
01987
01988
01992
01993
01994 bool StdMeshers_Quadrangle_2D::EvaluateQuadPref(SMESH_Mesh & aMesh,
01995 const TopoDS_Shape& aShape,
01996 std::vector<int>& aNbNodes,
01997 MapShapeNbElems& aResMap,
01998 bool IsQuadratic)
01999 {
02000
02001
02002
02003 bool OldVersion = false;
02004 if (myQuadType == QUAD_QUADRANGLE_PREF_REVERSED)
02005 OldVersion = true;
02006
02007 const TopoDS_Face& F = TopoDS::Face(aShape);
02008 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
02009
02010 int nb = aNbNodes[0];
02011 int nr = aNbNodes[1];
02012 int nt = aNbNodes[2];
02013 int nl = aNbNodes[3];
02014 int dh = abs(nb-nt);
02015 int dv = abs(nr-nl);
02016
02017 if (dh>=dv) {
02018 if (nt>nb) {
02019
02020 }
02021 else {
02022
02023 nb = aNbNodes[2];
02024 nr = aNbNodes[3];
02025 nt = aNbNodes[0];
02026 nl = aNbNodes[1];
02027 }
02028 }
02029 else {
02030 if (nr>nl) {
02031
02032 nb = aNbNodes[3];
02033 nr = aNbNodes[0];
02034 nt = aNbNodes[1];
02035 nl = aNbNodes[2];
02036 }
02037 else {
02038
02039 nb = aNbNodes[1];
02040 nr = aNbNodes[2];
02041 nt = aNbNodes[3];
02042 nl = aNbNodes[0];
02043 }
02044 }
02045
02046 dh = abs(nb-nt);
02047 dv = abs(nr-nl);
02048 int nbh = Max(nb,nt);
02049 int nbv = Max(nr,nl);
02050 int addh = 0;
02051 int addv = 0;
02052
02053 if (dh>dv) {
02054 addv = (dh-dv)/2;
02055 nbv = nbv + addv;
02056 }
02057 else {
02058 addh = (dv-dh)/2;
02059 nbh = nbh + addh;
02060 }
02061
02062 int dl,dr;
02063 if (OldVersion) {
02064
02065
02066 dr = nbv - nr;
02067
02068 dl = nbv - nl;
02069 }
02070
02071 int nnn = Min(nr,nl);
02072
02073 int nbNodes = 0;
02074 int nbFaces = 0;
02075 if (OldVersion) {
02076
02077 if (dl>0) {
02078 nbNodes += dl*(nl-1);
02079 nbFaces += dl*(nl-1);
02080 }
02081
02082 if (dr>0) {
02083 nbNodes += dr*(nr-1);
02084 nbFaces += dr*(nr-1);
02085 }
02086
02087 nbNodes += (nb-2)*(nnn-1) + (nbv-nnn-1)*(nb-2);
02088 nbFaces += (nb-1)*(nbv-1);
02089 }
02090 else {
02091 nbNodes += (nnn-2)*(nb-2);
02092 nbFaces += (nnn-2)*(nb-1);
02093 int drl = abs(nr-nl);
02094 nbNodes += drl*(nb-1) + addv*nb;
02095 nbFaces += (drl+addv)*(nb-1) + (nt-1);
02096 }
02097
02098 std::vector<int> aVec(SMDSEntity_Last);
02099 for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
02100 if (IsQuadratic) {
02101 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces;
02102 aVec[SMDSEntity_Node] = nbNodes + nbFaces*4;
02103 if (aNbNodes.size()==5) {
02104 aVec[SMDSEntity_Quad_Triangle] = aNbNodes[3] - 1;
02105 aVec[SMDSEntity_Quad_Quadrangle] = nbFaces - aNbNodes[3] + 1;
02106 }
02107 }
02108 else {
02109 aVec[SMDSEntity_Node] = nbNodes;
02110 aVec[SMDSEntity_Quadrangle] = nbFaces;
02111 if (aNbNodes.size()==5) {
02112 aVec[SMDSEntity_Triangle] = aNbNodes[3] - 1;
02113 aVec[SMDSEntity_Quadrangle] = nbFaces - aNbNodes[3] + 1;
02114 }
02115 }
02116 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
02117 aResMap.insert(std::make_pair(sm,aVec));
02118
02119 return true;
02120 }
02121
02122
02123
02127
02128 void StdMeshers_Quadrangle_2D::SplitQuad(SMESHDS_Mesh *theMeshDS,
02129 int theFaceID,
02130 const SMDS_MeshNode* theNode1,
02131 const SMDS_MeshNode* theNode2,
02132 const SMDS_MeshNode* theNode3,
02133 const SMDS_MeshNode* theNode4)
02134 {
02135 gp_Pnt a(theNode1->X(),theNode1->Y(),theNode1->Z());
02136 gp_Pnt b(theNode2->X(),theNode2->Y(),theNode2->Z());
02137 gp_Pnt c(theNode3->X(),theNode3->Y(),theNode3->Z());
02138 gp_Pnt d(theNode4->X(),theNode4->Y(),theNode4->Z());
02139 SMDS_MeshFace* face;
02140 if (a.Distance(c) > b.Distance(d)){
02141 face = myHelper->AddFace(theNode2, theNode4 , theNode1);
02142 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
02143 face = myHelper->AddFace(theNode2, theNode3, theNode4);
02144 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
02145
02146 }
02147 else{
02148 face = myHelper->AddFace(theNode1, theNode2 ,theNode3);
02149 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
02150 face = myHelper->AddFace(theNode1, theNode3, theNode4);
02151 if (face) theMeshDS->SetMeshElementOnShape(face, theFaceID);
02152 }
02153 }
02154
02155
02159
02160 bool StdMeshers_Quadrangle_2D::ComputeReduced (SMESH_Mesh & aMesh,
02161 const TopoDS_Shape& aShape,
02162 FaceQuadStruct* quad)
02163 {
02164 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
02165 const TopoDS_Face& F = TopoDS::Face(aShape);
02166 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
02167 int i,j,geomFaceID = meshDS->ShapeToIndex(F);
02168
02169 int nb = quad->side[0]->NbPoints();
02170 int nr = quad->side[1]->NbPoints();
02171 int nt = quad->side[2]->NbPoints();
02172 int nl = quad->side[3]->NbPoints();
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190 bool MultipleReduce = false;
02191 {
02192 int nb1 = nb;
02193 int nr1 = nr;
02194 int nt1 = nt;
02195
02196 if (nr == nl) {
02197 if (nb < nt) {
02198 nt1 = nb;
02199 nb1 = nt;
02200 }
02201 }
02202 else if (nb == nt) {
02203 nr1 = nb;
02204 if (nl < nr) {
02205 nt1 = nl;
02206 nb1 = nr;
02207 }
02208 else {
02209 nt1 = nr;
02210 nb1 = nl;
02211 }
02212 }
02213 else {
02214 return false;
02215 }
02216
02217
02218 int nrows = nr1 - 1;
02219 int ncol_top = nt1 - 1;
02220 int ncol_bot = nb1 - 1;
02221
02222 int nrows_tree31 = int( log( ncol_bot / ncol_top ) / log( 3 ));
02223 if ( nrows < nrows_tree31 )
02224 MultipleReduce = true;
02225 }
02226
02227 if (MultipleReduce) {
02228
02229 int dh = abs(nb-nt);
02230 int dv = abs(nr-nl);
02231
02232 if (dh >= dv) {
02233 if (nt > nb) {
02234
02235 ShiftQuad(quad,0,true);
02236 }
02237 else {
02238
02239 ShiftQuad(quad,2,true);
02240 }
02241 }
02242 else {
02243 if (nr > nl) {
02244
02245 ShiftQuad(quad,1,true);
02246 }
02247 else {
02248
02249 ShiftQuad(quad,3,true);
02250 }
02251 }
02252
02253 nb = quad->side[0]->NbPoints();
02254 nr = quad->side[1]->NbPoints();
02255 nt = quad->side[2]->NbPoints();
02256 nl = quad->side[3]->NbPoints();
02257 dh = abs(nb-nt);
02258 dv = abs(nr-nl);
02259 int nbh = Max(nb,nt);
02260 int nbv = Max(nr,nl);
02261 int addh = 0;
02262 int addv = 0;
02263
02264 if (dh>dv) {
02265 addv = (dh-dv)/2;
02266 nbv = nbv + addv;
02267 }
02268 else {
02269 addh = (dv-dh)/2;
02270 nbh = nbh + addh;
02271 }
02272
02273 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
02274 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
02275 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
02276 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
02277
02278 if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
02279 return error(COMPERR_BAD_INPUT_MESH);
02280
02281 if ( myNeedSmooth )
02282 UpdateDegenUV( quad );
02283
02284
02285 TColStd_SequenceOfReal npb, npr, npt, npl;
02286 for (j = 0; j < nb; j++) {
02287 npb.Append(uv_eb[j].normParam);
02288 }
02289 for (i = 0; i < nr; i++) {
02290 npr.Append(uv_er[i].normParam);
02291 }
02292 for (j = 0; j < nt; j++) {
02293 npt.Append(uv_et[j].normParam);
02294 }
02295 for (i = 0; i < nl; i++) {
02296 npl.Append(uv_el[i].normParam);
02297 }
02298
02299 int dl,dr;
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315 dr = nbv - nr;
02316 double dpr = (npr.Value(2) - npr.Value(1))/(dr+1);
02317 for (i=1; i<=dr; i++) {
02318 npr.InsertAfter(1,npr.Value(2)-dpr);
02319 }
02320
02321 dl = nbv - nl;
02322 dpr = (npl.Value(2) - npl.Value(1))/(dl+1);
02323 for (i=1; i<=dl; i++) {
02324 npl.InsertAfter(1,npl.Value(2)-dpr);
02325 }
02326
02327 gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
02328 gp_XY a1 (uv_eb.back().u, uv_eb.back().v);
02329 gp_XY a2 (uv_et.back().u, uv_et.back().v);
02330 gp_XY a3 (uv_et.front().u, uv_et.front().v);
02331
02332 int nnn = Min(nr,nl);
02333
02334
02335
02336 TColgp_SequenceOfXY UVL;
02337 TColgp_SequenceOfXY UVR;
02338
02339
02340
02341 StdMeshers_Array2OfNode NodesL(1,dl+1,1,nl);
02342
02343 for (j=1; j<=nl; j++)
02344 NodesL.SetValue(1,j,uv_el[j-1].node);
02345 if (dl>0) {
02346
02347 for (i=1; i<=dl; i++)
02348 NodesL.SetValue(i+1,nl,uv_et[i].node);
02349
02350 TColgp_SequenceOfXY UVtmp;
02351 for (i=1; i<=dl; i++) {
02352 double x0 = npt.Value(i+1);
02353 double x1 = x0;
02354
02355 double y0 = npl.Value(i+1);
02356 double y1 = npr.Value(i+1);
02357 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
02358 gp_Pnt P = S->Value(UV.X(),UV.Y());
02359 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
02360 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
02361 NodesL.SetValue(i+1,1,N);
02362 if (UVL.Length()<nbv-nnn-1) UVL.Append(UV);
02363
02364 for (j=2; j<nl; j++) {
02365 double y0 = npl.Value(dl+j);
02366 double y1 = npr.Value(dl+j);
02367 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
02368 gp_Pnt P = S->Value(UV.X(),UV.Y());
02369 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
02370 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
02371 NodesL.SetValue(i+1,j,N);
02372 if (i==dl) UVtmp.Append(UV);
02373 }
02374 }
02375 for (i=1; i<=UVtmp.Length() && UVL.Length()<nbv-nnn-1; i++) {
02376 UVL.Append(UVtmp.Value(i));
02377 }
02378
02379 for (i=1; i<=dl; i++) {
02380 for (j=1; j<nl; j++) {
02381 SMDS_MeshFace* F =
02382 myHelper->AddFace(NodesL.Value(i,j), NodesL.Value(i+1,j),
02383 NodesL.Value(i+1,j+1), NodesL.Value(i,j+1));
02384 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
02385 }
02386 }
02387 }
02388 else {
02389
02390 for (i=1; i<npl.Length() && UVL.Length()<nbv-nnn-1; i++) {
02391 UVL.Append(gp_UV (uv_el[i].u, uv_el[i].v));
02392 }
02393 }
02394
02395
02396 StdMeshers_Array2OfNode NodesR(1,dr+1,1,nr);
02397
02398 for (j=1; j<=nr; j++)
02399 NodesR.SetValue(1,j,uv_er[nr-j].node);
02400 if (dr>0) {
02401
02402 for (i=1; i<=dr; i++)
02403 NodesR.SetValue(i+1,1,uv_et[nt-1-i].node);
02404
02405 TColgp_SequenceOfXY UVtmp;
02406 for (i=1; i<=dr; i++) {
02407 double x0 = npt.Value(nt-i);
02408 double x1 = x0;
02409
02410 double y0 = npl.Value(i+1);
02411 double y1 = npr.Value(i+1);
02412 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
02413 gp_Pnt P = S->Value(UV.X(),UV.Y());
02414 SMDS_MeshNode * N = meshDS->AddNode(P.X(), P.Y(), P.Z());
02415 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
02416 NodesR.SetValue(i+1,nr,N);
02417 if (UVR.Length()<nbv-nnn-1) UVR.Append(UV);
02418
02419 for (j=2; j<nr; j++) {
02420 double y0 = npl.Value(nbv-j+1);
02421 double y1 = npr.Value(nbv-j+1);
02422 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
02423 gp_Pnt P = S->Value(UV.X(),UV.Y());
02424 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
02425 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
02426 NodesR.SetValue(i+1,j,N);
02427 if (i==dr) UVtmp.Prepend(UV);
02428 }
02429 }
02430 for (i=1; i<=UVtmp.Length() && UVR.Length()<nbv-nnn-1; i++) {
02431 UVR.Append(UVtmp.Value(i));
02432 }
02433
02434 for (i=1; i<=dr; i++) {
02435 for (j=1; j<nr; j++) {
02436 SMDS_MeshFace* F =
02437 myHelper->AddFace(NodesR.Value(i,j), NodesR.Value(i+1,j),
02438 NodesR.Value(i+1,j+1), NodesR.Value(i,j+1));
02439 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
02440 }
02441 }
02442 }
02443 else {
02444
02445 for (i=1; i<npr.Length() && UVR.Length()<nbv-nnn-1; i++) {
02446 UVR.Append(gp_UV(uv_er[i].u, uv_er[i].v));
02447 }
02448 }
02449
02450
02451 StdMeshers_Array2OfNode NodesC(1,nb,1,nbv);
02452
02453 for (i=1; i<=dl+1; i++)
02454 NodesC.SetValue(1,i,NodesL(i,1));
02455 for (i=2; i<=nl; i++)
02456 NodesC.SetValue(1,dl+i,NodesL(dl+1,i));
02457
02458 for (i=1; i<=dr+1; i++)
02459 NodesC.SetValue(nb,i,NodesR(i,nr));
02460 for (i=1; i<nr; i++)
02461 NodesC.SetValue(nb,dr+i+1,NodesR(dr+1,nr-i));
02462
02463 for (i=dl+2; i<nbh-dr; i++)
02464 NodesC.SetValue(i-dl,nbv,uv_et[i-1].node);
02465
02466 for (i=2; i<nb; i++)
02467 NodesC.SetValue(i,1,uv_eb[i-1].node);
02468
02469
02470
02471 for (i=2; i<nb; i++) {
02472 double x0 = npt.Value(dl+i);
02473 double x1 = x0;
02474 for (j=1; j<nnn; j++) {
02475 double y0 = npl.Value(nbv-nnn+j);
02476 double y1 = npr.Value(nbv-nnn+j);
02477 gp_UV UV = CalcUV(x0, x1, y0, y1, quad, a0, a1, a2, a3);
02478 gp_Pnt P = S->Value(UV.X(),UV.Y());
02479 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
02480 meshDS->SetNodeOnFace(N, geomFaceID, UV.X(), UV.Y());
02481 NodesC.SetValue(i,nbv-nnn+j,N);
02482 }
02483 }
02484
02485 for (i=1; i<nbv-nnn; i++) {
02486 double du = UVR.Value(i).X() - UVL.Value(i).X();
02487 double dv = UVR.Value(i).Y() - UVL.Value(i).Y();
02488 for (j=2; j<nb; j++) {
02489 double u = UVL.Value(i).X() + du*npb.Value(j);
02490 double v = UVL.Value(i).Y() + dv*npb.Value(j);
02491 gp_Pnt P = S->Value(u,v);
02492 SMDS_MeshNode* N = meshDS->AddNode(P.X(), P.Y(), P.Z());
02493 meshDS->SetNodeOnFace(N, geomFaceID, u, v);
02494 NodesC.SetValue(j,i+1,N);
02495 }
02496 }
02497
02498 for (i=1; i<nb; i++) {
02499 for (j=1; j<nbv; j++) {
02500 SMDS_MeshFace* F =
02501 myHelper->AddFace(NodesC.Value(i,j), NodesC.Value(i+1,j),
02502 NodesC.Value(i+1,j+1), NodesC.Value(i,j+1));
02503 if (F) meshDS->SetMeshElementOnShape(F, geomFaceID);
02504 }
02505 }
02506
02507 }
02508 else {
02509
02510 if (nr == nl) {
02511 if (nt < nb) {
02512
02513
02514 }
02515 else {
02516
02517 ShiftQuad(quad,2,true);
02518 }
02519 }
02520 else {
02521 if (nl > nr) {
02522
02523 ShiftQuad(quad,1,true);
02524 }
02525 else {
02526
02527 ShiftQuad(quad,3,true);
02528 }
02529 }
02530
02531 nb = quad->side[0]->NbPoints();
02532 nr = quad->side[1]->NbPoints();
02533 nt = quad->side[2]->NbPoints();
02534 nl = quad->side[3]->NbPoints();
02535
02536
02537 int nrows = nr - 1;
02538 int ncol_top = nt - 1;
02539 int ncol_bot = nb - 1;
02540 int npair_top = ncol_top / 2;
02541
02542 int max_lin = ncol_top + npair_top * 2 * nrows;
02543
02544 int max_lin31 = ncol_top + ncol_top * 2 * nrows;
02545
02546 int max_tree42 = 0;
02547
02548 int nrows_tree42 = int( log2( ncol_bot / ncol_top ));
02549 if (ncol_top > npair_top * 2 && nrows_tree42 < nrows) {
02550 max_tree42 = npair_top * pow(2.0, nrows + 1);
02551 int delta = ncol_bot - int( max_tree42 );
02552 for (int irow = 1; irow < nrows; irow++) {
02553 int nfour = delta / 4;
02554 delta -= nfour * 2;
02555 }
02556 if (delta <= (ncol_top - npair_top * 2))
02557 max_tree42 = ncol_bot;
02558 }
02559
02560
02561 bool is_lin_31 = false;
02562 bool is_lin_42 = false;
02563 bool is_tree_31 = false;
02564 bool is_tree_42 = false;
02565 if (ncol_bot > max_lin) {
02566 if (ncol_bot <= max_lin31) {
02567 is_lin_31 = true;
02568 max_lin = max_lin31;
02569 }
02570 }
02571 else {
02572
02573 if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
02574 is_lin_31 = true;
02575 max_lin = max_lin31;
02576 }
02577 else {
02578 is_lin_42 = true;
02579 }
02580 }
02581 if (ncol_bot > max_lin) {
02582 is_tree_31 = (ncol_bot > max_tree42);
02583 if (ncol_bot <= max_tree42) {
02584 if ((ncol_bot/3)*3 == ncol_bot || (ncol_bot/2)*2 != ncol_bot) {
02585 is_tree_31 = true;
02586 }
02587 else {
02588 is_tree_42 = true;
02589 }
02590 }
02591 }
02592
02593 const vector<UVPtStruct>& uv_eb = quad->side[0]->GetUVPtStruct(true,0);
02594 const vector<UVPtStruct>& uv_er = quad->side[1]->GetUVPtStruct(false,1);
02595 const vector<UVPtStruct>& uv_et = quad->side[2]->GetUVPtStruct(true,1);
02596 const vector<UVPtStruct>& uv_el = quad->side[3]->GetUVPtStruct(false,0);
02597
02598 if (uv_eb.size() != nb || uv_er.size() != nr || uv_et.size() != nt || uv_el.size() != nl)
02599 return error(COMPERR_BAD_INPUT_MESH);
02600
02601
02602 TColStd_SequenceOfReal npb, npr, npt, npl;
02603 for (j = 0; j < nb; j++) {
02604 npb.Append(uv_eb[j].normParam);
02605 }
02606 for (i = 0; i < nr; i++) {
02607 npr.Append(uv_er[i].normParam);
02608 }
02609 for (j = 0; j < nt; j++) {
02610 npt.Append(uv_et[j].normParam);
02611 }
02612 for (i = 0; i < nl; i++) {
02613 npl.Append(uv_el[i].normParam);
02614 }
02615
02616
02617 if (!SetNormalizedGrid(aMesh, aShape, quad))
02618 return false;
02619
02620
02621 gp_XY a0 (uv_eb.front().u, uv_eb.front().v);
02622 gp_XY a1 (uv_eb.back().u, uv_eb.back().v);
02623 gp_XY a2 (uv_et.back().u, uv_et.back().v);
02624 gp_XY a3 (uv_et.front().u, uv_et.front().v);
02625
02626
02627 TColStd_SequenceOfInteger curr_base, next_base;
02628 TColStd_SequenceOfReal curr_par_u, curr_par_v;
02629 TColStd_SequenceOfReal next_par_u, next_par_v;
02630 StdMeshers_Array2OfNode NodesBRD (1,nb, 1,nr);
02631 for (j = 1; j <= nb; j++) {
02632 NodesBRD.SetValue(j, 1, uv_eb[j - 1].node);
02633 curr_base.Append(j);
02634 next_base.Append(-1);
02635 curr_par_u.Append(uv_eb[j-1].u);
02636 curr_par_v.Append(uv_eb[j-1].v);
02637 next_par_u.Append(0.);
02638 next_par_v.Append(0.);
02639 }
02640 for (j = 1; j <= nt; j++) {
02641 NodesBRD.SetValue(j, nr, uv_et[j - 1].node);
02642 }
02643
02644 int curr_base_len = nb;
02645 int next_base_len = 0;
02646
02647 if (is_tree_42) {
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669 for (i = 1; i < nr; i++) {
02670
02671 NodesBRD.SetValue(1, i+1, uv_el[i].node);
02672 next_base.SetValue(++next_base_len, 1);
02673
02674 NodesBRD.SetValue(nb, i+1, uv_er[i].node);
02675
02676 next_par_u.SetValue(next_base_len, uv_el[i].u);
02677 next_par_v.SetValue(next_base_len, uv_el[i].v);
02678
02679
02680 int delta = curr_base_len - nt;
02681
02682
02683
02684
02685
02686 int nb_four = (curr_base_len - 1) / 4;
02687 int nb_next = nb_four*2 + (curr_base_len - nb_four*4);
02688 if (nb_next < nt) nb_next = nt;
02689
02690 for (j = 1; j + 4 <= curr_base_len && delta > 0; j += 4, delta -= 2) {
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703 double u,v;
02704
02705
02706 const SMDS_MeshNode* Na;
02707 next_base_len++;
02708 next_base.SetValue(next_base_len, curr_base.Value(j + 2));
02709 if (i + 1 == nr) {
02710 Na = uv_et[next_base_len - 1].node;
02711 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
02712 u = uv_et[next_base_len - 1].u;
02713 v = uv_et[next_base_len - 1].v;
02714 }
02715 else {
02716
02717
02718
02719 {
02720 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
02721 int nearest_node_j = (int)rel;
02722 rel -= nearest_node_j;
02723 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
02724 double u1 = quad->uv_grid[ij].u;
02725 double v1 = quad->uv_grid[ij].v;
02726 double u2 = quad->uv_grid[ij + 1].u;
02727 double v2 = quad->uv_grid[ij + 1].v;
02728 double duj = (u2 - u1) * rel;
02729 double dvj = (v2 - v1) * rel;
02730 u = u1 + duj;
02731 v = v1 + dvj;
02732 }
02733
02734
02735 gp_Pnt P = S->Value(u,v);
02736 SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
02737 meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
02738 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
02739 Na = Na1;
02740 }
02741 next_par_u.SetValue(next_base_len, u);
02742 next_par_v.SetValue(next_base_len, v);
02743
02744
02745 const SMDS_MeshNode* Nb;
02746 next_base_len++;
02747 next_base.SetValue(next_base_len, curr_base.Value(j + 4));
02748 if (i + 1 == nr) {
02749 Nb = uv_et[next_base_len - 1].node;
02750 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
02751 u = uv_et[next_base_len - 1].u;
02752 v = uv_et[next_base_len - 1].v;
02753 }
02754 else if (j + 4 == curr_base_len) {
02755 Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
02756 u = uv_er[i].u;
02757 v = uv_er[i].v;
02758 }
02759 else {
02760
02761
02762
02763 {
02764 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
02765 int nearest_node_j = (int)rel;
02766 rel -= nearest_node_j;
02767 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
02768 double u1 = quad->uv_grid[ij].u;
02769 double v1 = quad->uv_grid[ij].v;
02770 double u2 = quad->uv_grid[ij + 1].u;
02771 double v2 = quad->uv_grid[ij + 1].v;
02772 double duj = (u2 - u1) * rel;
02773 double dvj = (v2 - v1) * rel;
02774 u = u1 + duj;
02775 v = v1 + dvj;
02776 }
02777
02778
02779 gp_Pnt P = S->Value(u,v);
02780 SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
02781 meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
02782 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
02783 Nb = Nb1;
02784 }
02785 next_par_u.SetValue(next_base_len, u);
02786 next_par_v.SetValue(next_base_len, v);
02787
02788
02789 u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
02790 v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
02791 gp_Pnt P = S->Value(u,v);
02792 SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
02793 meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
02794
02795
02796 u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
02797 v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
02798 P = S->Value(u,v);
02799 SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
02800 meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
02801
02802
02803 u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
02804 v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
02805 P = S->Value(u,v);
02806 SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
02807 meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
02808
02809
02810 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
02811 NodesBRD.Value(curr_base.Value(j + 1), i),
02812 Nc,
02813 NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
02814 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
02815
02816 SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
02817 NodesBRD.Value(curr_base.Value(j + 2), i),
02818 Nd, Nc);
02819 if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
02820
02821 SMDS_MeshFace* F3 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
02822 NodesBRD.Value(curr_base.Value(j + 3), i),
02823 Ne, Nd);
02824 if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
02825
02826 SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
02827 NodesBRD.Value(curr_base.Value(j + 4), i),
02828 Nb, Ne);
02829 if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
02830
02831 SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Nd, Na,
02832 NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
02833 if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
02834
02835 SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na);
02836 if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
02837 }
02838
02839
02840 for (; j < curr_base_len; j++) {
02841
02842 const SMDS_MeshNode* Nf;
02843 double u,v;
02844 next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
02845 if (i + 1 == nr) {
02846 Nf = uv_et[next_base_len - 1].node;
02847 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
02848 u = uv_et[next_base_len - 1].u;
02849 v = uv_et[next_base_len - 1].v;
02850 }
02851 else if (j + 1 == curr_base_len) {
02852 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
02853 u = uv_er[i].u;
02854 v = uv_er[i].v;
02855 }
02856 else {
02857
02858
02859
02860 {
02861 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
02862 int nearest_node_j = (int)rel;
02863 rel -= nearest_node_j;
02864 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
02865 double u1 = quad->uv_grid[ij].u;
02866 double v1 = quad->uv_grid[ij].v;
02867 double u2 = quad->uv_grid[ij + 1].u;
02868 double v2 = quad->uv_grid[ij + 1].v;
02869 double duj = (u2 - u1) * rel;
02870 double dvj = (v2 - v1) * rel;
02871 u = u1 + duj;
02872 v = v1 + dvj;
02873 }
02874
02875
02876 gp_Pnt P = S->Value(u,v);
02877 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
02878 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
02879 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
02880 Nf = Nf1;
02881 }
02882 next_par_u.SetValue(next_base_len, u);
02883 next_par_v.SetValue(next_base_len, v);
02884 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
02885 NodesBRD.Value(curr_base.Value(j + 1), i),
02886 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
02887 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
02888 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
02889 }
02890
02891 curr_base_len = next_base_len;
02892 curr_base = next_base;
02893 curr_par_u = next_par_u;
02894 curr_par_v = next_par_v;
02895 next_base_len = 0;
02896 }
02897 }
02898 else if (is_tree_31) {
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916 for (i = 1; i < nr; i++) {
02917
02918 NodesBRD.SetValue(1, i+1, uv_el[i].node);
02919 next_base.SetValue(++next_base_len, 1);
02920
02921 NodesBRD.SetValue(nb, i+1, uv_er[i].node);
02922
02923 next_par_u.SetValue(next_base_len, uv_el[i].u);
02924 next_par_v.SetValue(next_base_len, uv_el[i].v);
02925
02926
02927 int delta = curr_base_len - nt;
02928
02929
02930 int nb_three = (curr_base_len - 1) / 3;
02931 int nb_next = nb_three + (curr_base_len - nb_three*3);
02932 if (nb_next < nt) nb_next = nt;
02933
02934 for (j = 1; j + 3 <= curr_base_len && delta > 0; j += 3, delta -= 2) {
02935
02936
02937
02938
02939
02940
02941
02942
02943
02944
02945
02946
02947 double u,v;
02948
02949
02950 const SMDS_MeshNode* Nb;
02951 next_base_len++;
02952 next_base.SetValue(next_base_len, curr_base.Value(j + 3));
02953 if (i + 1 == nr) {
02954 Nb = uv_et[next_base_len - 1].node;
02955 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
02956 u = uv_et[next_base_len - 1].u;
02957 v = uv_et[next_base_len - 1].v;
02958 }
02959 else if (j + 3 == curr_base_len) {
02960 Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
02961 u = uv_er[i].u;
02962 v = uv_er[i].v;
02963 }
02964 else {
02965 {
02966 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
02967 int nearest_node_j = (int)rel;
02968 rel -= nearest_node_j;
02969 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
02970 double u1 = quad->uv_grid[ij].u;
02971 double v1 = quad->uv_grid[ij].v;
02972 double u2 = quad->uv_grid[ij + 1].u;
02973 double v2 = quad->uv_grid[ij + 1].v;
02974 double duj = (u2 - u1) * rel;
02975 double dvj = (v2 - v1) * rel;
02976 u = u1 + duj;
02977 v = v1 + dvj;
02978 }
02979 gp_Pnt P = S->Value(u,v);
02980 SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
02981 meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
02982 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
02983 Nb = Nb1;
02984 }
02985 next_par_u.SetValue(next_base_len, u);
02986 next_par_v.SetValue(next_base_len, v);
02987
02988
02989 double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
02990 double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
02991 double u3 = (u2 - u1) / 3.0;
02992
02993 double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
02994 double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
02995 double v3 = (v2 - v1) / 3.0;
02996
02997
02998 u = u1 + u3;
02999 v = v1 + v3;
03000 gp_Pnt P = S->Value(u,v);
03001 SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
03002 meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
03003
03004
03005 u = u1 + u3 + u3;
03006 v = v1 + v3 + v3;
03007 P = S->Value(u,v);
03008 SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
03009 meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
03010
03011
03012 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
03013 NodesBRD.Value(curr_base.Value(j + 1), i),
03014 Nc,
03015 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03016 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03017
03018 SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
03019 NodesBRD.Value(curr_base.Value(j + 2), i),
03020 Ne, Nc);
03021 if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
03022
03023 SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
03024 NodesBRD.Value(curr_base.Value(j + 3), i),
03025 Nb, Ne);
03026 if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
03027
03028 SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Ne, Nb,
03029 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03030 if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
03031 }
03032
03033
03034 for (; j < curr_base_len; j++) {
03035
03036 const SMDS_MeshNode* Nf;
03037 double u,v;
03038 next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
03039 if (i + 1 == nr) {
03040 Nf = uv_et[next_base_len - 1].node;
03041 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03042 u = uv_et[next_base_len - 1].u;
03043 v = uv_et[next_base_len - 1].v;
03044 }
03045 else if (j + 1 == curr_base_len) {
03046 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03047 u = uv_er[i].u;
03048 v = uv_er[i].v;
03049 }
03050 else {
03051 {
03052 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03053 int nearest_node_j = (int)rel;
03054 rel -= nearest_node_j;
03055 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03056 double u1 = quad->uv_grid[ij].u;
03057 double v1 = quad->uv_grid[ij].v;
03058 double u2 = quad->uv_grid[ij + 1].u;
03059 double v2 = quad->uv_grid[ij + 1].v;
03060 double duj = (u2 - u1) * rel;
03061 double dvj = (v2 - v1) * rel;
03062 u = u1 + duj;
03063 v = v1 + dvj;
03064 }
03065 gp_Pnt P = S->Value(u,v);
03066 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03067 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03068 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03069 Nf = Nf1;
03070 }
03071 next_par_u.SetValue(next_base_len, u);
03072 next_par_v.SetValue(next_base_len, v);
03073 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
03074 NodesBRD.Value(curr_base.Value(j + 1), i),
03075 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03076 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03077 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03078 }
03079
03080 curr_base_len = next_base_len;
03081 curr_base = next_base;
03082 curr_par_u = next_par_u;
03083 curr_par_v = next_par_v;
03084 next_base_len = 0;
03085 }
03086 }
03087 else if (is_lin_42) {
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103
03104
03105
03106
03107
03108
03109
03110
03111
03112
03113
03114
03115
03116 int delta_all = nb - nt;
03117 int delta_one_col = (nr - 1) * 2;
03118 int nb_col = delta_all / delta_one_col;
03119 int remainder = delta_all - nb_col * delta_one_col;
03120 if (remainder > 0) {
03121 nb_col++;
03122 }
03123 int free_left = ((nt - 1) - nb_col * 2) / 2;
03124 free_left += nr - 2;
03125 int free_middle = (nr - 2) * 2;
03126 if (remainder > 0 && nb_col == 1) {
03127 int nb_rows_short_col = remainder / 2;
03128 int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
03129 free_left -= nb_rows_thrown;
03130 }
03131
03132
03133
03134
03135
03136
03137
03138
03139
03140 for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) {
03141
03142 NodesBRD.SetValue(1, i+1, uv_el[i].node);
03143 next_base.SetValue(++next_base_len, 1);
03144
03145 NodesBRD.SetValue(nb, i+1, uv_er[i].node);
03146
03147
03148 next_par_u.SetValue(next_base_len, uv_el[i].u);
03149 next_par_v.SetValue(next_base_len, uv_el[i].v);
03150
03151
03152 int nb_next = curr_base_len - nb_col * 2;
03153 if (remainder > 0 && i > remainder / 2)
03154
03155 nb_next += 2;
03156 if (nb_next < nt) nb_next = nt;
03157
03158
03159 for (j = 1; j <= free_left; j++) {
03160
03161 const SMDS_MeshNode* Nf;
03162 double u,v;
03163 next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
03164 if (i + 1 == nr) {
03165 Nf = uv_et[next_base_len - 1].node;
03166 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03167 u = uv_et[next_base_len - 1].u;
03168 v = uv_et[next_base_len - 1].v;
03169 }
03170 else {
03171 {
03172 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03173 int nearest_node_j = (int)rel;
03174 rel -= nearest_node_j;
03175 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03176 double u1 = quad->uv_grid[ij].u;
03177 double v1 = quad->uv_grid[ij].v;
03178 double u2 = quad->uv_grid[ij + 1].u;
03179 double v2 = quad->uv_grid[ij + 1].v;
03180 double duj = (u2 - u1) * rel;
03181 double dvj = (v2 - v1) * rel;
03182 u = u1 + duj;
03183 v = v1 + dvj;
03184 }
03185 gp_Pnt P = S->Value(u,v);
03186 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03187 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03188 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03189 Nf = Nf1;
03190 }
03191 next_par_u.SetValue(next_base_len, u);
03192 next_par_v.SetValue(next_base_len, v);
03193 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
03194 NodesBRD.Value(curr_base.Value(j + 1), i),
03195 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03196 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03197 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03198 }
03199
03200 for (int icol = 1; icol <= nb_col; icol++) {
03201
03202 if (remainder > 0 && icol == nb_col && i > remainder / 2)
03203
03204 break;
03205
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218 double u,v;
03219
03220
03221 const SMDS_MeshNode* Na;
03222 next_base_len++;
03223 next_base.SetValue(next_base_len, curr_base.Value(j + 2));
03224 if (i + 1 == nr) {
03225 Na = uv_et[next_base_len - 1].node;
03226 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na);
03227 u = uv_et[next_base_len - 1].u;
03228 v = uv_et[next_base_len - 1].v;
03229 }
03230 else {
03231 {
03232 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03233 int nearest_node_j = (int)rel;
03234 rel -= nearest_node_j;
03235 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03236 double u1 = quad->uv_grid[ij].u;
03237 double v1 = quad->uv_grid[ij].v;
03238 double u2 = quad->uv_grid[ij + 1].u;
03239 double v2 = quad->uv_grid[ij + 1].v;
03240 double duj = (u2 - u1) * rel;
03241 double dvj = (v2 - v1) * rel;
03242 u = u1 + duj;
03243 v = v1 + dvj;
03244 }
03245 gp_Pnt P = S->Value(u,v);
03246 SMDS_MeshNode* Na1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03247 meshDS->SetNodeOnFace(Na1, geomFaceID, u, v);
03248 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Na1);
03249 Na = Na1;
03250 }
03251 next_par_u.SetValue(next_base_len, u);
03252 next_par_v.SetValue(next_base_len, v);
03253
03254
03255 const SMDS_MeshNode* Nb;
03256 next_base_len++;
03257 next_base.SetValue(next_base_len, curr_base.Value(j + 4));
03258 if (i + 1 == nr) {
03259 Nb = uv_et[next_base_len - 1].node;
03260 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
03261 u = uv_et[next_base_len - 1].u;
03262 v = uv_et[next_base_len - 1].v;
03263 }
03264 else if (j + 4 == curr_base_len) {
03265 Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03266 u = uv_er[i].u;
03267 v = uv_er[i].v;
03268 }
03269 else {
03270 {
03271 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03272 int nearest_node_j = (int)rel;
03273 rel -= nearest_node_j;
03274 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03275 double u1 = quad->uv_grid[ij].u;
03276 double v1 = quad->uv_grid[ij].v;
03277 double u2 = quad->uv_grid[ij + 1].u;
03278 double v2 = quad->uv_grid[ij + 1].v;
03279 double duj = (u2 - u1) * rel;
03280 double dvj = (v2 - v1) * rel;
03281 u = u1 + duj;
03282 v = v1 + dvj;
03283 }
03284 gp_Pnt P = S->Value(u,v);
03285 SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03286 meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
03287 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
03288 Nb = Nb1;
03289 }
03290 next_par_u.SetValue(next_base_len, u);
03291 next_par_v.SetValue(next_base_len, v);
03292
03293
03294 u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 2)) / 2.0;
03295 v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 2)) / 2.0;
03296 gp_Pnt P = S->Value(u,v);
03297 SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
03298 meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
03299
03300
03301 u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len - 1)) / 2.0;
03302 v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len - 1)) / 2.0;
03303 P = S->Value(u,v);
03304 SMDS_MeshNode* Nd = meshDS->AddNode(P.X(), P.Y(), P.Z());
03305 meshDS->SetNodeOnFace(Nd, geomFaceID, u, v);
03306
03307
03308 u = (curr_par_u.Value(j + 2) + next_par_u.Value(next_base_len)) / 2.0;
03309 v = (curr_par_v.Value(j + 2) + next_par_v.Value(next_base_len)) / 2.0;
03310 P = S->Value(u,v);
03311 SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
03312 meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
03313
03314
03315 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
03316 NodesBRD.Value(curr_base.Value(j + 1), i),
03317 Nc,
03318 NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
03319 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03320
03321 SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
03322 NodesBRD.Value(curr_base.Value(j + 2), i),
03323 Nd, Nc);
03324 if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
03325
03326 SMDS_MeshFace* F3 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
03327 NodesBRD.Value(curr_base.Value(j + 3), i),
03328 Ne, Nd);
03329 if (F3) meshDS->SetMeshElementOnShape(F3, geomFaceID);
03330
03331 SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 3), i),
03332 NodesBRD.Value(curr_base.Value(j + 4), i),
03333 Nb, Ne);
03334 if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
03335
03336 SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Nd, Na,
03337 NodesBRD.Value(next_base.Value(next_base_len - 2), i + 1));
03338 if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
03339
03340 SMDS_MeshFace* F6 = myHelper->AddFace(Nd, Ne, Nb, Na);
03341 if (F6) meshDS->SetMeshElementOnShape(F6, geomFaceID);
03342
03343 j += 4;
03344
03345
03346 if (icol < nb_col) {
03347 if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
03348
03349 break;
03350
03351 int free_add = free_middle;
03352 if (remainder > 0 && icol == nb_col - 1)
03353
03354 free_add -= (nr - 1) - (remainder / 2);
03355
03356 for (int imiddle = 1; imiddle <= free_add; imiddle++) {
03357
03358 const SMDS_MeshNode* Nf;
03359 double u,v;
03360 next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
03361 if (i + 1 == nr) {
03362 Nf = uv_et[next_base_len - 1].node;
03363 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03364 u = uv_et[next_base_len - 1].u;
03365 v = uv_et[next_base_len - 1].v;
03366 }
03367 else if (j + imiddle == curr_base_len) {
03368 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03369 u = uv_er[i].u;
03370 v = uv_er[i].v;
03371 }
03372 else {
03373 {
03374 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03375 int nearest_node_j = (int)rel;
03376 rel -= nearest_node_j;
03377 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03378 double u1 = quad->uv_grid[ij].u;
03379 double v1 = quad->uv_grid[ij].v;
03380 double u2 = quad->uv_grid[ij + 1].u;
03381 double v2 = quad->uv_grid[ij + 1].v;
03382 double duj = (u2 - u1) * rel;
03383 double dvj = (v2 - v1) * rel;
03384 u = u1 + duj;
03385 v = v1 + dvj;
03386 }
03387 gp_Pnt P = S->Value(u,v);
03388 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03389 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03390 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03391 Nf = Nf1;
03392 }
03393 next_par_u.SetValue(next_base_len, u);
03394 next_par_v.SetValue(next_base_len, v);
03395 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
03396 NodesBRD.Value(curr_base.Value(j + imiddle), i),
03397 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03398 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03399 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03400 }
03401 j += free_add;
03402 }
03403 }
03404
03405
03406 for (; j < curr_base_len; j++) {
03407
03408 const SMDS_MeshNode* Nf;
03409 double u,v;
03410 next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
03411 if (i + 1 == nr) {
03412 Nf = uv_et[next_base_len - 1].node;
03413 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03414 u = uv_et[next_base_len - 1].u;
03415 v = uv_et[next_base_len - 1].v;
03416 }
03417 else if (j + 1 == curr_base_len) {
03418 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03419 u = uv_er[i].u;
03420 v = uv_er[i].v;
03421 }
03422 else {
03423 {
03424 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03425 int nearest_node_j = (int)rel;
03426 rel -= nearest_node_j;
03427 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03428 double u1 = quad->uv_grid[ij].u;
03429 double v1 = quad->uv_grid[ij].v;
03430 double u2 = quad->uv_grid[ij + 1].u;
03431 double v2 = quad->uv_grid[ij + 1].v;
03432 double duj = (u2 - u1) * rel;
03433 double dvj = (v2 - v1) * rel;
03434 u = u1 + duj;
03435 v = v1 + dvj;
03436 }
03437 gp_Pnt P = S->Value(u,v);
03438 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03439 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03440 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03441 Nf = Nf1;
03442 }
03443 next_par_u.SetValue(next_base_len, u);
03444 next_par_v.SetValue(next_base_len, v);
03445 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
03446 NodesBRD.Value(curr_base.Value(j + 1), i),
03447 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03448 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03449 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03450 }
03451
03452 curr_base_len = next_base_len;
03453 curr_base = next_base;
03454 curr_par_u = next_par_u;
03455 curr_par_v = next_par_v;
03456 next_base_len = 0;
03457 }
03458 }
03459 else if (is_lin_31) {
03460
03461
03462
03463
03464
03465
03466
03467
03468
03469
03470
03471
03472
03473
03474
03475
03476
03477 int delta_all = nb - nt;
03478 int delta_one_col = (nr - 1) * 2;
03479 int nb_col = delta_all / delta_one_col;
03480 int remainder = delta_all - nb_col * delta_one_col;
03481 if (remainder > 0) {
03482 nb_col++;
03483 }
03484 int free_left = ((nt - 1) - nb_col) / 2;
03485 free_left += nr - 2;
03486 int free_middle = (nr - 2) * 2;
03487 if (remainder > 0 && nb_col == 1) {
03488 int nb_rows_short_col = remainder / 2;
03489 int nb_rows_thrown = (nr - 1) - nb_rows_short_col;
03490 free_left -= nb_rows_thrown;
03491 }
03492
03493 for (i = 1; i < nr; i++, free_middle -= 2, free_left -= 1) {
03494
03495 NodesBRD.SetValue(1, i+1, uv_el[i].node);
03496 next_base.SetValue(++next_base_len, 1);
03497
03498 NodesBRD.SetValue(nb, i+1, uv_er[i].node);
03499
03500
03501 next_par_u.SetValue(next_base_len, uv_el[i].u);
03502 next_par_v.SetValue(next_base_len, uv_el[i].v);
03503
03504
03505 int nb_next = curr_base_len - nb_col * 2;
03506 if (remainder > 0 && i > remainder / 2)
03507
03508 nb_next += 2;
03509 if (nb_next < nt) nb_next = nt;
03510
03511
03512 for (j = 1; j <= free_left; j++) {
03513
03514 const SMDS_MeshNode* Nf;
03515 double u,v;
03516 next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
03517 if (i + 1 == nr) {
03518 Nf = uv_et[next_base_len - 1].node;
03519 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03520 u = uv_et[next_base_len - 1].u;
03521 v = uv_et[next_base_len - 1].v;
03522 }
03523 else {
03524 {
03525 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03526 int nearest_node_j = (int)rel;
03527 rel -= nearest_node_j;
03528 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03529 double u1 = quad->uv_grid[ij].u;
03530 double v1 = quad->uv_grid[ij].v;
03531 double u2 = quad->uv_grid[ij + 1].u;
03532 double v2 = quad->uv_grid[ij + 1].v;
03533 double duj = (u2 - u1) * rel;
03534 double dvj = (v2 - v1) * rel;
03535 u = u1 + duj;
03536 v = v1 + dvj;
03537 }
03538 gp_Pnt P = S->Value(u,v);
03539 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03540 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03541 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03542 Nf = Nf1;
03543 }
03544 next_par_u.SetValue(next_base_len, u);
03545 next_par_v.SetValue(next_base_len, v);
03546 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
03547 NodesBRD.Value(curr_base.Value(j + 1), i),
03548 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03549 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03550 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03551 }
03552
03553 for (int icol = 1; icol <= nb_col; icol++) {
03554
03555 if (remainder > 0 && icol == nb_col && i > remainder / 2)
03556
03557 break;
03558
03559
03560
03561
03562
03563
03564
03565
03566
03567
03568
03569
03570
03571 double u,v;
03572
03573
03574 const SMDS_MeshNode* Nb;
03575 next_base_len++;
03576 next_base.SetValue(next_base_len, curr_base.Value(j + 3));
03577 if (i + 1 == nr) {
03578 Nb = uv_et[next_base_len - 1].node;
03579 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb);
03580 u = uv_et[next_base_len - 1].u;
03581 v = uv_et[next_base_len - 1].v;
03582 }
03583 else if (j + 3 == curr_base_len) {
03584 Nb = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03585 u = uv_er[i].u;
03586 v = uv_er[i].v;
03587 }
03588 else {
03589 {
03590 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03591 int nearest_node_j = (int)rel;
03592 rel -= nearest_node_j;
03593 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03594 double u1 = quad->uv_grid[ij].u;
03595 double v1 = quad->uv_grid[ij].v;
03596 double u2 = quad->uv_grid[ij + 1].u;
03597 double v2 = quad->uv_grid[ij + 1].v;
03598 double duj = (u2 - u1) * rel;
03599 double dvj = (v2 - v1) * rel;
03600 u = u1 + duj;
03601 v = v1 + dvj;
03602 }
03603 gp_Pnt P = S->Value(u,v);
03604 SMDS_MeshNode* Nb1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03605 meshDS->SetNodeOnFace(Nb1, geomFaceID, u, v);
03606 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nb1);
03607 Nb = Nb1;
03608 }
03609 next_par_u.SetValue(next_base_len, u);
03610 next_par_v.SetValue(next_base_len, v);
03611
03612
03613 double u1 = (curr_par_u.Value(j) + next_par_u.Value(next_base_len - 1)) / 2.0;
03614 double u2 = (curr_par_u.Value(j + 3) + next_par_u.Value(next_base_len)) / 2.0;
03615 double u3 = (u2 - u1) / 3.0;
03616
03617 double v1 = (curr_par_v.Value(j) + next_par_v.Value(next_base_len - 1)) / 2.0;
03618 double v2 = (curr_par_v.Value(j + 3) + next_par_v.Value(next_base_len)) / 2.0;
03619 double v3 = (v2 - v1) / 3.0;
03620
03621
03622 u = u1 + u3;
03623 v = v1 + v3;
03624 gp_Pnt P = S->Value(u,v);
03625 SMDS_MeshNode* Nc = meshDS->AddNode(P.X(), P.Y(), P.Z());
03626 meshDS->SetNodeOnFace(Nc, geomFaceID, u, v);
03627
03628
03629 u = u1 + u3 + u3;
03630 v = v1 + v3 + v3;
03631 P = S->Value(u,v);
03632 SMDS_MeshNode* Ne = meshDS->AddNode(P.X(), P.Y(), P.Z());
03633 meshDS->SetNodeOnFace(Ne, geomFaceID, u, v);
03634
03635
03636 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 0), i),
03637 NodesBRD.Value(curr_base.Value(j + 1), i),
03638 Nc,
03639 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03640 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03641
03642 SMDS_MeshFace* F2 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 1), i),
03643 NodesBRD.Value(curr_base.Value(j + 2), i),
03644 Ne, Nc);
03645 if (F2) meshDS->SetMeshElementOnShape(F2, geomFaceID);
03646
03647 SMDS_MeshFace* F4 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j + 2), i),
03648 NodesBRD.Value(curr_base.Value(j + 3), i),
03649 Nb, Ne);
03650 if (F4) meshDS->SetMeshElementOnShape(F4, geomFaceID);
03651
03652 SMDS_MeshFace* F5 = myHelper->AddFace(Nc, Ne, Nb,
03653 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03654 if (F5) meshDS->SetMeshElementOnShape(F5, geomFaceID);
03655
03656 j += 3;
03657
03658
03659 if (icol < nb_col) {
03660 if (remainder > 0 && icol == nb_col - 1 && i > remainder / 2)
03661
03662 break;
03663
03664 int free_add = free_middle;
03665 if (remainder > 0 && icol == nb_col - 1)
03666
03667 free_add -= (nr - 1) - (remainder / 2);
03668
03669 for (int imiddle = 1; imiddle <= free_add; imiddle++) {
03670
03671 const SMDS_MeshNode* Nf;
03672 double u,v;
03673 next_base.SetValue(++next_base_len, curr_base.Value(j + imiddle));
03674 if (i + 1 == nr) {
03675 Nf = uv_et[next_base_len - 1].node;
03676 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03677 u = uv_et[next_base_len - 1].u;
03678 v = uv_et[next_base_len - 1].v;
03679 }
03680 else if (j + imiddle == curr_base_len) {
03681 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03682 u = uv_er[i].u;
03683 v = uv_er[i].v;
03684 }
03685 else {
03686 {
03687 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03688 int nearest_node_j = (int)rel;
03689 rel -= nearest_node_j;
03690 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03691 double u1 = quad->uv_grid[ij].u;
03692 double v1 = quad->uv_grid[ij].v;
03693 double u2 = quad->uv_grid[ij + 1].u;
03694 double v2 = quad->uv_grid[ij + 1].v;
03695 double duj = (u2 - u1) * rel;
03696 double dvj = (v2 - v1) * rel;
03697 u = u1 + duj;
03698 v = v1 + dvj;
03699 }
03700 gp_Pnt P = S->Value(u,v);
03701 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03702 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03703 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03704 Nf = Nf1;
03705 }
03706 next_par_u.SetValue(next_base_len, u);
03707 next_par_v.SetValue(next_base_len, v);
03708 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j - 1 + imiddle), i),
03709 NodesBRD.Value(curr_base.Value(j + imiddle), i),
03710 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03711 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03712 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03713 }
03714 j += free_add;
03715 }
03716 }
03717
03718
03719 for (; j < curr_base_len; j++) {
03720
03721 const SMDS_MeshNode* Nf;
03722 double u,v;
03723 next_base.SetValue(++next_base_len, curr_base.Value(j + 1));
03724 if (i + 1 == nr) {
03725 Nf = uv_et[next_base_len - 1].node;
03726 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf);
03727 u = uv_et[next_base_len - 1].u;
03728 v = uv_et[next_base_len - 1].v;
03729 }
03730 else if (j + 1 == curr_base_len) {
03731 Nf = NodesBRD.Value(next_base.Value(next_base_len), i + 1);
03732 u = uv_er[i].u;
03733 v = uv_er[i].v;
03734 }
03735 else {
03736 {
03737 double rel = double(next_base_len - 1) * double(nt - 1) / double(nb_next - 1) + 1;
03738 int nearest_node_j = (int)rel;
03739 rel -= nearest_node_j;
03740 int ij = (i + 1 - 1) * nt + (nearest_node_j - 1);
03741 double u1 = quad->uv_grid[ij].u;
03742 double v1 = quad->uv_grid[ij].v;
03743 double u2 = quad->uv_grid[ij + 1].u;
03744 double v2 = quad->uv_grid[ij + 1].v;
03745 double duj = (u2 - u1) * rel;
03746 double dvj = (v2 - v1) * rel;
03747 u = u1 + duj;
03748 v = v1 + dvj;
03749 }
03750 gp_Pnt P = S->Value(u,v);
03751 SMDS_MeshNode* Nf1 = meshDS->AddNode(P.X(), P.Y(), P.Z());
03752 meshDS->SetNodeOnFace(Nf1, geomFaceID, u, v);
03753 NodesBRD.SetValue(next_base.Value(next_base_len), i + 1, Nf1);
03754 Nf = Nf1;
03755 }
03756 next_par_u.SetValue(next_base_len, u);
03757 next_par_v.SetValue(next_base_len, v);
03758 SMDS_MeshFace* F1 = myHelper->AddFace(NodesBRD.Value(curr_base.Value(j), i),
03759 NodesBRD.Value(curr_base.Value(j + 1), i),
03760 NodesBRD.Value(next_base.Value(next_base_len), i + 1),
03761 NodesBRD.Value(next_base.Value(next_base_len - 1), i + 1));
03762 if (F1) meshDS->SetMeshElementOnShape(F1, geomFaceID);
03763 }
03764
03765 curr_base_len = next_base_len;
03766 curr_base = next_base;
03767 curr_par_u = next_par_u;
03768 curr_par_v = next_par_v;
03769 next_base_len = 0;
03770 }
03771 }
03772 else {
03773 }
03774 }
03775
03776 bool isOk = true;
03777 return isOk;
03778 }
03779
03780
03781 namespace
03782 {
03783 struct TSmoothNode;
03784
03790 struct TTriangle
03791 {
03792 TSmoothNode* _n1;
03793 TSmoothNode* _n2;
03794 TTriangle( TSmoothNode* n1=0, TSmoothNode* n2=0 ): _n1(n1), _n2(n2) {}
03795
03796 inline bool IsForward( gp_UV uv ) const;
03797 };
03798
03802 struct TSmoothNode
03803 {
03804 gp_XY _uv;
03805 vector< TTriangle > _triangles;
03806 };
03807
03808 inline bool TTriangle::IsForward( gp_UV uv ) const
03809 {
03810 gp_Vec2d v1( uv, _n1->_uv ), v2( uv, _n2->_uv );
03811 double d = v1 ^ v2;
03812 return d > 1e-100;
03813 }
03814 }
03815
03816
03822
03823
03824 void StdMeshers_Quadrangle_2D::UpdateDegenUV(FaceQuadStruct* quad)
03825 {
03826 for ( unsigned i = 0; i < quad->side.size(); ++i )
03827 {
03828 StdMeshers_FaceSide* side = quad->side[i];
03829 const vector<UVPtStruct>& uvVec = side->GetUVPtStruct();
03830
03831
03832 int degenInd = -1;
03833 if ( myHelper->IsDegenShape( uvVec[0].node->getshapeId() ))
03834 degenInd = 0;
03835 else if ( myHelper->IsDegenShape( uvVec.back().node->getshapeId() ))
03836 degenInd = uvVec.size() - 1;
03837 else
03838 continue;
03839
03840
03841 bool isPrev = ( degenInd == 0 );
03842 if ( i >= TOP_SIDE )
03843 isPrev = !isPrev;
03844 int i2 = ( isPrev ? ( i + 3 ) : ( i + 1 )) % 4;
03845 StdMeshers_FaceSide* side2 = quad->side[ i2 ];
03846 const vector<UVPtStruct>& uvVec2 = side2->GetUVPtStruct();
03847 int degenInd2 = -1;
03848 if ( uvVec[ degenInd ].node == uvVec2[0].node )
03849 degenInd2 = 0;
03850 else if ( uvVec[ degenInd ].node == uvVec2.back().node )
03851 degenInd2 = uvVec2.size() - 1;
03852 else
03853 throw SALOME_Exception( LOCALIZED( "Logical error" ));
03854
03855
03856 uvPtStruct& uv1 = const_cast<uvPtStruct&>( uvVec [ degenInd ]);
03857 uvPtStruct& uv2 = const_cast<uvPtStruct&>( uvVec2[ degenInd2 ]);
03858 uv1.u = uv2.u = 0.5 * ( uv1.u + uv2.u );
03859 uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
03860 }
03861 }
03862
03863
03867
03868
03869 void StdMeshers_Quadrangle_2D::Smooth (FaceQuadStruct* quad)
03870 {
03871 if ( !myNeedSmooth ) return;
03872
03873
03874
03875 typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
03876 TNo2SmooNoMap smooNoMap;
03877
03878 const TopoDS_Face& geomFace = TopoDS::Face( myHelper->GetSubShape() );
03879 SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
03880 SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
03881 SMDS_NodeIteratorPtr nIt = fSubMesh->GetNodes();
03882 while ( nIt->more() )
03883 {
03884 const SMDS_MeshNode* node = nIt->next();
03885 TSmoothNode & sNode = smooNoMap[ node ];
03886 sNode._uv = myHelper->GetNodeUV( geomFace, node );
03887
03888
03889 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator( SMDSAbs_Face );
03890 while ( fIt->more() )
03891 {
03892 const SMDS_MeshElement* face = fIt->next();
03893 const int nbN = face->NbCornerNodes();
03894 const int nInd = face->GetNodeIndex( node );
03895 const int prevInd = myHelper->WrapIndex( nInd - 1, nbN );
03896 const int nextInd = myHelper->WrapIndex( nInd + 1, nbN );
03897 const SMDS_MeshNode* prevNode = face->GetNode( prevInd );
03898 const SMDS_MeshNode* nextNode = face->GetNode( nextInd );
03899 sNode._triangles.push_back( TTriangle( & smooNoMap[ prevNode ],
03900 & smooNoMap[ nextNode ]));
03901 }
03902 }
03903
03904 for ( unsigned i = 0; i < quad->side.size(); ++i )
03905 {
03906 const vector<UVPtStruct>& uvVec = quad->side[i]->GetUVPtStruct();
03907 for ( unsigned j = 0; j < uvVec.size(); ++j )
03908 {
03909 TSmoothNode & sNode = smooNoMap[ uvVec[j].node ];
03910 sNode._uv.SetCoord( uvVec[j].u, uvVec[j].v );
03911 }
03912 }
03913
03914
03915 TNo2SmooNoMap::iterator n2sn = smooNoMap.begin();
03916 for ( ; n2sn != smooNoMap.end(); ++n2sn )
03917 if ( !n2sn->second._triangles.empty() )
03918 break;
03919 if ( n2sn == smooNoMap.end() ) return;
03920 const TSmoothNode & sampleNode = n2sn->second;
03921 const bool refForward = ( sampleNode._triangles[0].IsForward( sampleNode._uv ));
03922
03923
03924
03925 for ( int iLoop = 0; iLoop < 5; ++iLoop )
03926 {
03927 for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
03928 {
03929 TSmoothNode& sNode = n2sn->second;
03930 if ( sNode._triangles.empty() )
03931 continue;
03932
03933
03934 gp_XY newUV (0,0);
03935 for ( unsigned i = 0; i < sNode._triangles.size(); ++i )
03936 newUV += sNode._triangles[i]._n1->_uv;
03937 newUV /= sNode._triangles.size();
03938
03939
03940 bool isValid = true;
03941 for ( unsigned i = 0; i < sNode._triangles.size() && isValid; ++i )
03942 isValid = ( sNode._triangles[i].IsForward( newUV ) == refForward );
03943
03944 if ( isValid )
03945 sNode._uv = newUV;
03946 }
03947 }
03948
03949
03950
03951 Handle(Geom_Surface) surface = BRep_Tool::Surface( geomFace );
03952
03953 for ( n2sn = smooNoMap.begin(); n2sn != smooNoMap.end(); ++n2sn )
03954 {
03955 TSmoothNode& sNode = n2sn->second;
03956 if ( sNode._triangles.empty() )
03957 continue;
03958
03959 SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( n2sn->first );
03960 gp_Pnt xyz = surface->Value( sNode._uv.X(), sNode._uv.Y() );
03961 meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
03962
03963
03964 node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( sNode._uv.X(), sNode._uv.Y() )));
03965 }
03966
03967
03968 if ( _quadraticMesh )
03969 {
03970 const TLinkNodeMap& links = myHelper->GetTLinkNodeMap();
03971 TLinkNodeMap::const_iterator linkIt = links.begin();
03972 for ( ; linkIt != links.end(); ++linkIt )
03973 {
03974 const SMESH_TLink& link = linkIt->first;
03975 SMDS_MeshNode* node = const_cast< SMDS_MeshNode*>( linkIt->second );
03976
03977 if ( node->getshapeId() != myHelper->GetSubShapeID() )
03978 continue;
03979
03980 gp_XY uv1 = myHelper->GetNodeUV( geomFace, link.node1(), node );
03981 gp_XY uv2 = myHelper->GetNodeUV( geomFace, link.node2(), node );
03982
03983 gp_XY uv = myHelper->GetMiddleUV( surface, uv1, uv2 );
03984 node->SetPosition( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
03985
03986 gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
03987 meshDS->MoveNode( node, xyz.X(), xyz.Y(), xyz.Z() );
03988 }
03989 }
03990 }