1 // Copyright (C) 2009-2019 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 int MYDEBUG = HEXA_NS::on_debug ();
94 static int MYDEBUG = 0;
98 static double HEXA_EPS = 1.0e-3; //1E-3;
99 static double HEXA_QUAD_WAY = M_PI/4.; //3.*PI/8.;
101 // ============================================================ string2shape
102 TopoDS_Shape string2shape( const std::string& brep )
105 std::istringstream streamBrep(brep);
106 BRep_Builder aBuilder;
107 BRepTools::Read(shape, streamBrep, aBuilder);
111 // ============================================================ shape2coord
112 bool shape2coord(const TopoDS_Shape& aShape, double& x, double& y, double& z)
114 if ( aShape.ShapeType() == TopAbs_VERTEX ){
115 TopoDS_Vertex aPoint;
116 aPoint = TopoDS::Vertex( aShape );
117 gp_Pnt aPnt = BRep_Tool::Pnt( aPoint );
126 // ============================================================= edge_length
127 double edge_length(const TopoDS_Edge & E)
129 double UMin = 0, UMax = 0;
130 if (BRep_Tool::Degenerated(E))
133 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
134 GeomAdaptor_Curve AdaptCurve(C);
135 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
138 // =============================================================== make_curve
139 BRepAdaptor_Curve* make_curve (gp_Pnt& pstart, gp_Pnt& pend,
140 double& length, double& u_start,
141 HEXA_NS::Edge& edge, int nro)
143 Hex::AssoEdge* asso = edge.getAssociation (nro);
144 BRepAdaptor_Curve* the_curve = asso->getCurve();
145 length = asso->length ();
149 const double* point1 = asso->getOrigin ();
150 const double* point2 = asso->getExtrem ();
151 pstart = gp_Pnt (point1[0], point1[1], point1[2]);
152 pend = gp_Pnt (point2[0], point2[1], point2[2]);
153 u_start = asso->getUstart();
156 // ============================================================== Constructeur
157 // SMESH_HexaBlocks::SMESH_HexaBlocks( SMESH_Mesh* theMesh ):
158 SMESH_HexaBlocks::SMESH_HexaBlocks(SMESH_Mesh& theMesh):
162 _computeVertexOK(false),
163 _computeEdgeOK(false),
164 _computeQuadOK(false),
165 _theMesh(&theMesh), //groups creation
166 _theMeshDS(theMesh.GetMeshDS()) //meshing
171 SMESH_HexaBlocks::~SMESH_HexaBlocks()
176 // --------------------------------------------------------------
178 // --------------------------------------------------------------
180 // --------------------------------------------------------------
182 // --------------------------------------------------------------
183 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
184 // ============================================================== computeVertex
185 bool SMESH_HexaBlocks::computeVertex(HEXA_NS::Vertex& vx)
188 vx.getAssoCoord (px, py, pz);
190 SMDS_MeshNode* new_node = _theMeshDS->AddNode (px, py, pz);
191 _vertex [new_node] = &vx;
192 _node [&vx] = new_node; //needed in computeEdge()
193 _computeVertexOK = true;
196 // --------------------------------------------------------------
198 // --------------------------------------------------------------
199 bool SMESH_HexaBlocks::computeEdge(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
203 ok = computeEdgeByAssoc( edge, law);
205 ok = computeEdgeByShortestWire( edge, law);
208 ok = computeEdgeByPlanWire( edge, law);
211 ok = computeEdgeByIsoWire( edge, law);
214 ok = computeEdgeBySegment( edge, law);
217 _computeEdgeOK = true;
224 bool SMESH_HexaBlocks::computeEdgeByAssoc( HEXA_NS::Edge& edge, HEXA_NS::Law& law )
226 if(MYDEBUG) MESSAGE("computeEdgeByAssoc(edgeID = "<<edge.getId()<<"): begin <<<<<<");
227 ASSERT( _computeVertexOK );
230 if (NOT edge.isAssociated())
234 HEXA_NS::Vertex* vx0 = NULL;
235 HEXA_NS::Vertex* vx1 = NULL;
237 // way of discretization
238 if (edge.getWay() == true){
239 vx0 = edge.getVertex(0);
240 vx1 = edge.getVertex(1);
242 vx0 = edge.getVertex(1);
243 vx1 = edge.getVertex(0);
246 SMDS_MeshNode* FIRST_NODE = _node[vx0];
247 SMDS_MeshNode* LAST_NODE = _node[vx1];
250 // A) Build myCurve_list
251 std::list< BRepAdaptor_Curve* > myCurve_list;
252 double myCurve_tot_len;
253 std::map< BRepAdaptor_Curve*, double> myCurve_lengths;
254 std::map< BRepAdaptor_Curve*, bool> myCurve_ways;
255 std::map< BRepAdaptor_Curve*, double> myCurve_starts;
257 gp_Pnt myCurve_pt_start(FIRST_NODE->X(), FIRST_NODE->Y(), FIRST_NODE->Z());
258 gp_Pnt myCurve_pt_end (LAST_NODE->X(), LAST_NODE->Y(), LAST_NODE->Z());
272 // B) Build nodes and edges on mesh from myCurve_list
273 SMDS_MeshNode* node_a = NULL;
274 SMDS_MeshNode* node_b = NULL;
275 SMDS_MeshEdge* edge_ab = NULL;
276 SMESHNodes nodesOnEdge;
277 SMESHEdges edgesOnEdge; //backup for group creation
281 nodesOnEdge.push_back(FIRST_NODE);
282 // nodesXxOnEdge.push_back(0.);
283 // _nodeXx[FIRST_NODE] = 0.;
287 double myCurve_start_u = 0.;
288 int nbNodes = law.getNodes(); //law of discretization
289 if (myCurve_list.size()==0)
291 PutData (edge.getName());
294 if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
295 for (int i = 0; i < nbNodes; ++i){
296 u = _Xx(i, law, nbNodes); //u between [0,1]
297 myCurve_u = u*myCurve_tot_len;
300 MESSAGE("myCurve_u -> "<<myCurve_u);
301 MESSAGE("myCurve_tot_len -> "<<myCurve_tot_len);
303 ptOnMyCurve = _getPtOnMyCurve( myCurve_u,
311 node_b = _theMeshDS->AddNode( ptOnMyCurve.X(), ptOnMyCurve.Y(), ptOnMyCurve.Z() );
312 edge_ab = _theMeshDS->AddEdge( node_a, node_b );
313 nodesOnEdge.push_back( node_b );
314 edgesOnEdge.push_back( edge_ab );
315 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
319 edge_ab = _theMeshDS->AddEdge( node_a, LAST_NODE );
320 nodesOnEdge.push_back( LAST_NODE );
321 edgesOnEdge.push_back( edge_ab );
322 _nodesOnEdge[&edge] = nodesOnEdge;
323 _edgesOnEdge[&edge] = edgesOnEdge;
326 if(MYDEBUG) MESSAGE("computeEdgeByAssoc() : end >>>>>>>>");
333 bool SMESH_HexaBlocks::computeEdgeByShortestWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
335 if(MYDEBUG) MESSAGE("computeEdgeByShortestWire() not implemented");
339 bool SMESH_HexaBlocks::computeEdgeByPlanWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
341 if(MYDEBUG) MESSAGE("computeEdgeByPlanWire() not implemented");
345 bool SMESH_HexaBlocks::computeEdgeByIsoWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
347 if(MYDEBUG) MESSAGE("computeEdgeByIsoWire() not implemented");
352 bool SMESH_HexaBlocks::computeEdgeBySegment(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
354 if(MYDEBUG) MESSAGE("computeEdgeBySegment() : : begin <<<<<<");
355 ASSERT( _computeVertexOK );
359 HEXA_NS::Vertex* vx0 = NULL;
360 HEXA_NS::Vertex* vx1 = NULL;
362 // way of discretization
363 if (edge.getWay() == true){
364 vx0 = edge.getVertex(0);
365 vx1 = edge.getVertex(1);
367 vx0 = edge.getVertex(1);
368 vx1 = edge.getVertex(0);
372 SMDS_MeshNode* FIRST_NODE = _node[vx0];
373 SMDS_MeshNode* LAST_NODE = _node[vx1];
374 SMDS_MeshNode* node_a = NULL; //new node (to be added)
375 SMDS_MeshNode* node_b = NULL; //new node (to be added)
377 // node and edge creation
378 SMESHNodes nodesOnEdge;
379 SMESHEdges edgesOnEdge;
382 double newNodeX, newNodeY, newNodeZ;
383 SMDS_MeshEdge* newEdge = NULL;
386 nodesOnEdge.push_back(FIRST_NODE);
388 //law of discretization
389 int nbNodes = law.getNodes();
390 if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
391 for (int i = 0; i < nbNodes; ++i){
392 u = _Xx(i, law, nbNodes);
393 newNodeX = FIRST_NODE->X() + u * ( LAST_NODE->X() - FIRST_NODE->X() );
394 newNodeY = FIRST_NODE->Y() + u * ( LAST_NODE->Y() - FIRST_NODE->Y() );
395 newNodeZ = FIRST_NODE->Z() + u * ( LAST_NODE->Z() - FIRST_NODE->Z() );
396 node_b = _theMeshDS->AddNode(newNodeX, newNodeY, newNodeZ);
397 newEdge = _theMeshDS->AddEdge(node_a, node_b);
398 edgesOnEdge.push_back(newEdge);
399 nodesOnEdge.push_back(node_b);
400 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
401 _nodeXx[ node_b ] = u;
402 if(MYDEBUG) MESSAGE("_nodeXx <-"<<u);
405 newEdge = _theMeshDS->AddEdge(node_a, LAST_NODE);
406 nodesOnEdge.push_back(LAST_NODE);
407 edgesOnEdge.push_back(newEdge);
409 _nodesOnEdge[&edge] = nodesOnEdge;
410 _edgesOnEdge[&edge] = edgesOnEdge;
412 if(MYDEBUG) MESSAGE("computeEdgeBySegment() : end >>>>>>>>");
417 // --------------------------------------------------------------
419 // --------------------------------------------------------------
420 std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document* doc )
422 std::map<HEXA_NS::Quad*, bool> quadWays;
423 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
424 std::list<HEXA_NS::Quad*> skinQuad;
425 std::list<HEXA_NS::Quad*> workingQuad;
426 HEXA_NS::Quad* first_q = NULL;
427 HEXA_NS::Quad* q = NULL;
428 HEXA_NS::Edge* e = NULL;
429 HEXA_NS::Vertex *e_0, *e_1 = NULL;
431 // FIRST STEP: eliminate free quad + internal quad
432 int nTotalQuad = doc->countUsedQuad();
433 for (int i=0; i < nTotalQuad; ++i ){
434 q = doc->getUsedQuad(i);
435 switch ( q->getNbrParents() ){ // parent == hexaedron
436 case 0: case 2: quadWays[q] = true; break;
437 case 1: skinQuad.push_back(q); break;
438 default: if ( q->getNbrParents() > 2 ) ASSERT(false);
442 // SECOND STEP: setting edges ways
443 while ( skinQuad.size()>0 ){
444 if(MYDEBUG) MESSAGE("SEARCHING INITIAL QUAD ..." );
445 for ( std::list<HEXA_NS::Quad*>::iterator it = skinQuad.begin(); it != skinQuad.end(); it++ ){
446 _searchInitialQuadWay( *it, e_0, e_1 );
447 if ( e_0 != NULL && e_1 != NULL ){
452 if ( e_0 == NULL && e_1 == NULL ) ASSERT(false);// should never happened,
453 if(MYDEBUG) MESSAGE("INITIAL QUAD FOUND!" );
454 for ( int j=0 ; j < 4 ; ++j ){
456 if ( ((e_0 == e->getVertex(0)) && (e_1 == e->getVertex(1)))
457 || ((e_0 == e->getVertex(1)) && (e_1 == e->getVertex(0))) ){
461 if(MYDEBUG) MESSAGE("INITIAL EDGE WAY FOUND!" );
463 edgeWays[e] = std::make_pair( e_0, e_1 );
464 workingQuad.push_back(q);
466 while ( workingQuad.size() > 0 ){
467 if(MYDEBUG) MESSAGE("COMPUTE QUAD WAY ... ID ="<< q->getId());
468 HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
470 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
471 while ( localEdgeWays.size() != 4 ){
472 HEXA_NS::Edge* e = q->getEdge(i%4);
473 if ( lastVertex == NULL ){
474 if ( edgeWays.count(e) == 1 ){
476 localEdgeWays[e] = std::make_pair( edgeWays[e].first, edgeWays[e].second );
478 localEdgeWays[e] = std::make_pair( edgeWays[e].second, edgeWays[e].first);
480 firstVertex = localEdgeWays[e].first;
481 lastVertex = localEdgeWays[e].second;
484 HEXA_NS::Vertex* e_0 = e->getVertex(0);
485 HEXA_NS::Vertex* e_1 = e->getVertex(1);
487 if ( lastVertex == e_0 ){
488 firstVertex = e_0; lastVertex = e_1;
489 } else if ( lastVertex == e_1 ){
490 firstVertex = e_1; lastVertex = e_0;
491 } else if ( firstVertex == e_0 ) {
492 firstVertex = e_1; lastVertex = e_0;
493 } else if ( firstVertex == e_1 ) {
494 firstVertex = e_0; lastVertex = e_1;
498 localEdgeWays[e] = std::make_pair( firstVertex, lastVertex );
499 if ( edgeWays.count(e) == 0 ){ // keep current value if present otherwise add it
500 edgeWays[e] = localEdgeWays[e];
507 //add new working quad
508 for ( int i=0 ; i < 4; ++i ){
509 HEXA_NS::Quad* next_q = NULL;
510 HEXA_NS::Edge* e = q->getEdge(i);
511 for ( int j=0 ; j < e->getNbrParents(); ++j ){
512 next_q = e->getParent(j);
516 int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
518 int fromWorkingQuad = std::count( workingQuad.begin(), workingQuad.end(), next_q );
519 if ( fromWorkingQuad == 0 ){
520 workingQuad.push_front( next_q );
527 HEXA_NS::Edge* e0 = q->getEdge(0);
528 HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
530 if ( e0_0 == localEdgeWays[ e0 ].first ){
532 } else if ( e0_0 == localEdgeWays[ e0 ].second ){
537 workingQuad.remove( q );
538 skinQuad.remove( q );
539 q = workingQuad.back();
550 // std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document& doc, std::map<HEXA_NS::Quad*, bool> initQuads )
552 // std::map<HEXA_NS::Quad*, bool> quadWays;
553 // // std::map<HEXA_NS::Edge*, bool> edgeWays;
554 // // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
555 // std::map<HEXA_NS::Quad*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > workingQuads;
557 // std::list<HEXA_NS::Quad*> skinQuad;
558 // std::list<HEXA_NS::Quad*> notSkinQuad;
559 // // std::list<HEXA_NS::Quad*> workingQuad;
560 // HEXA_NS::Quad* first_q = NULL;
561 // HEXA_NS::Quad* q = NULL;
562 // HEXA_NS::Edge* e = NULL;
563 // HEXA_NS::Vertex *e_0, *e_1 = NULL;
565 // // FIRST STEP: eliminate free quad + internal quad
566 // int nTotalQuad = doc.countQuad();
567 // for (int i=0; i < nTotalQuad; ++i ){
568 // q = doc.getQuad(i);
569 // switch ( q->getNbrParents() ){ // parent == hexaedron
570 // case 0: case 2: quadWays[q] = true; break;
571 // // case 0: case 2: notSkinQuad.push_back(q); break; //CS_TEST
572 // case 1: skinQuad.push_back(q); break;
573 // default: if ( q->getNbrParents() > 2 ) ASSERT(false);
579 // q = first_q = skinQuad.front();
580 // e = q->getUsedEdge(0);
581 // e_0 = e->getVertex(0);
582 // e_1 = e->getVertex(1);
584 // workingQuads[q] = std::make_pair( e_0, e_1 );
586 // while ( workingQuads.size() > 0 ){
587 // MESSAGE("CURRENT QUAD ------>"<< q->getId());
588 // HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
591 // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
592 // while ( localEdgeWays.size() != 4 ){
593 // HEXA_NS::Edge* e = q->getUsedEdge(i%4);
594 // if ( lastVertex == NULL ){
595 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
596 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
598 // if ( (workingQuads[q].first == e_0 and workingQuads[q].second == e_1)
599 // or (workingQuads[q].first == e_1 and workingQuads[q].second == e_0) ){
600 // if ( q == first_q ){
601 // localEdgeWays[e] = std::make_pair( workingQuads[q].first, workingQuads[q].second );
602 // MESSAGE("FIRST QUAD ");
604 // localEdgeWays[e] = std::make_pair( workingQuads[q].second, workingQuads[q].first);
605 // MESSAGE("NOT FIRST QUAD ");
607 // firstVertex = localEdgeWays[e].first;
608 // lastVertex = localEdgeWays[e].second;
611 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
612 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
613 // if ( lastVertex == e_0 ){
614 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
615 // firstVertex = e_0;
617 // } else if ( lastVertex == e_1 ){
618 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
619 // firstVertex = e_1;
621 // } else if ( firstVertex == e_0 ) {
622 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
623 // firstVertex = e_1;
625 // } else if ( firstVertex == e_1 ) {
626 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
627 // firstVertex = e_0;
637 // //add new working quad
638 // for ( int i=0 ; i < 4; ++i ){
639 // HEXA_NS::Quad* next_q = NULL;
640 // HEXA_NS::Edge* e = q->getUsedEdge(i);
641 // MESSAGE("NB PARENTS ="<< e->getNbrParents() );
642 // for ( int j=0 ; j < e->getNbrParents(); ++j ){
643 // next_q = e->getParent(j);
644 // if ( next_q == q ){
647 // int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
648 // if (fromSkin != 0){
649 // // int fromWorkingQuad = std::count( workingQuads.begin(), workingQuads.end(), next_q );
650 // int fromWorkingQuad = workingQuads.count( next_q );
651 // // MESSAGE("CHECK QUAD:"<< newWorkingQuad->getId());
652 // if ( fromWorkingQuad == 0 ){
653 // // workingQuads.push_front( next_q );
654 // workingQuads[ next_q ] = localEdgeWays[e];
655 // // MESSAGE("EDGE :"<<e->getId()<<"ADD QUAD :"<< newWorkingQuad->getId());
661 // //setting quad way
662 // HEXA_NS::Edge* e0 = q->getUsedEdge(0);
663 // HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
665 // if ( e0_0 == localEdgeWays[ e0 ].first ){
666 // quadWays[q] = true;
667 // } else if ( e0_0 == localEdgeWays[ e0 ].second ){
668 // quadWays[q] = false;
672 // MESSAGE("quadWays ID ="<< q->getId() << ", WAY = " << quadWays[q] );
674 // // workingQuad.remove( q );
675 // workingQuads.erase( q );
676 // skinQuad.remove( q );
677 // *workingQuads.begin();
678 // q = (*workingQuads.begin()).first;
684 bool SMESH_HexaBlocks::computeQuad( HEXA_NS::Quad& quad, bool way )
688 ok = computeQuadByAssoc(quad, way);
690 ok = computeQuadByFindingGeom(quad, way);
693 ok = computeQuadByLinearApproximation(quad, way);
696 _computeQuadOK = true;
702 bool SMESH_HexaBlocks::computeQuadByAssoc( HEXA_NS::Quad& quad, bool way )
704 // int id = quad.getId();
705 // if ( id != 11 ) return false; //7
707 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
708 MESSAGE("quadID = "<<quad.getId());
710 ASSERT( _computeEdgeOK );
713 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
714 SMESHFaces facesOnQuad;
715 SMDS_MeshFace* newFace = NULL;
716 std::vector<double> xx, yy;
718 // Elements for quad computation
719 SMDS_MeshNode *S1, *S2, *S4, *S3;
721 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
722 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
723 if ( initOk == false ){
727 int nbass = quad.countAssociation ();
730 if(MYDEBUG) MESSAGE("computeQuadByAssoc() : end >>>>>>>>");
734 TopoDS_Shape shapeOrCompound = getFaceShapes ( quad );
737 std::map<SMDS_MeshNode*, gp_Pnt> interpolatedPoints;
738 int iSize = nodesOnQuad.size();
739 int jSize = nodesOnQuad[0].size();
741 S1 = nodesOnQuad[0][0];
742 S2 = nodesOnQuad[iSize-1][0];
743 S4 = nodesOnQuad[0][jSize-1];
744 S3 = nodesOnQuad[iSize-1][jSize-1];
747 for (int j = 1; j < jSize; ++j){
748 for (int i = 1; i < iSize; ++i){
749 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
750 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
751 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
752 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
755 double newNodeX, newNodeY, newNodeZ;
756 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
757 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
758 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
759 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
763 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
764 gp_Pnt newPt = gp_Pnt( newNodeX, newNodeY, newNodeZ );//interpolated point
767 if ( interpolatedPoints.count(n1) > 0 ){
768 pt1 = interpolatedPoints[n1];
770 pt1 = gp_Pnt( n1->X(), n1->Y(), n1->Z() );
772 if ( interpolatedPoints.count(n3) > 0 ){
773 pt3 = interpolatedPoints[n3];
775 pt3 = gp_Pnt( n3->X(), n3->Y(), n3->Z() );
777 gp_Vec vec1( newPt, pt1 );
778 gp_Vec vec2( newPt, pt3 );
780 gp_Pnt ptOnShape = _intersect(newPt, vec1, vec2, shapeOrCompound);
781 newNodeX = ptOnShape.X();
782 newNodeY = ptOnShape.Y();
783 newNodeZ = ptOnShape.Z();
784 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
785 nodesOnQuad[i][j] = n4;
786 interpolatedPoints[ n4 ] = newPt;
789 MESSAGE("u parameter is "<<u);
790 MESSAGE("v parameter is "<<v);
791 MESSAGE("point interpolated ("<<newPt.X()<<","<<newPt.Y()<<","<<newPt.Z()<<" )");
792 MESSAGE("point on shape ("<<newNodeX<<","<<newNodeY<<","<<newNodeZ<<" )");
797 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
798 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
799 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
800 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
804 if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
805 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
807 if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
808 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
810 facesOnQuad.push_back(newFace);
813 _quadNodes[ &quad ] = nodesOnQuad;
814 _facesOnQuad[&quad] = facesOnQuad;
816 if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
821 bool SMESH_HexaBlocks::computeQuadByFindingGeom( HEXA_NS::Quad& quad, bool way )
823 if(MYDEBUG) MESSAGE("computeQuadByFindingGeom() not implemented");
827 bool SMESH_HexaBlocks::_computeQuadInit(
829 ArrayOfSMESHNodes& nodesOnQuad,
830 std::vector<double>& xx, std::vector<double>& yy)
832 if(MYDEBUG) MESSAGE("_computeQuadInit() : begin ---------------");
835 SMDS_MeshNode *S1, *S2, *S4, *S3;
836 HEXA_NS::Edge *eh, *eb, *eg, *ed;
837 HEXA_NS::Edge *e1, *e2, *e3, *e4;
838 HEXA_NS::Vertex *e1_0, *e1_1, *e2_0, *e2_1, *e3_0, *e3_1, *e4_0, *e4_1;
840 e1 = quad.getEdge(0);
841 e2 = quad.getEdge(1);
842 e3 = quad.getEdge(2);
843 e4 = quad.getEdge(3);
845 e1_0 = e1->getVertex(0); e1_1 = e1->getVertex(1);
846 e2_0 = e2->getVertex(0); e2_1 = e2->getVertex(1);
847 e3_0 = e3->getVertex(0); e3_1 = e3->getVertex(1);
848 e4_0 = e4->getVertex(0); e4_1 = e4->getVertex(1);
851 S1 = _node[e1_0]; S2 = _node[e1_1];
857 } else if ( e1_0 == e2_1 ){
860 } else if ( e1_0 == e4_0 ){
863 } else if ( e1_0 == e4_1 ){
870 if ( S4 == _node[e3_0] ){
872 } else if ( S4 == _node[e3_1] ){
878 SMESHNodes hNodes = _nodesOnEdge[eh];
879 SMESHNodes bNodes = _nodesOnEdge[eb];
880 SMESHNodes gNodes = _nodesOnEdge[eg];
881 SMESHNodes dNodes = _nodesOnEdge[ed];
882 nodesOnQuad.resize( bNodes.size(), SMESHNodes(gNodes.size(), static_cast<SMDS_MeshNode*>(NULL)) );
886 // int &b_i = i, &h_i = i, &g_j = j, &d_j = j;
887 int *b_i = &i, *h_i = &i, *g_j = &j, *d_j = &j;
888 bool uWay = true, vWay = true;
890 if ( bNodes[0] != S1 ){
893 ASSERT( bNodes[bNodes.size()-1] == S1 );
895 ASSERT( bNodes[0] == S1);
897 if ( hNodes[0] != S4 ){
900 ASSERT( hNodes[0] == S4 );
902 if ( gNodes[0] != S1 ){
906 ASSERT( gNodes[0] == S1 );
908 if ( dNodes[0] != S2 ){
911 ASSERT( dNodes[0] == S2 );
916 for (i = 0, _i = bNodes.size()-1; i < bNodes.size(); ++i, --_i){
917 nodesOnQuad[i][0] = bNodes[*b_i];
918 nodesOnQuad[i][gNodes.size()-1 ] = hNodes[*h_i];
920 u = _nodeXx[ bNodes[*b_i] ];
927 if ( S1 != nodesOnQuad[0][0] ){
928 if(MYDEBUG) MESSAGE("ZZZZZZZZZZZZZZZZ quadID = "<<quad.getId());
930 // ASSERT( S1 == nodesOnQuad[0][0] );
934 for (j = 0, _j = gNodes.size()-1; j < gNodes.size(); ++j, --_j){
935 nodesOnQuad[0][j] = gNodes[*g_j];
936 if ( S1 != nodesOnQuad[0][0] ){
937 if(MYDEBUG) MESSAGE("XXXXXXXXXXXXXXXX quadID = "<<quad.getId());
939 // ASSERT( S1 == nodesOnQuad[0][0] );
940 nodesOnQuad[bNodes.size()-1][j] = dNodes[*d_j];
941 v = _nodeXx[ gNodes[*g_j] ];
950 int iSize = nodesOnQuad.size();
951 int jSize = nodesOnQuad[0].size();
952 ASSERT( iSize = bNodes.size() );
953 ASSERT( jSize = gNodes.size() );
955 // ASSERT( S1 == nodesOnQuad[0][0] );
956 // ASSERT( S2 == nodesOnQuad[iSize-1][0]);
957 // ASSERT( S4 == nodesOnQuad[0][jSize-1]);
958 // ASSERT( S3 == nodesOnQuad[iSize-1][jSize-1]);
964 bool SMESH_HexaBlocks::computeQuadByLinearApproximation( HEXA_NS::Quad& quad, bool way )
966 // int id = quad.getId();
967 // if ( quad.getId() != 66 ) return false; //7, 41
969 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
970 MESSAGE("quadID = "<<quad.getId());
972 ASSERT( _computeEdgeOK );
975 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
976 SMESHFaces facesOnQuad;
977 SMDS_MeshFace* newFace = NULL;
978 std::vector<double> xx, yy;
980 // Elements for quad computation
981 SMDS_MeshNode *S1, *S2, *S4, *S3;
983 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
984 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
985 if ( initOk == false ){
989 int iSize = nodesOnQuad.size();
990 int jSize = nodesOnQuad[0].size();
992 S1 = nodesOnQuad[0][0];
993 // S2 = nodesOnQuad[bNodes.size()-1][0];
994 S2 = nodesOnQuad[iSize-1][0];
995 S4 = nodesOnQuad[0][jSize-1];
996 S3 = nodesOnQuad[iSize-1][jSize-1];
998 for (int j = 1; j < jSize; ++j){
999 for (int i = 1; i < iSize; ++i){
1000 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
1001 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
1002 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
1003 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
1006 double newNodeX, newNodeY, newNodeZ;
1007 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
1008 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
1009 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
1010 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
1014 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
1015 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
1016 nodesOnQuad[i][j] = n4;
1020 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
1021 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
1022 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
1023 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
1027 if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
1028 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
1030 if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
1031 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
1033 facesOnQuad.push_back(newFace);
1036 _quadNodes[ &quad ] = nodesOnQuad;
1037 _facesOnQuad[&quad] = facesOnQuad;
1039 if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
1044 // --------------------------------------------------------------
1046 // --------------------------------------------------------------
1047 bool SMESH_HexaBlocks::computeHexa( HEXA_NS::Document* doc )
1049 if(MYDEBUG) MESSAGE("computeHexa() : : begin <<<<<<");
1052 SMESH_MesherHelper aHelper(*_theMesh);
1053 TopoDS_Shape shape = _theMesh->GetShapeToMesh();
1054 aHelper.SetSubShape( shape );
1055 aHelper.SetElementsOnShape( true );
1057 SMESH_Gen* gen = _theMesh->GetGen();
1058 SMESH_HexaFromSkin_3D algo( 0, gen, doc );
1059 algo.InitComputeError();
1061 ok = algo.Compute( *_theMesh, &aHelper, _volumesOnHexa, _node );
1063 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D error!!! ");
1066 MESSAGE("SMESH_HexaFromSkin_3D.comment = "<<algo.GetComputeError()->myComment);
1067 MESSAGE("computeHexa() : end >>>>>>>>");
1074 // --------------------------------------------------------------
1075 // Document computing
1076 // --------------------------------------------------------------
1077 bool SMESH_HexaBlocks::computeDoc( HEXA_NS::Document* doc )
1079 if(MYDEBUG) MESSAGE("computeDoc() : : begin <<<<<<");
1082 // A) Vertex computation
1085 int nVertex = doc->countUsedVertex();
1086 HEXA_NS::Vertex* vertex = NULL;
1088 for (int j=0; j <nVertex; ++j ){ //Computing each vertex of the document
1089 vertex = doc->getUsedVertex(j);
1090 ok = computeVertex(*vertex);
1093 // B) Edges computation
1095 HEXA_NS::Propagation* propa = NULL;
1096 HEXA_NS::Law* law = NULL;
1097 HEXA_NS::Edges edges;
1099 nbPropa = doc->countPropagation();
1100 for (int j=0; j < nbPropa; ++j ){//Computing each edge's propagations of the document
1101 propa = doc->getPropagation(j);
1102 edges = propa->getEdges();
1103 law = propa->getLaw();
1106 law = doc->getLaw(0); // default law
1108 for( HEXA_NS::Edges::const_iterator iter = edges.begin();
1109 iter != edges.end();
1111 ok = computeEdge(**iter, *law);
1114 // C) Quad computation
1115 std::map<HEXA_NS::Quad*, bool> quadWays = computeQuadWays(doc);
1116 int nQuad = doc->countUsedQuad();
1117 HEXA_NS::Quad* q = NULL;
1118 for (int j=0; j <nQuad; ++j ){ //Computing each quad of the document
1119 q = doc->getUsedQuad(j);
1120 int id = q->getId();
1121 if ( quadWays.count(q) > 0 )
1122 ok = computeQuad( *q, quadWays[q] );
1124 if(MYDEBUG) MESSAGE("NO QUAD WAY ID = "<<id);
1128 // D) Hexa computation: Calling HexaFromSkin algo
1129 ok = computeHexa(doc);
1131 if(MYDEBUG) MESSAGE("computeDoc() : end >>>>>>>>");
1137 void SMESH_HexaBlocks::buildGroups(HEXA_NS::Document* doc)
1140 MESSAGE("_addGroups() : : begin <<<<<<");
1141 MESSAGE("_addGroups() : : nb. hexas= " << doc->countUsedHexa());
1142 MESSAGE("_addGroups() : : nb. quads= " << doc->countUsedQuad());
1143 MESSAGE("_addGroups() : : nb. edges= " << doc->countUsedEdge());
1144 MESSAGE("_addGroups() : : nb. nodes= " << doc->countUsedVertex());
1146 // Looping on each groups of the document
1147 for ( int i=0; i < doc->countGroup(); i++ ){
1148 _fillGroup( doc->getGroup(i) );
1151 if(MYDEBUG) MESSAGE("_addGroups() : end >>>>>>>>");
1154 // --------------------------------------------------------------
1156 // --------------------------------------------------------------
1157 double SMESH_HexaBlocks::_Xx( double i, HEXA_NS::Law law, double nbNodes) //, double pos0 )
1162 HEXA_NS::KindLaw k = law.getKind();
1163 double coeff = law.getCoefficient();
1165 case HEXA_NS::Uniform:
1166 result = (i+1)/(nbNodes+1);
1167 if(MYDEBUG) MESSAGE( "_Xx():" << " Uniform u("<<i<< ")"<< " = " << result);
1169 case HEXA_NS::Arithmetic:
1170 u0 = 1./(nbNodes + 1.) - (coeff*nbNodes)/2.;
1172 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Arithmetic discretization : check coefficient")));
1176 result = (i + 1.)*u0 + coeff*i*(i+1.)/2.;
1178 if(MYDEBUG) MESSAGE( "_Xx():" << " Arithmetic u("<<i<< ")"<< " = " << result);
1180 case HEXA_NS::Geometric:
1181 u0 = (1.-coeff)/(1.-pow(coeff, nbNodes + 1) ) ;
1183 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Geometric discretization : check coefficient")));
1187 result = u0 * (1.- pow(coeff, i + 1) )/(1.-coeff) ;
1189 if(MYDEBUG) MESSAGE( "_Xx():" << " Geometric u("<<i<< ")"<< " = " << result);
1196 // ============================================================= _edgeLength
1197 double SMESH_HexaBlocks::_edgeLength(const TopoDS_Edge & E)
1199 if(MYDEBUG) MESSAGE("_edgeLength() : : begin <<<<<<");
1200 double UMin = 0, UMax = 0;
1201 if (BRep_Tool::Degenerated(E))
1204 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
1205 GeomAdaptor_Curve AdaptCurve(C);
1206 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
1207 if(MYDEBUG) MESSAGE("_edgeLength() : end >>>>>>>>");
1210 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1211 // ============================================================== _buildMyCurve
1212 // === construire ma courbe a moi
1213 void SMESH_HexaBlocks::_buildMyCurve(
1214 const gp_Pnt& myCurve_pt_start, //IN
1215 const gp_Pnt& myCurve_pt_end, //IN
1216 std::list< BRepAdaptor_Curve* >& myCurve_list, //INOUT
1217 double& myCurve_tot_len, //INOUT
1218 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//INOUT
1219 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //INOUT
1220 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //INOUT
1221 HEXA_NS::Edge& edge) // For error diagnostic
1223 if(MYDEBUG) MESSAGE("_buildMyCurve() : : begin <<<<<<");
1224 bool current_way = true;
1225 myCurve_tot_len = 0.;
1226 BRepAdaptor_Curve* thePreviousCurve = NULL;
1227 BRepAdaptor_Curve* theCurve = NULL;
1229 gp_Pnt curv_start, curv_end;
1230 gp_Pnt p_curv_start , p_curv_end;
1232 double curv_length = 0;
1234 int nbass = edge.countAssociation();
1235 for (int nro = 0 ; nro < nbass ; ++nro)
1237 theCurve = make_curve (curv_start, curv_end, curv_length, u_start,
1239 if (theCurve != NULL)
1242 if ( thePreviousCurve == NULL )
1244 // setting current_way and first curve way
1245 if ( myCurve_pt_start.IsEqual(curv_start, HEXA_EPS) )
1250 else if ( myCurve_pt_start.IsEqual(curv_end, HEXA_EPS) )
1255 else if ( myCurve_pt_end.IsEqual(curv_end, HEXA_EPS) )
1257 current_way = false;
1260 else if ( myCurve_pt_end.IsEqual(curv_start, HEXA_EPS) )
1262 current_way = false;
1268 MESSAGE("SOMETHING WRONG on edge association... Bad script?");
1270 throw (SALOME_Exception (LOCALIZED("Edge association : check association parameters ( first, last ) between HEXA model and CAO")));
1274 // it is not the first or last curve.
1275 // ways are calculated between previous and new one.
1278 if ( p_curv_end.IsEqual( curv_end, HEXA_EPS )
1279 || p_curv_start.IsEqual( curv_start, HEXA_EPS ) )
1281 sens = NOT myCurve_ways[thePreviousCurve];// opposite WAY
1283 else if ( p_curv_end.IsEqual( curv_start, HEXA_EPS )
1284 || p_curv_start.IsEqual( curv_end, HEXA_EPS ) )
1286 sens = myCurve_ways[thePreviousCurve];// same WAY
1291 MESSAGE("SOMETHING WRONG on edge association... bad script?");
1294 throw (SALOME_Exception(LOCALIZED("Edge association : Check association parameters ( first, last ) between HEXA model and CAO")));
1298 myCurve_list.push_back( theCurve );
1299 myCurve_ways [theCurve] = sens;
1300 myCurve_starts [theCurve] = u_start;
1301 myCurve_lengths[theCurve] = curv_length;
1302 myCurve_tot_len += curv_length;
1304 thePreviousCurve = theCurve;
1305 p_curv_start = curv_start;
1306 p_curv_end = curv_end;
1307 }//if ( theCurveLength > 0 )
1311 if ( NOT current_way ){
1312 std::list< BRepAdaptor_Curve* > tmp( myCurve_list.size() );
1313 std::copy( myCurve_list.rbegin(), myCurve_list.rend(), tmp.begin() );
1318 MESSAGE("current_way was :" << current_way);
1319 MESSAGE("_buildMyCurve() : end >>>>>>>>");
1326 gp_Pnt SMESH_HexaBlocks::_getPtOnMyCurve(
1327 const double& myCurve_u, //IN
1328 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //IN
1329 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//IN
1330 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //IN
1331 std::list< BRepAdaptor_Curve* >& myCurve_list, //INOUT
1332 double& myCurve_start ) //INOUT
1333 // std::map< BRepAdaptor_Curve*, double>& myCurve_firsts,
1334 // std::map< BRepAdaptor_Curve*, double>& myCurve_lasts,
1336 if(MYDEBUG) MESSAGE("_getPtOnMyCurve() : : begin <<<<<<");
1339 // looking for curve which contains parameter myCurve_u
1340 BRepAdaptor_Curve* curve = myCurve_list.front();
1341 double curve_start = myCurve_start;
1342 double curve_end = curve_start + myCurve_lengths[curve];
1344 GCPnts_AbscissaPoint discret;
1347 MESSAGE("looking for curve = "<<(long) curve);
1348 MESSAGE("looking for curve: curve_u = "<<myCurve_u);
1349 MESSAGE("looking for curve: curve_start = "<<curve_start);
1350 MESSAGE("looking for curve: curve_end = "<<curve_end);
1351 MESSAGE("looking for curve: curve_lenght = "<<myCurve_lengths[curve]);
1352 MESSAGE("looking for curve: curve.size _lenght= "<<myCurve_list.size());
1354 while ( !( (myCurve_u >= curve_start) && (myCurve_u <= curve_end) ) ) {
1356 ASSERT( myCurve_list.size() != 0 );
1357 myCurve_list.pop_front();
1358 curve = myCurve_list.front();
1359 curve_start = curve_end;
1360 curve_end = curve_start + myCurve_lengths[curve];
1362 MESSAGE("go next curve: curve_lenght = "<<myCurve_lengths[curve]);
1363 MESSAGE("go next curve: curve_start = "<<curve_start);
1364 MESSAGE("go next curve: curve_end = "<<curve_end);
1365 MESSAGE("go next curve: myCurve_u = "<<myCurve_u);
1368 myCurve_start = curve_start;
1371 if ( myCurve_ways[curve] ){
1372 discret = GCPnts_AbscissaPoint( *curve, (myCurve_u - curve_start), myCurve_starts[curve] );
1374 discret = GCPnts_AbscissaPoint( *curve, myCurve_lengths[curve]- (myCurve_u - curve_start), myCurve_starts[curve] );
1376 // PutData (discret);
1377 ASSERT(discret.IsDone());
1378 curve_u = discret.Parameter();
1379 ptOnMyCurve = curve->Value( curve_u );
1382 MESSAGE("curve found!");
1383 MESSAGE("curve_u = "<< curve_u);
1384 MESSAGE("curve way = "<< myCurve_ways[curve]);
1385 MESSAGE("_getPtOnMyCurve() : end >>>>>>>>");
1393 void SMESH_HexaBlocks::_nodeInterpolationUV(double u, double v,
1394 SMDS_MeshNode* Pg, SMDS_MeshNode* Pd, SMDS_MeshNode* Ph, SMDS_MeshNode* Pb,
1395 SMDS_MeshNode* S1, SMDS_MeshNode* S2, SMDS_MeshNode* S3, SMDS_MeshNode* S4,
1396 double& xOut, double& yOut, double& zOut )
1399 MESSAGE("_nodeInterpolationUV() IN:");
1400 MESSAGE("u ( "<< u <<" )");
1401 MESSAGE("v ( "<< v <<" )");
1403 MESSAGE("S1 (" << S1->X() << "," << S1->Y() << "," << S1->Z() << ")");
1404 MESSAGE("S2 (" << S2->X() << "," << S2->Y() << "," << S2->Z() << ")");
1405 MESSAGE("S4 (" << S4->X() << "," << S4->Y() << "," << S4->Z() << ")");
1406 MESSAGE("S3 (" << S3->X() << "," << S3->Y() << "," << S3->Z() << ")");
1408 MESSAGE("Pg (" << Pg->X() << "," << Pg->Y() << "," << Pg->Z() << ")");
1409 MESSAGE("Pd (" << Pd->X() << "," << Pd->Y() << "," << Pd->Z() << ")");
1410 MESSAGE("Ph (" << Ph->X() << "," << Ph->Y() << "," << Ph->Z() << ")");
1411 MESSAGE("Pb (" << Pb->X() << "," << Pb->Y() << "," << Pb->Z() << ")");
1414 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();
1415 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();
1416 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();
1419 MESSAGE("_nodeInterpolationUV() OUT("<<xOut<<","<<yOut<<","<<zOut<<" )");
1423 // =========================================================== getFaceShapes
1424 TopoDS_Shape SMESH_HexaBlocks::getFaceShapes (Hex::Quad& quad)
1426 int nbass = quad.countAssociation ();
1427 Hex::FaceShape* face = quad.getAssociation (0);
1429 return face->getShape ();
1431 TopoDS_Compound compound;
1432 BRep_Builder builder;
1433 builder.MakeCompound (compound);
1435 for (int nro=0 ; nro < nbass ; nro++)
1437 face = quad.getAssociation (nro);
1438 TopoDS_Shape shape = face->getShape ();
1439 builder.Add (compound, shape);
1446 // ================================================== carre
1447 inline double carre (double val)
1451 // ================================================== dist2
1452 inline double dist2 (const gp_Pnt& pt1, const gp_Pnt& pt2)
1454 double dist = carre (pt2.X()-pt1.X()) + carre (pt2.Y()-pt1.Y())
1455 + carre (pt2.Z()-pt1.Z());
1458 // ================================================== _intersect
1459 gp_Pnt SMESH_HexaBlocks::_intersect( const gp_Pnt& Pt,
1460 const gp_Vec& u, const gp_Vec& v,
1461 const TopoDS_Shape& shape,
1466 gp_Vec normale = u^v;
1467 gp_Dir dir(normale);
1468 gp_Lin li( Pt, dir );
1470 Standard_Real s = -Precision::Infinite();
1471 Standard_Real e = +Precision::Infinite();
1473 IntCurvesFace_ShapeIntersector inter;
1474 inter.Load(shape, tol);
1475 // inter.Load(S, tol);
1476 inter.Perform(li, s, e);//inter.PerformNearest(li, s, e);
1478 /*********************************************** Abu 2011-11-04 */
1479 if ( inter.IsDone() )
1481 result = inter.Pnt(1);//first
1482 int nbrpts = inter.NbPnt();
1485 double d0 = dist2 (result, Pt);
1486 for (int i=2; i <= inter.NbPnt(); ++i )
1488 double d1 = dist2 (Pt, inter.Pnt(i));
1492 result = inter.Pnt (i);
1496 /********************************************** Abu 2011-11-04 (getEnd()) */
1498 MESSAGE("_intersect() : OK");
1499 for ( int i=1; i <= inter.NbPnt(); ++i ){
1500 gp_Pnt tmp = inter.Pnt(i);
1501 MESSAGE("_intersect() : pnt("<<i<<") = ("<<tmp.X()<<","<<tmp.Y()<<","<<tmp.Z()<<" )");
1506 if(MYDEBUG) MESSAGE("_intersect() : KO");
1515 // parameters q : IN, v0: INOUT, v1: INOUT
1516 void SMESH_HexaBlocks::_searchInitialQuadWay( HEXA_NS::Quad* q, HEXA_NS::Vertex*& v0, HEXA_NS::Vertex*& v1 )
1518 if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : begin");
1519 v0 = NULL; v1 = NULL;
1520 if ( q->getNbrParents() != 1 ) return; // q must be a skin quad
1522 HEXA_NS::Vertex* qA = q->getVertex(0);
1523 HEXA_NS::Vertex* qB = q->getVertex(1);
1524 HEXA_NS::Vertex* qC = q->getVertex(2);
1525 HEXA_NS::Vertex* qD = q->getVertex(3);
1527 // searching for vertex on opposed quad
1528 HEXA_NS::Vertex *qAA = NULL, *qBB = NULL, *qCC = NULL, *qDD = NULL;
1529 HEXA_NS::Hexa* h = q->getParent(0);
1530 for( int i=0; i < h->countEdge(); ++i ){
1531 HEXA_NS::Edge* e = h->getEdge(i);
1532 HEXA_NS::Vertex* e0 = e->getVertex(0);
1533 HEXA_NS::Vertex* e1 = e->getVertex(1);
1535 if ( e0 == qA && e1 != qB && e1 != qC && e1 != qD ){
1537 } else if ( e1 == qA && e0 != qB && e0 != qC && e0 != qD ){
1539 } else if ( e0 == qB && e1 != qA && e1 != qC && e1 != qD ){
1541 } else if ( e1 == qB && e0 != qA && e0 != qC && e0 != qD ){
1543 } else if ( e0 == qC && e1 != qA && e1 != qB && e1 != qD ){
1545 } else if ( e1 == qC && e0 != qA && e0 != qB && e0 != qD ){
1547 } else if ( e0 == qD && e1 != qA && e1 != qB && e1 != qC ){
1549 } else if ( e1 == qD && e0 != qA && e0 != qB && e0 != qC ){
1554 // working on final value ( point on CAO ), not on model
1555 SMDS_MeshNode *nA = _node[qA], *nAA = _node[qAA];
1556 SMDS_MeshNode *nB = _node[qB], *nBB = _node[qBB];
1557 SMDS_MeshNode *nC = _node[qC], *nCC = _node[qCC];
1558 SMDS_MeshNode *nD = _node[qD], *nDD = _node[qDD];
1560 gp_Pnt pA( nA->X(), nA->Y(), nA->Z() );
1561 gp_Pnt pB( nB->X(), nB->Y(), nB->Z() );
1562 gp_Pnt pC( nC->X(), nC->Y(), nC->Z() );
1563 gp_Pnt pD( nD->X(), nD->Y(), nD->Z() );
1565 gp_Pnt pAA( nAA->X(), nAA->Y(), nAA->Z() );
1566 gp_Pnt pBB( nBB->X(), nBB->Y(), nBB->Z() );
1567 gp_Pnt pCC( nCC->X(), nCC->Y(), nCC->Z() );
1568 gp_Pnt pDD( nDD->X(), nDD->Y(), nDD->Z() );
1570 gp_Vec AB( pA, pB );
1571 gp_Vec AC( pA, pC );
1572 gp_Vec normP = AB^AC;
1573 gp_Dir dirP( normP );
1575 // building plane for point projection
1576 gp_Pln plnP( gp_Pnt(nA->X(), nA->Y(), nA->Z()), dirP);
1577 TopoDS_Shape sPlnP = BRepBuilderAPI_MakeFace(plnP).Face();
1579 // PAAA is the result of PAA projection
1580 gp_Pnt pAAA = _intersect( pAA, AB, AC, sPlnP );
1581 gp_Pnt pBBB = _intersect( pBB, AB, AC, sPlnP );
1582 gp_Pnt pCCC = _intersect( pCC, AB, AC, sPlnP );
1583 gp_Pnt pDDD = _intersect( pDD, AB, AC, sPlnP );
1585 gp_Dir AA( gp_Vec(pAA, pAAA) );
1586 gp_Dir BB( gp_Vec(pBB, pBBB) );
1587 gp_Dir CC( gp_Vec(pCC, pCCC) );
1588 gp_Dir DD( gp_Vec(pDD, pDDD) );
1590 // eventually, we are able to know if the input quad is a good client!
1591 // exit the fonction otherwise
1592 if ( AA.IsOpposite(BB, HEXA_QUAD_WAY) ) return;
1593 if ( BB.IsOpposite(CC, HEXA_QUAD_WAY) ) return;
1594 if ( CC.IsOpposite(DD, HEXA_QUAD_WAY) ) return;
1596 // ok, give the input quad the good orientation by
1598 if ( !dirP.IsOpposite(AA, HEXA_QUAD_WAY) ) { //OK
1604 if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : end");
1607 SMESH_Group* SMESH_HexaBlocks::_createGroup(HEXA_NS::Group& grHex)
1609 if(MYDEBUG) MESSAGE("_createGroup() : : begin <<<<<<");
1611 std::string aGrName = grHex.getName();
1612 HEXA_NS::EnumGroup grHexKind = grHex.getKind();
1614 if(MYDEBUG) MESSAGE("aGrName"<<aGrName);
1616 SMDSAbs_ElementType aGrType;
1617 switch ( grHexKind ){
1618 case HEXA_NS::HexaCell : aGrType = SMDSAbs_Volume; break;
1619 case HEXA_NS::QuadCell : aGrType = SMDSAbs_Face ; break;
1620 case HEXA_NS::EdgeCell : aGrType = SMDSAbs_Edge ; break;
1621 case HEXA_NS::HexaNode : aGrType = SMDSAbs_Node ; break;
1622 case HEXA_NS::QuadNode : aGrType = SMDSAbs_Node ; break;
1623 case HEXA_NS::EdgeNode : aGrType = SMDSAbs_Node ; break;
1624 case HEXA_NS::VertexNode : aGrType = SMDSAbs_Node ; break;
1625 default : ASSERT(false);
1628 SMESH_Group* aGr = _theMesh->AddGroup(aGrType, aGrName.c_str());
1630 if(MYDEBUG) MESSAGE("_createGroup() : end >>>>>>>>");
1634 void SMESH_HexaBlocks::_fillGroup(HEXA_NS::Group* grHex )
1636 if(MYDEBUG) MESSAGE("_fillGroup() : : begin <<<<<<");
1638 SMESH_Group* aGr = _createGroup( *grHex );
1639 HEXA_NS::EltBase* grHexElt = NULL;
1640 HEXA_NS::EnumGroup grHexKind = grHex->getKind();
1641 int grHexNbElt = grHex->countElement();
1643 if(MYDEBUG) MESSAGE("_fillGroup() : kind = " << grHexKind);
1644 if(MYDEBUG) MESSAGE("_fillGroup() : count= " << grHexNbElt);
1646 // A)Looking for elements ID
1647 std::vector<const SMDS_MeshElement*> aGrEltIDs;
1649 for ( int n=0; n<grHexNbElt; ++n ){
1650 grHexElt = grHex->getElement(n);
1652 switch ( grHexKind ){
1653 case HEXA_NS::HexaCell:
1655 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1656 if ( _volumesOnHexa.count(h)>0 ){
1657 SMESHVolumes volumes = _volumesOnHexa[h];
1658 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1659 aGrEltIDs.push_back(*aVolume);
1662 if(MYDEBUG) MESSAGE("GROUP OF VOLUME: volume for hexa (id = "<<h->getId()<<") not found");
1666 case HEXA_NS::QuadCell:
1668 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1669 if ( _facesOnQuad.count(q)>0 ){
1670 SMESHFaces faces = _facesOnQuad[q];
1671 for ( SMESHFaces::iterator aFace = faces.begin(); aFace != faces.end(); ++aFace ){
1672 aGrEltIDs.push_back(*aFace);
1675 if(MYDEBUG) MESSAGE("GROUP OF FACE: face for quad (id = "<<q->getId()<<") not found");
1679 case HEXA_NS::EdgeCell:
1681 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1682 if ( _edgesOnEdge.count(e)>0 ){
1683 SMESHEdges edges = _edgesOnEdge[e];
1684 for ( SMESHEdges::iterator anEdge = edges.begin(); anEdge != edges.end(); ++anEdge ){
1685 aGrEltIDs.push_back(*anEdge);
1688 if(MYDEBUG) MESSAGE("GROUP OF Edge: edge for edge (id = "<<e->getId()<<") not found");
1692 case HEXA_NS::HexaNode:
1694 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1695 if ( _volumesOnHexa.count(h)>0 ){
1696 SMESHVolumes volumes = _volumesOnHexa[h];
1697 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1698 SMDS_ElemIteratorPtr aNodeIter = (*aVolume)->nodesIterator();
1699 while( aNodeIter->more() ){
1700 const SMDS_MeshNode* aNode =
1701 dynamic_cast<const SMDS_MeshNode*>( aNodeIter->next() );
1703 aGrEltIDs.push_back(aNode);
1708 if(MYDEBUG) MESSAGE("GROUP OF HEXA NODES: nodes on hexa (id = "<<h->getId()<<") not found");
1712 case HEXA_NS::QuadNode:
1714 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1715 if ( _quadNodes.count(q)>0 ){
1716 ArrayOfSMESHNodes nodesOnQuad = _quadNodes[q];
1717 for ( ArrayOfSMESHNodes::iterator nodes = nodesOnQuad.begin(); nodes != nodesOnQuad.end(); ++nodes){
1718 for ( SMESHNodes::iterator aNode = nodes->begin(); aNode != nodes->end(); ++aNode){
1719 aGrEltIDs.push_back(*aNode);
1723 if(MYDEBUG) MESSAGE("GROUP OF QUAD NODES: nodes on quad (id = "<<q->getId()<<") not found");
1727 case HEXA_NS::EdgeNode:
1729 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1730 if ( _nodesOnEdge.count(e)>0 ){
1731 SMESHNodes nodes = _nodesOnEdge[e];
1732 for ( SMESHNodes::iterator aNode = nodes.begin(); aNode != nodes.end(); ++aNode){
1733 aGrEltIDs.push_back(*aNode);
1736 if(MYDEBUG) MESSAGE("GROUP OF EDGE NODES: nodes on edge (id = "<<e->getId()<<") not found");
1740 case HEXA_NS::VertexNode:
1742 HEXA_NS::Vertex* v = reinterpret_cast<HEXA_NS::Vertex*>(grHexElt);
1743 if ( _node.count(v)>0 ){
1744 aGrEltIDs.push_back(_node[v]);
1746 if(MYDEBUG) MESSAGE("GROUP OF VERTEX NODES: nodes for vertex (id = "<<v->getId()<<") not found");
1750 default : ASSERT(false);
1754 // B)Filling the group on SMESH
1755 SMESHDS_Group* aGroupDS = dynamic_cast<SMESHDS_Group*>( aGr->GetGroupDS() );
1757 for ( int i=0; i < aGrEltIDs.size(); i++ ) {
1758 aGroupDS->SMDSGroup().Add( aGrEltIDs[i] );
1761 if(MYDEBUG) MESSAGE("_fillGroup() : end >>>>>>>>");