1 // SMESH SMESH : implementaion of SMESH idl descriptions
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
24 // File : SMESH_Quadrangle_2D.cxx
25 // Author : Paul RASCLE, EDF
31 #include "SMESH_Quadrangle_2D.hxx"
32 #include "SMESH_Gen.hxx"
33 #include "SMESH_Mesh.hxx"
35 #include "SMDS_MeshElement.hxx"
36 #include "SMDS_MeshNode.hxx"
37 #include "SMDS_EdgePosition.hxx"
38 #include "SMDS_FacePosition.hxx"
40 #include <BRep_Tool.hxx>
41 #include <BRepTools.hxx>
42 #include <BRepTools_WireExplorer.hxx>
43 #include <Geom_Surface.hxx>
44 #include <Geom_Curve.hxx>
45 #include <Geom2d_Curve.hxx>
46 #include <Handle_Geom2d_Curve.hxx>
47 #include <Handle_Geom_Curve.hxx>
48 #include <gp_Pnt2d.hxx>
49 #include <TColStd_ListIteratorOfListOfInteger.hxx>
51 #include "utilities.h"
53 //=============================================================================
57 //=============================================================================
59 SMESH_Quadrangle_2D::SMESH_Quadrangle_2D(int hypId,
62 : SMESH_2D_Algo(hypId, studyId, gen)
64 MESSAGE("SMESH_Quadrangle_2D::SMESH_Quadrangle_2D");
65 _name = "Quadrangle_2D";
66 // _shapeType = TopAbs_FACE;
67 _shapeType = (1<<TopAbs_FACE);
70 //=============================================================================
74 //=============================================================================
76 SMESH_Quadrangle_2D::~SMESH_Quadrangle_2D()
78 MESSAGE("SMESH_Quadrangle_2D::~SMESH_Quadrangle_2D");
81 //=============================================================================
85 //=============================================================================
87 bool SMESH_Quadrangle_2D::CheckHypothesis(SMESH_Mesh& aMesh,
88 const TopoDS_Shape& aShape)
90 //MESSAGE("SMESH_Quadrangle_2D::CheckHypothesis");
100 //=============================================================================
104 //=============================================================================
106 bool SMESH_Quadrangle_2D::Compute(SMESH_Mesh& aMesh,
107 const TopoDS_Shape& aShape)
108 throw (SALOME_Exception)
110 //MESSAGE("SMESH_Quadrangle_2D::Compute");
111 const Handle(SMESHDS_Mesh)& meshDS = aMesh.GetMeshDS();
112 SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
114 FaceQuadStruct* quad = CheckAnd2Dcompute(aMesh, aShape);
115 if (!quad) return false;
117 // --- compute 3D values on points, store points & quadrangles
119 int nbdown = quad->nbPts[0];
120 int nbright = quad->nbPts[1];
121 int nbVertices = nbdown*nbright;
122 int nbQuad = (nbdown-1)*(nbright-1);
123 //SCRUTE(nbVertices);
126 // const TopoDS_Face& FF = TopoDS::Face(aShape);
127 // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
128 // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
129 const TopoDS_Face& F = TopoDS::Face(aShape);
130 bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
131 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
133 for (int i=1; i<nbdown-1; i++)
134 for (int j=1; j<nbright-1; j++) // internal points
136 int ij = j*nbdown +i;
137 double u = quad->uv_grid[ij].u;
138 double v = quad->uv_grid[ij].v;
139 gp_Pnt P = S->Value(u,v);
140 int nodeId = meshDS->AddNode(P.X(), P.Y(), P.Z());
141 //MESSAGE("point "<< nodeId<<" "<<" "<<P.X()<<" "<<P.Y()<<" "<<P.Z());
142 Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
143 Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
144 meshDS->SetNodeOnFace(node, F);
145 quad->uv_grid[ij].nodeId = nodeId;
146 // Handle (SMDS_FacePosition) fpos
147 // = new SMDS_FacePosition(theSubMesh->GetId(),i,j); // easier than u,v
148 // node->SetPosition(fpos);
149 Handle (SMDS_FacePosition) fpos
150 = Handle (SMDS_FacePosition)::DownCast(node->GetPosition());
151 fpos->SetUParameter(i);
152 fpos->SetVParameter(j);
155 // bool isQuadForward = ( faceIsForward == quad->isEdgeForward[0]);
156 for (int i=0; i<nbdown-1; i++)
157 for (int j=0; j<nbright-1; j++) // faces
159 int a = quad->uv_grid[ j *nbdown +i ].nodeId;
160 int b = quad->uv_grid[ j *nbdown +i+1].nodeId;
161 int c = quad->uv_grid[(j+1)*nbdown +i+1].nodeId;
162 int d = quad->uv_grid[(j+1)*nbdown +i ].nodeId;
164 // if (isQuadForward) faceId = meshDS->AddFace(a,b,c,d);
165 // else faceId = meshDS->AddFace(a,d,c,b);
166 faceId = meshDS->AddFace(a,b,c,d);
167 Handle (SMDS_MeshElement) elt = meshDS->FindElement(faceId);
168 meshDS->SetMeshElementOnShape(elt, F);
176 //=============================================================================
180 //=============================================================================
183 SMESH_Quadrangle_2D::CheckAnd2Dcompute(SMESH_Mesh& aMesh,
184 const TopoDS_Shape& aShape)
185 throw (SALOME_Exception)
187 //MESSAGE("SMESH_Quadrangle_2D::ComputeWithoutStore");
189 SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
191 // const TopoDS_Face& FF = TopoDS::Face(aShape);
192 // bool faceIsForward = (FF.Orientation() == TopAbs_FORWARD);
193 // TopoDS_Face F = TopoDS::Face(FF.Oriented(TopAbs_FORWARD));
194 const TopoDS_Face& F = TopoDS::Face(aShape);
195 bool faceIsForward = (F.Orientation() == TopAbs_FORWARD);
197 // verify 1 wire only, with 4 edges, same number of points on opposite edges
199 if (NumberOfWires (F) != 1)
201 MESSAGE("only 1 wire by face (quadrangles)");
203 //throw SALOME_Exception(LOCALIZED("only 1 wire by face (quadrangles)"));
205 // const TopoDS_Wire WW = BRepTools::OuterWire(F);
206 // TopoDS_Wire W = TopoDS::Wire(WW.Oriented(TopAbs_FORWARD));
207 const TopoDS_Wire& W = BRepTools::OuterWire(F);
208 BRepTools_WireExplorer wexp(W,F);
210 FaceQuadStruct* quad = new FaceQuadStruct;
211 for (int i=0; i<4; i++) quad->uv_edges[i] = 0;
215 for (wexp.Init(W,F);wexp.More(); wexp.Next())
217 // const TopoDS_Edge& EE = wexp.Current();
218 // TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
219 const TopoDS_Edge& E = wexp.Current();
220 int nb = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
223 quad->edge[nbEdges] = E;
224 quad->nbPts[nbEdges] = nb +2; // internal points + 2 extrema
231 MESSAGE("face must have 4 edges /quadrangles");
234 //throw SALOME_Exception(LOCALIZED("face must have 4 edges /quadrangles"));
237 if (quad->nbPts[0] != quad->nbPts[2])
239 MESSAGE("different point number-opposed edge");
242 //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
245 if (quad->nbPts[1] != quad->nbPts[3])
247 MESSAGE("different point number-opposed edge");
250 //throw SALOME_Exception(LOCALIZED("different point number-opposed edge"));
253 // set normalized grid on unit square in parametric domain
255 SetNormalizedGrid(aMesh, F, quad);
261 //=============================================================================
265 //=============================================================================
267 void SMESH_Quadrangle_2D::QuadDelete(FaceQuadStruct* quad)
269 //MESSAGE("SMESH_Quadrangle_2D::QuadDelete");
272 for (int i=0; i<4; i++)
274 if (quad->uv_edges[i]) delete [] quad->uv_edges[i];
275 quad->edge[i].Nullify();
277 if (quad->uv_grid) delete [] quad->uv_grid;
282 //=============================================================================
286 //=============================================================================
288 void SMESH_Quadrangle_2D::SetNormalizedGrid(SMESH_Mesh& aMesh,
289 const TopoDS_Shape& aShape,
290 FaceQuadStruct* quad)
291 throw (SALOME_Exception)
293 // Algorithme décrit dans "Génération automatique de maillages"
294 // P.L. GEORGE, MASSON, § 6.4.1 p. 84-85
295 // traitement dans le domaine paramétrique 2d u,v
296 // transport - projection sur le carré unité
298 const TopoDS_Face& F = TopoDS::Face(aShape);
300 // 1 --- find orientation of the 4 edges, by test on extrema
303 // |<----north-2-------^ a3 -------------> a2
305 // west-3 east-1 =right | |
309 // v----south-0--------> a0 -------------> a1
314 Handle (Geom2d_Curve) c2d[4];
317 for (int i=0; i<4; i++)
319 c2d[i] = BRep_Tool::CurveOnSurface(quad->edge[i],
323 pf[i] = c2d[i]->Value(quad->first[i]);
324 pl[i] = c2d[i]->Value(quad->last[i]);
325 quad->isEdgeForward[i] = false;
328 double eps2d = 1.e-3; // *** utiliser plutot TopExp::CommonVertex, puis
329 // distances si piece fausse
331 if ((pf[1].Distance(pl[0]) < eps2d) || (pl[1].Distance(pl[0]) < eps2d))
333 quad->isEdgeForward[0] = true;
337 double tmp =quad->first[0];
338 quad->first[0] = quad->last[0];
340 pf[0] = c2d[0]->Value(quad->first[0]);
341 pl[0] = c2d[0]->Value(quad->last[0]);
343 for (int i=1; i<4; i++)
345 quad->isEdgeForward[i] = (pf[i].Distance(pl[i-1]) < eps2d);
346 if (! quad->isEdgeForward[i])
348 double tmp =quad->first[i];
349 quad->first[i] = quad->last[i];
351 pf[i] = c2d[i]->Value(quad->first[i]);
352 pl[i] = c2d[i]->Value(quad->last[i]);
353 //SCRUTE(pf[i].Distance(pl[i-1]));
354 ASSERT(pf[i].Distance(pl[i-1]) < eps2d);
357 //SCRUTE(pf[0].Distance(pl[3]));
358 ASSERT(pf[0].Distance(pl[3]) < eps2d);
360 // for (int i=0; i<4; i++)
362 // SCRUTE(quad->isEdgeForward[i]);
363 // MESSAGE(" -first "<<i<<" "<<pf[i].X()<<" "<<pf[i].Y());
364 // MESSAGE(" -last "<<i<<" "<<pl[i].X()<<" "<<pl[i].Y());
367 // 2 --- load 2d edge points (u,v) with orientation and value on unit square
369 for (int i=0; i<2; i++)
371 quad->uv_edges[i] = LoadEdgePoints(aMesh, F,
376 // quad->isEdgeForward[i]);
378 for (int i=2; i<4; i++)
380 quad->uv_edges[i] = LoadEdgePoints(aMesh, F,
385 // !quad->isEdgeForward[i]);
388 // 3 --- 2D normalized values on unit square [0..1][0..1]
390 int nbdown = quad->nbPts[0];
391 int nbright = quad->nbPts[1];
392 quad->uv_grid = new UVPtStruct[nbright*nbdown];
394 UVPtStruct* uv_grid = quad->uv_grid;
395 UVPtStruct* uv_e0 = quad->uv_edges[0];
396 UVPtStruct* uv_e1 = quad->uv_edges[1];
397 UVPtStruct* uv_e2 = quad->uv_edges[2];
398 UVPtStruct* uv_e3 = quad->uv_edges[3];
407 for (int i=0; i<nbdown; i++)
409 int ij = j*nbdown +i;
410 uv_grid[ij].nodeId = uv_e0[i].nodeId;
413 for (int j=0; j<nbright; j++)
415 int ij = j*nbdown +i;
416 uv_grid[ij].nodeId = uv_e1[j].nodeId;
419 for (int i=0; i<nbdown; i++)
421 int ij = j*nbdown +i;
422 uv_grid[ij].nodeId = uv_e2[i].nodeId;
425 for (int j=0; j<nbright; j++)
427 int ij = j*nbdown +i;
428 uv_grid[ij].nodeId = uv_e3[j].nodeId;
431 // normalized 2d values on grid
433 for (int i=0; i<nbdown; i++)
434 for (int j=0; j<nbright; j++)
436 int ij = j*nbdown +i;
437 // --- droite i cste : x = x0 + y(x1-x0)
438 double x0 = uv_e0[i].normParam; // bas - sud
439 double x1 = uv_e2[i].normParam; // haut - nord
440 // --- droite j cste : y = y0 + x(y1-y0)
441 double y0 = uv_e3[j].normParam; // gauche-ouest
442 double y1 = uv_e1[j].normParam; // droite - est
443 // --- intersection : x=x0+(y0+x(y1-y0))(x1-x0)
444 double x=(x0+y0*(x1-x0))/(1-(y1-y0)*(x1-x0));
445 double y=y0+x*(y1-y0);
448 //MESSAGE("-xy-01 "<<x0<<" "<<x1<<" "<<y0<<" "<<y1);
449 //MESSAGE("-xy-norm "<<i<<" "<<j<<" "<<x<<" "<<y);
452 // 4 --- projection on 2d domain (u,v)
454 for (int i=0; i<nbdown; i++)
455 for (int j=0; j<nbright; j++)
457 int ij = j*nbdown +i;
458 double x = uv_grid[ij].x;
459 double y = uv_grid[ij].y;
460 double param_0 = uv_e0[0].param
461 + x*(uv_e0[nbdown-1].param -uv_e0[0].param); // sud
462 double param_2 = uv_e2[0].param
463 + x*(uv_e2[nbdown-1].param -uv_e2[0].param); // nord
464 double param_1 = uv_e1[0].param
465 + y*(uv_e1[nbright-1].param -uv_e1[0].param); // est
466 double param_3 = uv_e3[0].param
467 + y*(uv_e3[nbright-1].param -uv_e3[0].param); // ouest
469 //MESSAGE("params "<<param_0<<" "<<param_1<<" "<<param_2<<" "<<param_3);
470 gp_Pnt2d p0 = c2d[0]->Value(param_0);
471 gp_Pnt2d p1 = c2d[1]->Value(param_1);
472 gp_Pnt2d p2 = c2d[2]->Value(param_2);
473 gp_Pnt2d p3 = c2d[3]->Value(param_3);
475 double u = (1-y)*p0.X() + x*p1.X() + y*p2.X() + (1-x)*p3.X();
476 double v = (1-y)*p0.Y() + x*p1.Y() + y*p2.Y() + (1-x)*p3.Y();
478 u -= (1-x)*(1-y)*a0.X() + x*(1-y)*a1.X() + x*y*a2.X() + (1-x)*y*a3.X();
479 v -= (1-x)*(1-y)*a0.Y() + x*(1-y)*a1.Y() + x*y*a2.Y() + (1-x)*y*a3.Y();
484 //MESSAGE("-uv- "<<i<<" "<<j<<" "<<uv_grid[ij].u<<" "<<uv_grid[ij].v);
488 //=============================================================================
492 //=============================================================================
494 UVPtStruct* SMESH_Quadrangle_2D::LoadEdgePoints(SMESH_Mesh& aMesh,
495 const TopoDS_Face& F,
496 const TopoDS_Edge& E,
501 //MESSAGE("SMESH_Quadrangle_2D::LoadEdgePoints");
503 Handle (SMDS_Mesh) meshDS = aMesh.GetMeshDS();
505 // --- IDNodes of first and last Vertex
507 TopoDS_Vertex VFirst, VLast;
508 TopExp::Vertices(E, VFirst, VLast); // corresponds to f and l
510 ASSERT(!VFirst.IsNull());
511 SMESH_subMesh* firstSubMesh = aMesh.GetSubMesh(VFirst);
512 const TColStd_ListOfInteger& lidf
513 = firstSubMesh->GetSubMeshDS()->GetIDNodes();
514 int idFirst= lidf.First();
517 ASSERT(!VLast.IsNull());
518 SMESH_subMesh* lastSubMesh = aMesh.GetSubMesh(VLast);
519 const TColStd_ListOfInteger& lidl
520 = lastSubMesh->GetSubMeshDS()->GetIDNodes();
521 int idLast= lidl.First();
524 // --- edge internal IDNodes (relies on good order storage, not checked)
526 int nbPoints = aMesh.GetSubMesh(E)->GetSubMeshDS()->NbNodes();
528 UVPtStruct * uvslf = new UVPtStruct[nbPoints+2];
531 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E,F,f,l);
533 const TColStd_ListOfInteger& indElt
534 = aMesh.GetSubMesh(E)->GetSubMeshDS()->GetIDNodes();
535 TColStd_ListIteratorOfListOfInteger ite(indElt);
537 //SCRUTE(indElt.Extent());
538 ASSERT(nbPoints == indElt.Extent());
540 map<double,int> params;
541 for (; ite.More(); ite.Next())
543 int nodeId = ite.Value();
544 Handle (SMDS_MeshElement) elt = meshDS->FindNode(nodeId);
545 Handle (SMDS_MeshNode) node = meshDS->GetNode(1, elt);
546 Handle (SMDS_EdgePosition) epos
547 = Handle (SMDS_EdgePosition)::DownCast(node->GetPosition());
548 double param = epos->GetUParameter();
549 params[param] = nodeId;
552 bool isForward = (((l-f)*(last-first)) > 0);
559 gp_Pnt2d p = C2d->Value(f); // first point = Vertex Forward
563 uvslf [0].nodeId = idFirst;
564 //MESSAGE("__ f "<<f<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
565 map<double,int>::iterator itp = params.begin();
566 for (int i = 1; i <= nbPoints; i++) // nbPoints internal
568 double param = (*itp).first;
569 int nodeId = (*itp).second;
570 gp_Pnt2d p = C2d->Value(param);
573 uvslf[i].param = param;
574 uvslf[i].nodeId = nodeId;
575 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
578 p = C2d->Value(l); // last point = Vertex Reversed
579 uvslf [nbPoints+1].x = p.X();
580 uvslf [nbPoints+1].y = p.Y();
581 uvslf [nbPoints+1].param = l;
582 uvslf [nbPoints+1].nodeId = idLast;
583 //MESSAGE("__ l "<<l<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
589 gp_Pnt2d p = C2d->Value(l); // first point = Vertex Reversed
593 uvslf [0].nodeId = idLast;
594 //MESSAGE("__ l "<<l<<" "<<uvslf[0].x <<" "<<uvslf[0].y);
595 map<double,int>::reverse_iterator itp = params.rbegin();
596 for (int j = nbPoints; j >= 1; j--) // nbPoints internal
598 double param = (*itp).first;
599 int nodeId = (*itp).second;
600 int i = nbPoints +1 -j;
601 gp_Pnt2d p = C2d->Value(param);
604 uvslf[i].param = param;
605 uvslf[i].nodeId = nodeId;
606 //MESSAGE("__ "<<i<<" "<<param<<" "<<uvslf[i].x <<" "<<uvslf[i].y);
609 p = C2d->Value(f); // last point = Vertex Forward
610 uvslf [nbPoints+1].x = p.X();
611 uvslf [nbPoints+1].y = p.Y();
612 uvslf [nbPoints+1].param = f;
613 uvslf [nbPoints+1].nodeId = idFirst;
614 //MESSAGE("__ f "<<f<<" "<<uvslf[nbPoints+1].x <<" "<<uvslf[nbPoints+1].y);
617 ASSERT(paramin != paramax);
618 for (int i = 0; i< nbPoints+2; i++)
620 uvslf[i].normParam = (uvslf[i].param -paramin)/(paramax -paramin);
621 //SCRUTE(uvslf[i].normParam);
627 //=============================================================================
631 //=============================================================================
633 ostream & SMESH_Quadrangle_2D::SaveTo(ostream & save)
638 //=============================================================================
642 //=============================================================================
644 istream & SMESH_Quadrangle_2D::LoadFrom(istream & load)
646 return load >> (*this);
649 //=============================================================================
653 //=============================================================================
655 ostream & operator << (ostream & save, SMESH_Quadrangle_2D & hyp)
660 //=============================================================================
664 //=============================================================================
666 istream & operator >> (istream & load, SMESH_Quadrangle_2D & hyp)