1 // Copyright (C) 2009-2023 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : SMESH_HexaBlocks.cxx
29 #include <Precision.hxx>
30 #include <BRep_Tool.hxx>
31 #include <BRepTools.hxx>
32 #include <BRep_Builder.hxx>
33 #include <BRepAdaptor_Curve.hxx>
34 #include <BRepBuilderAPI_MakeFace.hxx>
36 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
37 #include <GCPnts_AbscissaPoint.hxx>
38 #include <TopoDS_Wire.hxx>
41 #include <TopoDS_Shape.hxx>
42 #include <TopoDS_Compound.hxx>
47 #include <IntCurvesFace_ShapeIntersector.hxx>
50 #include "SMDS_MeshNode.hxx"
51 #include "SMDS_MeshVolume.hxx"
52 #include "SMESH_Gen.hxx"
53 #include "SMESH_MesherHelper.hxx"
54 #include "SMESHDS_Group.hxx"
57 #include "HexDocument.hxx"
58 #include "HexVertex.hxx"
59 #include "HexEdge.hxx"
60 #include "HexQuad.hxx"
61 #include "HexHexa.hxx"
62 #include "HexPropagation.hxx"
63 #include "HexGroup.hxx"
64 #include "hexa_base.hxx"
65 #include "HexAssoEdge.hxx"
67 // HEXABLOCKPLUGIN includes
68 #include "HEXABLOCKPlugin_mesh.hxx"
69 #include "HEXABLOCKPlugin_FromSkin_3D.hxx"
72 #include "Basics_Utils.hxx"
73 #include "utilities.h"
84 #define EXCEPTION(TYPE, MSG) {\
85 std::ostringstream aStream;\
86 aStream<<__FILE__<<"["<<__LINE__<<"]::"<<MSG;\
87 throw TYPE(aStream.str());\
92 static double HEXA_EPS = 1.0e-3; //1E-3;
93 static double HEXA_QUAD_WAY = M_PI/4.; //3.*PI/8.;
95 // ============================================================ string2shape
96 TopoDS_Shape string2shape( const std::string& brep )
99 std::istringstream streamBrep(brep);
100 BRep_Builder aBuilder;
101 BRepTools::Read(shape, streamBrep, aBuilder);
105 // ============================================================ shape2coord
106 bool shape2coord(const TopoDS_Shape& aShape, double& x, double& y, double& z)
108 if ( aShape.ShapeType() == TopAbs_VERTEX ){
109 TopoDS_Vertex aPoint;
110 aPoint = TopoDS::Vertex( aShape );
111 gp_Pnt aPnt = BRep_Tool::Pnt( aPoint );
120 // ============================================================= edge_length
121 double edge_length(const TopoDS_Edge & E)
123 double UMin = 0, UMax = 0;
124 if (BRep_Tool::Degenerated(E))
127 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
128 GeomAdaptor_Curve AdaptCurve(C);
129 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
132 // =============================================================== make_curve
133 BRepAdaptor_Curve* make_curve (gp_Pnt& pstart, gp_Pnt& pend,
134 double& length, double& u_start,
135 HEXA_NS::Edge& edge, int nro)
137 Hex::AssoEdge* asso = edge.getAssociation (nro);
138 BRepAdaptor_Curve* the_curve = asso->getCurve();
139 length = asso->length ();
143 const double* point1 = asso->getOrigin ();
144 const double* point2 = asso->getExtrem ();
145 pstart = gp_Pnt (point1[0], point1[1], point1[2]);
146 pend = gp_Pnt (point2[0], point2[1], point2[2]);
147 u_start = asso->getUstart();
150 // ============================================================== Constructeur
151 // SMESH_HexaBlocks::SMESH_HexaBlocks( SMESH_Mesh* theMesh ):
152 SMESH_HexaBlocks::SMESH_HexaBlocks(SMESH_Mesh& theMesh):
156 _computeVertexOK(false),
157 _computeEdgeOK(false),
158 _computeQuadOK(false),
159 _theMesh(&theMesh), //groups creation
160 _theMeshDS(theMesh.GetMeshDS()) //meshing
165 SMESH_HexaBlocks::~SMESH_HexaBlocks()
170 // --------------------------------------------------------------
172 // --------------------------------------------------------------
174 // --------------------------------------------------------------
176 // --------------------------------------------------------------
177 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
178 // ============================================================== computeVertex
179 bool SMESH_HexaBlocks::computeVertex(HEXA_NS::Vertex& vx)
182 vx.getAssoCoord (px, py, pz);
184 SMDS_MeshNode* new_node = _theMeshDS->AddNode (px, py, pz);
185 _vertex [new_node] = &vx;
186 _node [&vx] = new_node; //needed in computeEdge()
187 _computeVertexOK = true;
190 // --------------------------------------------------------------
192 // --------------------------------------------------------------
193 bool SMESH_HexaBlocks::computeEdge(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
197 ok = computeEdgeByAssoc( edge, law);
199 ok = computeEdgeByShortestWire( edge, law);
202 ok = computeEdgeByPlanWire( edge, law);
205 ok = computeEdgeByIsoWire( edge, law);
208 ok = computeEdgeBySegment( edge, law);
211 _computeEdgeOK = true;
218 bool SMESH_HexaBlocks::computeEdgeByAssoc( HEXA_NS::Edge& edge, HEXA_NS::Law& law )
220 MESSAGE("computeEdgeByAssoc(edgeID = "<<edge.getId()<<"): begin <<<<<<");
221 ASSERT( _computeVertexOK );
224 if (NOT edge.isAssociated())
228 HEXA_NS::Vertex* vx0 = NULL;
229 HEXA_NS::Vertex* vx1 = NULL;
231 // way of discretization
232 if (edge.getWay() == true){
233 vx0 = edge.getVertex(0);
234 vx1 = edge.getVertex(1);
236 vx0 = edge.getVertex(1);
237 vx1 = edge.getVertex(0);
240 SMDS_MeshNode* FIRST_NODE = _node[vx0];
241 SMDS_MeshNode* LAST_NODE = _node[vx1];
244 // A) Build myCurve_list
245 std::list< BRepAdaptor_Curve* > myCurve_list;
246 double myCurve_tot_len;
247 std::map< BRepAdaptor_Curve*, double> myCurve_lengths;
248 std::map< BRepAdaptor_Curve*, bool> myCurve_ways;
249 std::map< BRepAdaptor_Curve*, double> myCurve_starts;
251 gp_Pnt myCurve_pt_start(FIRST_NODE->X(), FIRST_NODE->Y(), FIRST_NODE->Z());
252 gp_Pnt myCurve_pt_end (LAST_NODE->X(), LAST_NODE->Y(), LAST_NODE->Z());
266 // B) Build nodes and edges on mesh from myCurve_list
267 SMDS_MeshNode* node_a = NULL;
268 SMDS_MeshNode* node_b = NULL;
269 SMDS_MeshEdge* edge_ab = NULL;
270 SMESHNodes nodesOnEdge;
271 SMESHEdges edgesOnEdge; //backup for group creation
275 nodesOnEdge.push_back(FIRST_NODE);
276 // nodesXxOnEdge.push_back(0.);
277 // _nodeXx[FIRST_NODE] = 0.;
281 double myCurve_start_u = 0.;
282 int nbNodes = law.getNodes(); //law of discretization
283 if (myCurve_list.size()==0)
285 PutData (edge.getName());
288 MESSAGE("nbNodes -> "<<nbNodes);
289 for (int i = 0; i < nbNodes; ++i){
290 u = _Xx(i, law, nbNodes); //u between [0,1]
291 myCurve_u = u*myCurve_tot_len;
294 MESSAGE("myCurve_u -> "<<myCurve_u);
295 MESSAGE("myCurve_tot_len -> "<<myCurve_tot_len);
297 ptOnMyCurve = _getPtOnMyCurve( myCurve_u,
305 node_b = _theMeshDS->AddNode( ptOnMyCurve.X(), ptOnMyCurve.Y(), ptOnMyCurve.Z() );
306 edge_ab = _theMeshDS->AddEdge( node_a, node_b );
307 nodesOnEdge.push_back( node_b );
308 edgesOnEdge.push_back( edge_ab );
309 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
313 edge_ab = _theMeshDS->AddEdge( node_a, LAST_NODE );
314 nodesOnEdge.push_back( LAST_NODE );
315 edgesOnEdge.push_back( edge_ab );
316 _nodesOnEdge[&edge] = nodesOnEdge;
317 _edgesOnEdge[&edge] = edgesOnEdge;
320 MESSAGE("computeEdgeByAssoc() : end >>>>>>>>");
327 bool SMESH_HexaBlocks::computeEdgeByShortestWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
329 MESSAGE("computeEdgeByShortestWire() not implemented");
333 bool SMESH_HexaBlocks::computeEdgeByPlanWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
335 MESSAGE("computeEdgeByPlanWire() not implemented");
339 bool SMESH_HexaBlocks::computeEdgeByIsoWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
341 MESSAGE("computeEdgeByIsoWire() not implemented");
346 bool SMESH_HexaBlocks::computeEdgeBySegment(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
348 MESSAGE("computeEdgeBySegment() : : begin <<<<<<");
349 ASSERT( _computeVertexOK );
353 HEXA_NS::Vertex* vx0 = NULL;
354 HEXA_NS::Vertex* vx1 = NULL;
356 // way of discretization
357 if (edge.getWay() == true){
358 vx0 = edge.getVertex(0);
359 vx1 = edge.getVertex(1);
361 vx0 = edge.getVertex(1);
362 vx1 = edge.getVertex(0);
366 SMDS_MeshNode* FIRST_NODE = _node[vx0];
367 SMDS_MeshNode* LAST_NODE = _node[vx1];
368 SMDS_MeshNode* node_a = NULL; //new node (to be added)
369 SMDS_MeshNode* node_b = NULL; //new node (to be added)
371 // node and edge creation
372 SMESHNodes nodesOnEdge;
373 SMESHEdges edgesOnEdge;
376 double newNodeX, newNodeY, newNodeZ;
377 SMDS_MeshEdge* newEdge = NULL;
380 nodesOnEdge.push_back(FIRST_NODE);
382 //law of discretization
383 int nbNodes = law.getNodes();
384 MESSAGE("nbNodes -> "<<nbNodes);
385 for (int i = 0; i < nbNodes; ++i){
386 u = _Xx(i, law, nbNodes);
387 newNodeX = FIRST_NODE->X() + u * ( LAST_NODE->X() - FIRST_NODE->X() );
388 newNodeY = FIRST_NODE->Y() + u * ( LAST_NODE->Y() - FIRST_NODE->Y() );
389 newNodeZ = FIRST_NODE->Z() + u * ( LAST_NODE->Z() - FIRST_NODE->Z() );
390 node_b = _theMeshDS->AddNode(newNodeX, newNodeY, newNodeZ);
391 newEdge = _theMeshDS->AddEdge(node_a, node_b);
392 edgesOnEdge.push_back(newEdge);
393 nodesOnEdge.push_back(node_b);
394 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
395 _nodeXx[ node_b ] = u;
396 MESSAGE("_nodeXx <-"<<u);
399 newEdge = _theMeshDS->AddEdge(node_a, LAST_NODE);
400 nodesOnEdge.push_back(LAST_NODE);
401 edgesOnEdge.push_back(newEdge);
403 _nodesOnEdge[&edge] = nodesOnEdge;
404 _edgesOnEdge[&edge] = edgesOnEdge;
406 MESSAGE("computeEdgeBySegment() : end >>>>>>>>");
411 // --------------------------------------------------------------
413 // --------------------------------------------------------------
414 std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document* doc )
416 std::map<HEXA_NS::Quad*, bool> quadWays;
417 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
418 std::list<HEXA_NS::Quad*> skinQuad;
419 std::list<HEXA_NS::Quad*> workingQuad;
420 HEXA_NS::Quad* first_q = NULL;
421 HEXA_NS::Quad* q = NULL;
422 HEXA_NS::Edge* e = NULL;
423 HEXA_NS::Vertex *e_0, *e_1 = NULL;
425 // FIRST STEP: eliminate free quad + internal quad
426 int nTotalQuad = doc->countUsedQuad();
427 for (int i=0; i < nTotalQuad; ++i ){
428 q = doc->getUsedQuad(i);
429 switch ( q->getNbrParents() ){ // parent == hexaedron
430 case 0: case 2: quadWays[q] = true; break;
431 case 1: skinQuad.push_back(q); break;
432 default: if ( q->getNbrParents() > 2 ) ASSERT(false);
436 // SECOND STEP: setting edges ways
437 while ( skinQuad.size()>0 ){
438 MESSAGE("SEARCHING INITIAL QUAD ..." );
439 for ( std::list<HEXA_NS::Quad*>::iterator it = skinQuad.begin(); it != skinQuad.end(); it++ ){
440 _searchInitialQuadWay( *it, e_0, e_1 );
441 if ( e_0 != NULL && e_1 != NULL ){
446 if ( e_0 == NULL && e_1 == NULL ) ASSERT(false);// should never happened,
447 MESSAGE("INITIAL QUAD FOUND!" );
448 for ( int j=0 ; j < 4 ; ++j ){
450 if ( ((e_0 == e->getVertex(0)) && (e_1 == e->getVertex(1)))
451 || ((e_0 == e->getVertex(1)) && (e_1 == e->getVertex(0))) ){
455 MESSAGE("INITIAL EDGE WAY FOUND!" );
457 edgeWays[e] = std::make_pair( e_0, e_1 );
458 workingQuad.push_back(q);
460 while ( workingQuad.size() > 0 ){
461 MESSAGE("COMPUTE QUAD WAY ... ID ="<< q->getId());
462 HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
464 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
465 while ( localEdgeWays.size() != 4 ){
466 HEXA_NS::Edge* e = q->getEdge(i%4);
467 if ( lastVertex == NULL ){
468 if ( edgeWays.count(e) == 1 ){
470 localEdgeWays[e] = std::make_pair( edgeWays[e].first, edgeWays[e].second );
472 localEdgeWays[e] = std::make_pair( edgeWays[e].second, edgeWays[e].first);
474 firstVertex = localEdgeWays[e].first;
475 lastVertex = localEdgeWays[e].second;
478 HEXA_NS::Vertex* e_0 = e->getVertex(0);
479 HEXA_NS::Vertex* e_1 = e->getVertex(1);
481 if ( lastVertex == e_0 ){
482 firstVertex = e_0; lastVertex = e_1;
483 } else if ( lastVertex == e_1 ){
484 firstVertex = e_1; lastVertex = e_0;
485 } else if ( firstVertex == e_0 ) {
486 firstVertex = e_1; lastVertex = e_0;
487 } else if ( firstVertex == e_1 ) {
488 firstVertex = e_0; lastVertex = e_1;
492 localEdgeWays[e] = std::make_pair( firstVertex, lastVertex );
493 if ( edgeWays.count(e) == 0 ){ // keep current value if present otherwise add it
494 edgeWays[e] = localEdgeWays[e];
501 //add new working quad
502 for ( int i=0 ; i < 4; ++i ){
503 HEXA_NS::Quad* next_q = NULL;
504 HEXA_NS::Edge* e = q->getEdge(i);
505 for ( int j=0 ; j < e->getNbrParents(); ++j ){
506 next_q = e->getParent(j);
510 int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
512 int fromWorkingQuad = std::count( workingQuad.begin(), workingQuad.end(), next_q );
513 if ( fromWorkingQuad == 0 ){
514 workingQuad.push_front( next_q );
521 HEXA_NS::Edge* e0 = q->getEdge(0);
522 HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
524 if ( e0_0 == localEdgeWays[ e0 ].first ){
526 } else if ( e0_0 == localEdgeWays[ e0 ].second ){
531 workingQuad.remove( q );
532 skinQuad.remove( q );
533 q = workingQuad.back();
544 // std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document& doc, std::map<HEXA_NS::Quad*, bool> initQuads )
546 // std::map<HEXA_NS::Quad*, bool> quadWays;
547 // // std::map<HEXA_NS::Edge*, bool> edgeWays;
548 // // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
549 // std::map<HEXA_NS::Quad*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > workingQuads;
551 // std::list<HEXA_NS::Quad*> skinQuad;
552 // std::list<HEXA_NS::Quad*> notSkinQuad;
553 // // std::list<HEXA_NS::Quad*> workingQuad;
554 // HEXA_NS::Quad* first_q = NULL;
555 // HEXA_NS::Quad* q = NULL;
556 // HEXA_NS::Edge* e = NULL;
557 // HEXA_NS::Vertex *e_0, *e_1 = NULL;
559 // // FIRST STEP: eliminate free quad + internal quad
560 // int nTotalQuad = doc.countQuad();
561 // for (int i=0; i < nTotalQuad; ++i ){
562 // q = doc.getQuad(i);
563 // switch ( q->getNbrParents() ){ // parent == hexaedron
564 // case 0: case 2: quadWays[q] = true; break;
565 // // case 0: case 2: notSkinQuad.push_back(q); break; //CS_TEST
566 // case 1: skinQuad.push_back(q); break;
567 // default: if ( q->getNbrParents() > 2 ) ASSERT(false);
573 // q = first_q = skinQuad.front();
574 // e = q->getUsedEdge(0);
575 // e_0 = e->getVertex(0);
576 // e_1 = e->getVertex(1);
578 // workingQuads[q] = std::make_pair( e_0, e_1 );
580 // while ( workingQuads.size() > 0 ){
581 // MESSAGE("CURRENT QUAD ------>"<< q->getId());
582 // HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
585 // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
586 // while ( localEdgeWays.size() != 4 ){
587 // HEXA_NS::Edge* e = q->getUsedEdge(i%4);
588 // if ( lastVertex == NULL ){
589 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
590 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
592 // if ( (workingQuads[q].first == e_0 and workingQuads[q].second == e_1)
593 // or (workingQuads[q].first == e_1 and workingQuads[q].second == e_0) ){
594 // if ( q == first_q ){
595 // localEdgeWays[e] = std::make_pair( workingQuads[q].first, workingQuads[q].second );
596 // MESSAGE("FIRST QUAD ");
598 // localEdgeWays[e] = std::make_pair( workingQuads[q].second, workingQuads[q].first);
599 // MESSAGE("NOT FIRST QUAD ");
601 // firstVertex = localEdgeWays[e].first;
602 // lastVertex = localEdgeWays[e].second;
605 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
606 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
607 // if ( lastVertex == e_0 ){
608 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
609 // firstVertex = e_0;
611 // } else if ( lastVertex == e_1 ){
612 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
613 // firstVertex = e_1;
615 // } else if ( firstVertex == e_0 ) {
616 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
617 // firstVertex = e_1;
619 // } else if ( firstVertex == e_1 ) {
620 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
621 // firstVertex = e_0;
631 // //add new working quad
632 // for ( int i=0 ; i < 4; ++i ){
633 // HEXA_NS::Quad* next_q = NULL;
634 // HEXA_NS::Edge* e = q->getUsedEdge(i);
635 // MESSAGE("NB PARENTS ="<< e->getNbrParents() );
636 // for ( int j=0 ; j < e->getNbrParents(); ++j ){
637 // next_q = e->getParent(j);
638 // if ( next_q == q ){
641 // int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
642 // if (fromSkin != 0){
643 // // int fromWorkingQuad = std::count( workingQuads.begin(), workingQuads.end(), next_q );
644 // int fromWorkingQuad = workingQuads.count( next_q );
645 // // MESSAGE("CHECK QUAD:"<< newWorkingQuad->getId());
646 // if ( fromWorkingQuad == 0 ){
647 // // workingQuads.push_front( next_q );
648 // workingQuads[ next_q ] = localEdgeWays[e];
649 // // MESSAGE("EDGE :"<<e->getId()<<"ADD QUAD :"<< newWorkingQuad->getId());
655 // //setting quad way
656 // HEXA_NS::Edge* e0 = q->getUsedEdge(0);
657 // HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
659 // if ( e0_0 == localEdgeWays[ e0 ].first ){
660 // quadWays[q] = true;
661 // } else if ( e0_0 == localEdgeWays[ e0 ].second ){
662 // quadWays[q] = false;
666 // MESSAGE("quadWays ID ="<< q->getId() << ", WAY = " << quadWays[q] );
668 // // workingQuad.remove( q );
669 // workingQuads.erase( q );
670 // skinQuad.remove( q );
671 // *workingQuads.begin();
672 // q = (*workingQuads.begin()).first;
678 bool SMESH_HexaBlocks::computeQuad( HEXA_NS::Quad& quad, bool way )
682 ok = computeQuadByAssoc(quad, way);
684 ok = computeQuadByFindingGeom(quad, way);
687 ok = computeQuadByLinearApproximation(quad, way);
690 _computeQuadOK = true;
696 bool SMESH_HexaBlocks::computeQuadByAssoc( HEXA_NS::Quad& quad, bool way )
698 // int id = quad.getId();
699 // if ( id != 11 ) return false; //7
701 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
702 MESSAGE("quadID = "<<quad.getId());
704 ASSERT( _computeEdgeOK );
707 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
708 SMESHFaces facesOnQuad;
709 SMDS_MeshFace* newFace = NULL;
710 std::vector<double> xx, yy;
712 // Elements for quad computation
713 SMDS_MeshNode *S1, *S2, *S4, *S3;
715 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
716 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
717 if ( initOk == false ){
721 int nbass = quad.countAssociation ();
724 MESSAGE("computeQuadByAssoc() : end >>>>>>>>");
728 TopoDS_Shape shapeOrCompound = getFaceShapes ( quad );
731 std::map<SMDS_MeshNode*, gp_Pnt> interpolatedPoints;
732 int iSize = nodesOnQuad.size();
733 int jSize = nodesOnQuad[0].size();
735 S1 = nodesOnQuad[0][0];
736 S2 = nodesOnQuad[iSize-1][0];
737 S4 = nodesOnQuad[0][jSize-1];
738 S3 = nodesOnQuad[iSize-1][jSize-1];
741 for (int j = 1; j < jSize; ++j){
742 for (int i = 1; i < iSize; ++i){
743 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
744 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
745 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
746 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
749 double newNodeX, newNodeY, newNodeZ;
750 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
751 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
752 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
753 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
757 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
758 gp_Pnt newPt = gp_Pnt( newNodeX, newNodeY, newNodeZ );//interpolated point
761 if ( interpolatedPoints.count(n1) > 0 ){
762 pt1 = interpolatedPoints[n1];
764 pt1 = gp_Pnt( n1->X(), n1->Y(), n1->Z() );
766 if ( interpolatedPoints.count(n3) > 0 ){
767 pt3 = interpolatedPoints[n3];
769 pt3 = gp_Pnt( n3->X(), n3->Y(), n3->Z() );
771 gp_Vec vec1( newPt, pt1 );
772 gp_Vec vec2( newPt, pt3 );
774 gp_Pnt ptOnShape = _intersect(newPt, vec1, vec2, shapeOrCompound);
775 newNodeX = ptOnShape.X();
776 newNodeY = ptOnShape.Y();
777 newNodeZ = ptOnShape.Z();
778 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
779 nodesOnQuad[i][j] = n4;
780 interpolatedPoints[ n4 ] = newPt;
782 MESSAGE("u parameter is "<<u);
783 MESSAGE("v parameter is "<<v);
784 MESSAGE("point interpolated ("<<newPt.X()<<","<<newPt.Y()<<","<<newPt.Z()<<" )");
785 MESSAGE("point on shape ("<<newNodeX<<","<<newNodeY<<","<<newNodeZ<<" )");
788 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
789 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
790 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
791 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
794 MESSAGE("AddFace( n1, n2, n3, n4 )");
795 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
797 MESSAGE("AddFace( n4, n3, n2, n1 )");
798 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
800 facesOnQuad.push_back(newFace);
803 _quadNodes[ &quad ] = nodesOnQuad;
804 _facesOnQuad[&quad] = facesOnQuad;
806 MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
811 bool SMESH_HexaBlocks::computeQuadByFindingGeom( HEXA_NS::Quad& quad, bool way )
813 MESSAGE("computeQuadByFindingGeom() not implemented");
817 bool SMESH_HexaBlocks::_computeQuadInit(
819 ArrayOfSMESHNodes& nodesOnQuad,
820 std::vector<double>& xx, std::vector<double>& yy)
822 MESSAGE("_computeQuadInit() : begin ---------------");
825 SMDS_MeshNode *S1, *S2, *S4, *S3;
826 HEXA_NS::Edge *eh, *eb, *eg, *ed;
827 HEXA_NS::Edge *e1, *e2, *e3, *e4;
828 HEXA_NS::Vertex *e1_0, *e1_1, *e2_0, *e2_1, *e3_0, *e3_1, *e4_0, *e4_1;
830 e1 = quad.getEdge(0);
831 e2 = quad.getEdge(1);
832 e3 = quad.getEdge(2);
833 e4 = quad.getEdge(3);
835 e1_0 = e1->getVertex(0); e1_1 = e1->getVertex(1);
836 e2_0 = e2->getVertex(0); e2_1 = e2->getVertex(1);
837 e3_0 = e3->getVertex(0); e3_1 = e3->getVertex(1);
838 e4_0 = e4->getVertex(0); e4_1 = e4->getVertex(1);
841 S1 = _node[e1_0]; S2 = _node[e1_1];
847 } else if ( e1_0 == e2_1 ){
850 } else if ( e1_0 == e4_0 ){
853 } else if ( e1_0 == e4_1 ){
860 if ( S4 == _node[e3_0] ){
862 } else if ( S4 == _node[e3_1] ){
868 SMESHNodes hNodes = _nodesOnEdge[eh];
869 SMESHNodes bNodes = _nodesOnEdge[eb];
870 SMESHNodes gNodes = _nodesOnEdge[eg];
871 SMESHNodes dNodes = _nodesOnEdge[ed];
872 nodesOnQuad.resize( bNodes.size(), SMESHNodes(gNodes.size(), static_cast<SMDS_MeshNode*>(NULL)) );
876 // int &b_i = i, &h_i = i, &g_j = j, &d_j = j;
877 int *b_i = &i, *h_i = &i, *g_j = &j, *d_j = &j;
878 bool uWay = true, vWay = true;
880 if ( bNodes[0] != S1 ){
883 ASSERT( bNodes[bNodes.size()-1] == S1 );
885 ASSERT( bNodes[0] == S1);
887 if ( hNodes[0] != S4 ){
890 ASSERT( hNodes[0] == S4 );
892 if ( gNodes[0] != S1 ){
896 ASSERT( gNodes[0] == S1 );
898 if ( dNodes[0] != S2 ){
901 ASSERT( dNodes[0] == S2 );
906 for (i = 0, _i = bNodes.size()-1; i < bNodes.size(); ++i, --_i){
907 nodesOnQuad[i][0] = bNodes[*b_i];
908 nodesOnQuad[i][gNodes.size()-1 ] = hNodes[*h_i];
910 u = _nodeXx[ bNodes[*b_i] ];
917 if ( S1 != nodesOnQuad[0][0] ){
918 MESSAGE("ZZZZZZZZZZZZZZZZ quadID = "<<quad.getId());
920 // ASSERT( S1 == nodesOnQuad[0][0] );
924 for (j = 0, _j = gNodes.size()-1; j < gNodes.size(); ++j, --_j){
925 nodesOnQuad[0][j] = gNodes[*g_j];
926 if ( S1 != nodesOnQuad[0][0] ){
927 MESSAGE("XXXXXXXXXXXXXXXX quadID = "<<quad.getId());
929 // ASSERT( S1 == nodesOnQuad[0][0] );
930 nodesOnQuad[bNodes.size()-1][j] = dNodes[*d_j];
931 v = _nodeXx[ gNodes[*g_j] ];
940 int iSize = nodesOnQuad.size();
941 int jSize = nodesOnQuad[0].size();
942 ASSERT( iSize = bNodes.size() );
943 ASSERT( jSize = gNodes.size() );
945 // ASSERT( S1 == nodesOnQuad[0][0] );
946 // ASSERT( S2 == nodesOnQuad[iSize-1][0]);
947 // ASSERT( S4 == nodesOnQuad[0][jSize-1]);
948 // ASSERT( S3 == nodesOnQuad[iSize-1][jSize-1]);
954 bool SMESH_HexaBlocks::computeQuadByLinearApproximation( HEXA_NS::Quad& quad, bool way )
956 // int id = quad.getId();
957 // if ( quad.getId() != 66 ) return false; //7, 41
959 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
960 MESSAGE("quadID = "<<quad.getId());
962 ASSERT( _computeEdgeOK );
965 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
966 SMESHFaces facesOnQuad;
967 SMDS_MeshFace* newFace = NULL;
968 std::vector<double> xx, yy;
970 // Elements for quad computation
971 SMDS_MeshNode *S1, *S2, *S4, *S3;
973 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
974 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
975 if ( initOk == false ){
979 int iSize = nodesOnQuad.size();
980 int jSize = nodesOnQuad[0].size();
982 S1 = nodesOnQuad[0][0];
983 // S2 = nodesOnQuad[bNodes.size()-1][0];
984 S2 = nodesOnQuad[iSize-1][0];
985 S4 = nodesOnQuad[0][jSize-1];
986 S3 = nodesOnQuad[iSize-1][jSize-1];
988 for (int j = 1; j < jSize; ++j){
989 for (int i = 1; i < iSize; ++i){
990 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
991 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
992 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
993 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
996 double newNodeX, newNodeY, newNodeZ;
997 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
998 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
999 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
1000 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
1004 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
1005 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
1006 nodesOnQuad[i][j] = n4;
1009 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
1010 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
1011 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
1012 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
1015 MESSAGE("AddFace( n1, n2, n3, n4 )");
1016 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
1018 MESSAGE("AddFace( n4, n3, n2, n1 )");
1019 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
1021 facesOnQuad.push_back(newFace);
1024 _quadNodes[ &quad ] = nodesOnQuad;
1025 _facesOnQuad[&quad] = facesOnQuad;
1027 MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
1032 // --------------------------------------------------------------
1034 // --------------------------------------------------------------
1035 bool SMESH_HexaBlocks::computeHexa( HEXA_NS::Document* doc )
1037 MESSAGE("computeHexa() : : begin <<<<<<");
1040 SMESH_MesherHelper aHelper(*_theMesh);
1041 TopoDS_Shape shape = _theMesh->GetShapeToMesh();
1042 aHelper.SetSubShape( shape );
1043 aHelper.SetElementsOnShape( true );
1045 SMESH_Gen* gen = _theMesh->GetGen();
1046 SMESH_HexaFromSkin_3D algo( 0, gen, doc );
1047 algo.InitComputeError();
1049 ok = algo.Compute( *_theMesh, &aHelper, _volumesOnHexa, _node );
1051 MESSAGE("SMESH_HexaFromSkin_3D error!!! ");
1054 MESSAGE("SMESH_HexaFromSkin_3D.comment = "<<algo.GetComputeError()->myComment);
1055 MESSAGE("computeHexa() : end >>>>>>>>");
1062 // --------------------------------------------------------------
1063 // Document computing
1064 // --------------------------------------------------------------
1065 bool SMESH_HexaBlocks::computeDoc( HEXA_NS::Document* doc )
1067 MESSAGE("computeDoc() : : begin <<<<<<");
1070 // A) Vertex computation
1073 int nVertex = doc->countUsedVertex();
1074 HEXA_NS::Vertex* vertex = NULL;
1076 for (int j=0; j <nVertex; ++j ){ //Computing each vertex of the document
1077 vertex = doc->getUsedVertex(j);
1078 ok = computeVertex(*vertex);
1081 // B) Edges computation
1083 HEXA_NS::Propagation* propa = NULL;
1084 HEXA_NS::Law* law = NULL;
1085 HEXA_NS::Edges edges;
1087 nbPropa = doc->countPropagation();
1088 for (int j=0; j < nbPropa; ++j ){//Computing each edge's propagations of the document
1089 propa = doc->getPropagation(j);
1090 edges = propa->getEdges();
1091 law = propa->getLaw();
1094 law = doc->getLaw(0); // default law
1096 for( HEXA_NS::Edges::const_iterator iter = edges.begin();
1097 iter != edges.end();
1099 ok = computeEdge(**iter, *law);
1102 // C) Quad computation
1103 std::map<HEXA_NS::Quad*, bool> quadWays = computeQuadWays(doc);
1104 int nQuad = doc->countUsedQuad();
1105 HEXA_NS::Quad* q = NULL;
1106 for (int j=0; j <nQuad; ++j ){ //Computing each quad of the document
1107 q = doc->getUsedQuad(j);
1108 int id = q->getId();
1109 if ( quadWays.count(q) > 0 )
1110 ok = computeQuad( *q, quadWays[q] );
1112 MESSAGE("NO QUAD WAY ID = "<<id);
1116 // D) Hexa computation: Calling HexaFromSkin algo
1117 ok = computeHexa(doc);
1119 MESSAGE("computeDoc() : end >>>>>>>>");
1125 void SMESH_HexaBlocks::buildGroups(HEXA_NS::Document* doc)
1127 MESSAGE("_addGroups() : : begin <<<<<<");
1128 MESSAGE("_addGroups() : : nb. hexas= " << doc->countUsedHexa());
1129 MESSAGE("_addGroups() : : nb. quads= " << doc->countUsedQuad());
1130 MESSAGE("_addGroups() : : nb. edges= " << doc->countUsedEdge());
1131 MESSAGE("_addGroups() : : nb. nodes= " << doc->countUsedVertex());
1133 // Looping on each groups of the document
1134 for ( int i=0; i < doc->countGroup(); i++ ){
1135 _fillGroup( doc->getGroup(i) );
1138 MESSAGE("_addGroups() : end >>>>>>>>");
1141 // --------------------------------------------------------------
1143 // --------------------------------------------------------------
1144 double SMESH_HexaBlocks::_Xx( double i, HEXA_NS::Law law, double nbNodes) //, double pos0 )
1149 HEXA_NS::KindLaw k = law.getKind();
1150 double coeff = law.getCoefficient();
1152 case HEXA_NS::Uniform:
1153 result = (i+1)/(nbNodes+1);
1154 MESSAGE( "_Xx():" << " Uniform u("<<i<< ")"<< " = " << result);
1156 case HEXA_NS::Arithmetic:
1157 u0 = 1./(nbNodes + 1.) - (coeff*nbNodes)/2.;
1159 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Arithmetic discretization : check coefficient")));
1163 result = (i + 1.)*u0 + coeff*i*(i+1.)/2.;
1165 MESSAGE( "_Xx():" << " Arithmetic u("<<i<< ")"<< " = " << result);
1167 case HEXA_NS::Geometric:
1168 u0 = (1.-coeff)/(1.-pow(coeff, nbNodes + 1) ) ;
1170 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Geometric discretization : check coefficient")));
1174 result = u0 * (1.- pow(coeff, i + 1) )/(1.-coeff) ;
1176 MESSAGE( "_Xx():" << " Geometric u("<<i<< ")"<< " = " << result);
1183 // ============================================================= _edgeLength
1184 double SMESH_HexaBlocks::_edgeLength(const TopoDS_Edge & E)
1186 MESSAGE("_edgeLength() : : begin <<<<<<");
1187 double UMin = 0, UMax = 0;
1188 if (BRep_Tool::Degenerated(E))
1191 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
1192 GeomAdaptor_Curve AdaptCurve(C);
1193 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
1194 MESSAGE("_edgeLength() : end >>>>>>>>");
1197 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1198 // ============================================================== _buildMyCurve
1199 // === construire ma courbe a moi
1200 void SMESH_HexaBlocks::_buildMyCurve(
1201 const gp_Pnt& myCurve_pt_start, //IN
1202 const gp_Pnt& myCurve_pt_end, //IN
1203 std::list< BRepAdaptor_Curve* >& myCurve_list, //INOUT
1204 double& myCurve_tot_len, //INOUT
1205 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//INOUT
1206 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //INOUT
1207 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //INOUT
1208 HEXA_NS::Edge& edge) // For error diagnostic
1210 MESSAGE("_buildMyCurve() : : begin <<<<<<");
1211 bool current_way = true;
1212 myCurve_tot_len = 0.;
1213 BRepAdaptor_Curve* thePreviousCurve = NULL;
1214 BRepAdaptor_Curve* theCurve = NULL;
1216 gp_Pnt curv_start, curv_end;
1217 gp_Pnt p_curv_start , p_curv_end;
1219 double curv_length = 0;
1221 int nbass = edge.countAssociation();
1222 for (int nro = 0 ; nro < nbass ; ++nro)
1224 theCurve = make_curve (curv_start, curv_end, curv_length, u_start,
1226 if (theCurve != NULL)
1229 if ( thePreviousCurve == NULL )
1231 // setting current_way and first curve way
1232 if ( myCurve_pt_start.IsEqual(curv_start, HEXA_EPS) )
1237 else if ( myCurve_pt_start.IsEqual(curv_end, HEXA_EPS) )
1242 else if ( myCurve_pt_end.IsEqual(curv_end, HEXA_EPS) )
1244 current_way = false;
1247 else if ( myCurve_pt_end.IsEqual(curv_start, HEXA_EPS) )
1249 current_way = false;
1254 MESSAGE("SOMETHING WRONG on edge association... Bad script?");
1256 throw (SALOME_Exception (LOCALIZED("Edge association : check association parameters ( first, last ) between HEXA model and CAO")));
1260 // it is not the first or last curve.
1261 // ways are calculated between previous and new one.
1264 if ( p_curv_end.IsEqual( curv_end, HEXA_EPS )
1265 || p_curv_start.IsEqual( curv_start, HEXA_EPS ) )
1267 sens = NOT myCurve_ways[thePreviousCurve];// opposite WAY
1269 else if ( p_curv_end.IsEqual( curv_start, HEXA_EPS )
1270 || p_curv_start.IsEqual( curv_end, HEXA_EPS ) )
1272 sens = myCurve_ways[thePreviousCurve];// same WAY
1276 MESSAGE("SOMETHING WRONG on edge association... bad script?");
1279 throw (SALOME_Exception(LOCALIZED("Edge association : Check association parameters ( first, last ) between HEXA model and CAO")));
1283 myCurve_list.push_back( theCurve );
1284 myCurve_ways [theCurve] = sens;
1285 myCurve_starts [theCurve] = u_start;
1286 myCurve_lengths[theCurve] = curv_length;
1287 myCurve_tot_len += curv_length;
1289 thePreviousCurve = theCurve;
1290 p_curv_start = curv_start;
1291 p_curv_end = curv_end;
1292 }//if ( theCurveLength > 0 )
1296 if ( NOT current_way ){
1297 std::list< BRepAdaptor_Curve* > tmp( myCurve_list.size() );
1298 std::copy( myCurve_list.rbegin(), myCurve_list.rend(), tmp.begin() );
1302 MESSAGE("current_way was :" << current_way);
1303 MESSAGE("_buildMyCurve() : end >>>>>>>>");
1309 gp_Pnt SMESH_HexaBlocks::_getPtOnMyCurve(
1310 const double& myCurve_u, //IN
1311 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //IN
1312 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//IN
1313 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //IN
1314 std::list< BRepAdaptor_Curve* >& myCurve_list, //INOUT
1315 double& myCurve_start ) //INOUT
1316 // std::map< BRepAdaptor_Curve*, double>& myCurve_firsts,
1317 // std::map< BRepAdaptor_Curve*, double>& myCurve_lasts,
1319 MESSAGE("_getPtOnMyCurve() : : begin <<<<<<");
1322 // looking for curve which contains parameter myCurve_u
1323 BRepAdaptor_Curve* curve = myCurve_list.front();
1324 double curve_start = myCurve_start;
1325 double curve_end = curve_start + myCurve_lengths[curve];
1327 GCPnts_AbscissaPoint discret;
1329 MESSAGE("looking for curve = "<<(long) curve);
1330 MESSAGE("looking for curve: curve_u = "<<myCurve_u);
1331 MESSAGE("looking for curve: curve_start = "<<curve_start);
1332 MESSAGE("looking for curve: curve_end = "<<curve_end);
1333 MESSAGE("looking for curve: curve_lenght = "<<myCurve_lengths[curve]);
1334 MESSAGE("looking for curve: curve.size _lenght= "<<myCurve_list.size());
1336 while ( !( (myCurve_u >= curve_start) && (myCurve_u <= curve_end) ) ) {
1338 ASSERT( myCurve_list.size() != 0 );
1339 myCurve_list.pop_front();
1340 curve = myCurve_list.front();
1341 curve_start = curve_end;
1342 curve_end = curve_start + myCurve_lengths[curve];
1344 MESSAGE("go next curve: curve_lenght = "<<myCurve_lengths[curve]);
1345 MESSAGE("go next curve: curve_start = "<<curve_start);
1346 MESSAGE("go next curve: curve_end = "<<curve_end);
1347 MESSAGE("go next curve: myCurve_u = "<<myCurve_u);
1349 myCurve_start = curve_start;
1352 if ( myCurve_ways[curve] ){
1353 discret = GCPnts_AbscissaPoint( *curve, (myCurve_u - curve_start), myCurve_starts[curve] );
1355 discret = GCPnts_AbscissaPoint( *curve, myCurve_lengths[curve]- (myCurve_u - curve_start), myCurve_starts[curve] );
1357 // PutData (discret);
1358 ASSERT(discret.IsDone());
1359 curve_u = discret.Parameter();
1360 ptOnMyCurve = curve->Value( curve_u );
1362 MESSAGE("curve found!");
1363 MESSAGE("curve_u = "<< curve_u);
1364 MESSAGE("curve way = "<< myCurve_ways[curve]);
1365 MESSAGE("_getPtOnMyCurve() : end >>>>>>>>");
1373 void SMESH_HexaBlocks::_nodeInterpolationUV(double u, double v,
1374 SMDS_MeshNode* Pg, SMDS_MeshNode* Pd, SMDS_MeshNode* Ph, SMDS_MeshNode* Pb,
1375 SMDS_MeshNode* S1, SMDS_MeshNode* S2, SMDS_MeshNode* S3, SMDS_MeshNode* S4,
1376 double& xOut, double& yOut, double& zOut )
1378 MESSAGE("_nodeInterpolationUV() IN:");
1379 MESSAGE("u ( "<< u <<" )");
1380 MESSAGE("v ( "<< v <<" )");
1382 MESSAGE("S1 (" << S1->X() << "," << S1->Y() << "," << S1->Z() << ")");
1383 MESSAGE("S2 (" << S2->X() << "," << S2->Y() << "," << S2->Z() << ")");
1384 MESSAGE("S4 (" << S4->X() << "," << S4->Y() << "," << S4->Z() << ")");
1385 MESSAGE("S3 (" << S3->X() << "," << S3->Y() << "," << S3->Z() << ")");
1387 MESSAGE("Pg (" << Pg->X() << "," << Pg->Y() << "," << Pg->Z() << ")");
1388 MESSAGE("Pd (" << Pd->X() << "," << Pd->Y() << "," << Pd->Z() << ")");
1389 MESSAGE("Ph (" << Ph->X() << "," << Ph->Y() << "," << Ph->Z() << ")");
1390 MESSAGE("Pb (" << Pb->X() << "," << Pb->Y() << "," << Pb->Z() << ")");
1392 xOut = ((1.-u)*Pg->X() + v*Ph->X() + u*Pd->X() + (1.-v)*Pb->X()) - (1.-u)*(1.-v)*S1->X() - u*(1.-v)*S2->X() - u*v*S3->X() - (1.-u)*v*S4->X();
1393 yOut = ((1.-u)*Pg->Y() + v*Ph->Y() + u*Pd->Y() + (1.-v)*Pb->Y()) - (1.-u)*(1.-v)*S1->Y() - u*(1.-v)*S2->Y() - u*v*S3->Y() - (1.-u)*v*S4->Y();
1394 zOut = ((1.-u)*Pg->Z() + v*Ph->Z() + u*Pd->Z() + (1.-v)*Pb->Z()) - (1.-u)*(1.-v)*S1->Z() - u*(1.-v)*S2->Z() - u*v*S3->Z() - (1.-u)*v*S4->Z();
1396 MESSAGE("_nodeInterpolationUV() OUT("<<xOut<<","<<yOut<<","<<zOut<<" )");
1399 // =========================================================== getFaceShapes
1400 TopoDS_Shape SMESH_HexaBlocks::getFaceShapes (Hex::Quad& quad)
1402 int nbass = quad.countAssociation ();
1403 Hex::FaceShape* face = quad.getAssociation (0);
1405 return face->getShape ();
1407 TopoDS_Compound compound;
1408 BRep_Builder builder;
1409 builder.MakeCompound (compound);
1411 for (int nro=0 ; nro < nbass ; nro++)
1413 face = quad.getAssociation (nro);
1414 TopoDS_Shape shape = face->getShape ();
1415 builder.Add (compound, shape);
1422 // ================================================== carre
1423 inline double carre (double val)
1427 // ================================================== dist2
1428 inline double dist2 (const gp_Pnt& pt1, const gp_Pnt& pt2)
1430 double dist = carre (pt2.X()-pt1.X()) + carre (pt2.Y()-pt1.Y())
1431 + carre (pt2.Z()-pt1.Z());
1434 // ================================================== _intersect
1435 gp_Pnt SMESH_HexaBlocks::_intersect( const gp_Pnt& Pt,
1436 const gp_Vec& u, const gp_Vec& v,
1437 const TopoDS_Shape& shape,
1442 gp_Vec normale = u^v;
1443 gp_Dir dir(normale);
1444 gp_Lin li( Pt, dir );
1446 Standard_Real s = -Precision::Infinite();
1447 Standard_Real e = +Precision::Infinite();
1449 IntCurvesFace_ShapeIntersector inter;
1450 inter.Load(shape, tol);
1451 // inter.Load(S, tol);
1452 inter.Perform(li, s, e);//inter.PerformNearest(li, s, e);
1454 /*********************************************** Abu 2011-11-04 */
1455 if ( inter.IsDone() )
1457 result = inter.Pnt(1);//first
1458 int nbrpts = inter.NbPnt();
1461 double d0 = dist2 (result, Pt);
1462 for (int i=2; i <= inter.NbPnt(); ++i )
1464 double d1 = dist2 (Pt, inter.Pnt(i));
1468 result = inter.Pnt (i);
1472 /********************************************** Abu 2011-11-04 (getEnd()) */
1473 if (SALOME::VerbosityActivated()){
1474 MESSAGE("_intersect() : OK");
1475 for ( int i=1; i <= inter.NbPnt(); ++i ){
1476 gp_Pnt tmp = inter.Pnt(i);
1477 MESSAGE("_intersect() : pnt("<<i<<") = ("<<tmp.X()<<","<<tmp.Y()<<","<<tmp.Z()<<" )");
1482 MESSAGE("_intersect() : KO");
1491 // parameters q : IN, v0: INOUT, v1: INOUT
1492 void SMESH_HexaBlocks::_searchInitialQuadWay( HEXA_NS::Quad* q, HEXA_NS::Vertex*& v0, HEXA_NS::Vertex*& v1 )
1494 MESSAGE("_searchInitialQuadWay() : begin");
1495 v0 = NULL; v1 = NULL;
1496 if ( q->getNbrParents() != 1 ) return; // q must be a skin quad
1498 HEXA_NS::Vertex* qA = q->getVertex(0);
1499 HEXA_NS::Vertex* qB = q->getVertex(1);
1500 HEXA_NS::Vertex* qC = q->getVertex(2);
1501 HEXA_NS::Vertex* qD = q->getVertex(3);
1503 // searching for vertex on opposed quad
1504 HEXA_NS::Vertex *qAA = NULL, *qBB = NULL, *qCC = NULL, *qDD = NULL;
1505 HEXA_NS::Hexa* h = q->getParent(0);
1506 for( int i=0; i < h->countEdge(); ++i ){
1507 HEXA_NS::Edge* e = h->getEdge(i);
1508 HEXA_NS::Vertex* e0 = e->getVertex(0);
1509 HEXA_NS::Vertex* e1 = e->getVertex(1);
1511 if ( e0 == qA && e1 != qB && e1 != qC && e1 != qD ){
1513 } else if ( e1 == qA && e0 != qB && e0 != qC && e0 != qD ){
1515 } else if ( e0 == qB && e1 != qA && e1 != qC && e1 != qD ){
1517 } else if ( e1 == qB && e0 != qA && e0 != qC && e0 != qD ){
1519 } else if ( e0 == qC && e1 != qA && e1 != qB && e1 != qD ){
1521 } else if ( e1 == qC && e0 != qA && e0 != qB && e0 != qD ){
1523 } else if ( e0 == qD && e1 != qA && e1 != qB && e1 != qC ){
1525 } else if ( e1 == qD && e0 != qA && e0 != qB && e0 != qC ){
1530 // working on final value ( point on CAO ), not on model
1531 SMDS_MeshNode *nA = _node[qA], *nAA = _node[qAA];
1532 SMDS_MeshNode *nB = _node[qB], *nBB = _node[qBB];
1533 SMDS_MeshNode *nC = _node[qC], *nCC = _node[qCC];
1534 SMDS_MeshNode *nD = _node[qD], *nDD = _node[qDD];
1536 gp_Pnt pA( nA->X(), nA->Y(), nA->Z() );
1537 gp_Pnt pB( nB->X(), nB->Y(), nB->Z() );
1538 gp_Pnt pC( nC->X(), nC->Y(), nC->Z() );
1539 gp_Pnt pD( nD->X(), nD->Y(), nD->Z() );
1541 gp_Pnt pAA( nAA->X(), nAA->Y(), nAA->Z() );
1542 gp_Pnt pBB( nBB->X(), nBB->Y(), nBB->Z() );
1543 gp_Pnt pCC( nCC->X(), nCC->Y(), nCC->Z() );
1544 gp_Pnt pDD( nDD->X(), nDD->Y(), nDD->Z() );
1546 gp_Vec AB( pA, pB );
1547 gp_Vec AC( pA, pC );
1548 gp_Vec normP = AB^AC;
1549 gp_Dir dirP( normP );
1551 // building plane for point projection
1552 gp_Pln plnP( gp_Pnt(nA->X(), nA->Y(), nA->Z()), dirP);
1553 TopoDS_Shape sPlnP = BRepBuilderAPI_MakeFace(plnP).Face();
1555 // PAAA is the result of PAA projection
1556 gp_Pnt pAAA = _intersect( pAA, AB, AC, sPlnP );
1557 gp_Pnt pBBB = _intersect( pBB, AB, AC, sPlnP );
1558 gp_Pnt pCCC = _intersect( pCC, AB, AC, sPlnP );
1559 gp_Pnt pDDD = _intersect( pDD, AB, AC, sPlnP );
1561 gp_Dir AA( gp_Vec(pAA, pAAA) );
1562 gp_Dir BB( gp_Vec(pBB, pBBB) );
1563 gp_Dir CC( gp_Vec(pCC, pCCC) );
1564 gp_Dir DD( gp_Vec(pDD, pDDD) );
1566 // eventually, we are able to know if the input quad is a good client!
1567 // exit the fonction otherwise
1568 if ( AA.IsOpposite(BB, HEXA_QUAD_WAY) ) return;
1569 if ( BB.IsOpposite(CC, HEXA_QUAD_WAY) ) return;
1570 if ( CC.IsOpposite(DD, HEXA_QUAD_WAY) ) return;
1572 // ok, give the input quad the good orientation by
1574 if ( !dirP.IsOpposite(AA, HEXA_QUAD_WAY) ) { //OK
1580 MESSAGE("_searchInitialQuadWay() : end");
1583 SMESH_Group* SMESH_HexaBlocks::_createGroup(HEXA_NS::Group& grHex)
1585 MESSAGE("_createGroup() : : begin <<<<<<");
1587 std::string aGrName = grHex.getName();
1588 HEXA_NS::EnumGroup grHexKind = grHex.getKind();
1590 MESSAGE("aGrName"<<aGrName);
1592 SMDSAbs_ElementType aGrType;
1593 switch ( grHexKind ){
1594 case HEXA_NS::HexaCell : aGrType = SMDSAbs_Volume; break;
1595 case HEXA_NS::QuadCell : aGrType = SMDSAbs_Face ; break;
1596 case HEXA_NS::EdgeCell : aGrType = SMDSAbs_Edge ; break;
1597 case HEXA_NS::HexaNode : aGrType = SMDSAbs_Node ; break;
1598 case HEXA_NS::QuadNode : aGrType = SMDSAbs_Node ; break;
1599 case HEXA_NS::EdgeNode : aGrType = SMDSAbs_Node ; break;
1600 case HEXA_NS::VertexNode : aGrType = SMDSAbs_Node ; break;
1601 default : ASSERT(false);
1604 SMESH_Group* aGr = _theMesh->AddGroup(aGrType, aGrName.c_str());
1606 MESSAGE("_createGroup() : end >>>>>>>>");
1610 void SMESH_HexaBlocks::_fillGroup(HEXA_NS::Group* grHex )
1612 MESSAGE("_fillGroup() : : begin <<<<<<");
1614 SMESH_Group* aGr = _createGroup( *grHex );
1615 HEXA_NS::EltBase* grHexElt = NULL;
1616 HEXA_NS::EnumGroup grHexKind = grHex->getKind();
1617 int grHexNbElt = grHex->countElement();
1619 MESSAGE("_fillGroup() : kind = " << grHexKind);
1620 MESSAGE("_fillGroup() : count= " << grHexNbElt);
1622 // A)Looking for elements ID
1623 std::vector<const SMDS_MeshElement*> aGrEltIDs;
1625 for ( int n=0; n<grHexNbElt; ++n ){
1626 grHexElt = grHex->getElement(n);
1628 switch ( grHexKind ){
1629 case HEXA_NS::HexaCell:
1631 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1632 if ( _volumesOnHexa.count(h)>0 ){
1633 SMESHVolumes volumes = _volumesOnHexa[h];
1634 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1635 aGrEltIDs.push_back(*aVolume);
1638 MESSAGE("GROUP OF VOLUME: volume for hexa (id = "<<h->getId()<<") not found");
1642 case HEXA_NS::QuadCell:
1644 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1645 if ( _facesOnQuad.count(q)>0 ){
1646 SMESHFaces faces = _facesOnQuad[q];
1647 for ( SMESHFaces::iterator aFace = faces.begin(); aFace != faces.end(); ++aFace ){
1648 aGrEltIDs.push_back(*aFace);
1651 MESSAGE("GROUP OF FACE: face for quad (id = "<<q->getId()<<") not found");
1655 case HEXA_NS::EdgeCell:
1657 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1658 if ( _edgesOnEdge.count(e)>0 ){
1659 SMESHEdges edges = _edgesOnEdge[e];
1660 for ( SMESHEdges::iterator anEdge = edges.begin(); anEdge != edges.end(); ++anEdge ){
1661 aGrEltIDs.push_back(*anEdge);
1664 MESSAGE("GROUP OF Edge: edge for edge (id = "<<e->getId()<<") not found");
1668 case HEXA_NS::HexaNode:
1670 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1671 if ( _volumesOnHexa.count(h)>0 ){
1672 SMESHVolumes volumes = _volumesOnHexa[h];
1673 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1674 SMDS_ElemIteratorPtr aNodeIter = (*aVolume)->nodesIterator();
1675 while( aNodeIter->more() ){
1676 const SMDS_MeshNode* aNode =
1677 dynamic_cast<const SMDS_MeshNode*>( aNodeIter->next() );
1679 aGrEltIDs.push_back(aNode);
1684 MESSAGE("GROUP OF HEXA NODES: nodes on hexa (id = "<<h->getId()<<") not found");
1688 case HEXA_NS::QuadNode:
1690 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1691 if ( _quadNodes.count(q)>0 ){
1692 ArrayOfSMESHNodes nodesOnQuad = _quadNodes[q];
1693 for ( ArrayOfSMESHNodes::iterator nodes = nodesOnQuad.begin(); nodes != nodesOnQuad.end(); ++nodes){
1694 for ( SMESHNodes::iterator aNode = nodes->begin(); aNode != nodes->end(); ++aNode){
1695 aGrEltIDs.push_back(*aNode);
1699 MESSAGE("GROUP OF QUAD NODES: nodes on quad (id = "<<q->getId()<<") not found");
1703 case HEXA_NS::EdgeNode:
1705 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1706 if ( _nodesOnEdge.count(e)>0 ){
1707 SMESHNodes nodes = _nodesOnEdge[e];
1708 for ( SMESHNodes::iterator aNode = nodes.begin(); aNode != nodes.end(); ++aNode){
1709 aGrEltIDs.push_back(*aNode);
1712 MESSAGE("GROUP OF EDGE NODES: nodes on edge (id = "<<e->getId()<<") not found");
1716 case HEXA_NS::VertexNode:
1718 HEXA_NS::Vertex* v = reinterpret_cast<HEXA_NS::Vertex*>(grHexElt);
1719 if ( _node.count(v)>0 ){
1720 aGrEltIDs.push_back(_node[v]);
1722 MESSAGE("GROUP OF VERTEX NODES: nodes for vertex (id = "<<v->getId()<<") not found");
1726 default : ASSERT(false);
1730 // B)Filling the group on SMESH
1731 SMESHDS_Group* aGroupDS = dynamic_cast<SMESHDS_Group*>( aGr->GetGroupDS() );
1733 for ( int i=0; i < aGrEltIDs.size(); i++ ) {
1734 aGroupDS->SMDSGroup().Add( aGrEltIDs[i] );
1737 MESSAGE("_fillGroup() : end >>>>>>>>");