1 // Copyright (C) 2009-2012 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.
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
28 #include <AIS_Shape.hxx>
30 #include <Precision.hxx>
31 #include <BRep_Tool.hxx>
32 #include <BRepTools.hxx>
33 #include <BRep_Builder.hxx>
34 #include <BRepAdaptor_Curve.hxx>
35 #include <BRepBuilderAPI_MakeFace.hxx>
37 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
38 #include <GCPnts_AbscissaPoint.hxx>
39 #include <TopoDS_Wire.hxx>
42 #include <TopoDS_Shape.hxx>
43 #include <TopoDS_Compound.hxx>
48 #include <IntCurvesFace_ShapeIntersector.hxx>
51 #include "SMDS_MeshNode.hxx"
52 #include "SMDS_MeshVolume.hxx"
53 #include "SMESH_Gen.hxx"
54 #include "SMESH_MesherHelper.hxx"
55 #include "SMESHDS_Group.hxx"
58 #include "HexDocument.hxx"
59 #include "HexVertex.hxx"
60 #include "HexEdge.hxx"
61 #include "HexQuad.hxx"
62 #include "HexHexa.hxx"
63 #include "HexPropagation.hxx"
64 #include "HexShape.hxx"
65 #include "HexGroup.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_EPSILON = 1E-3; //1E-3;
99 static double HEXA_QUAD_WAY = PI/4.; //3.*PI/8.;
100 // static double HEXA_QUAD_WAY2 = 499999999.*PI/1000000000.;
103 TopoDS_Shape string2shape( const string& brep )
106 istringstream streamBrep(brep);
107 BRep_Builder aBuilder;
108 BRepTools::Read(shape, streamBrep, aBuilder);
113 bool shape2coord(const TopoDS_Shape& aShape, double& x, double& y, double& z)
115 if ( aShape.ShapeType() == TopAbs_VERTEX ){
116 TopoDS_Vertex aPoint;
117 aPoint = TopoDS::Vertex( aShape );
118 gp_Pnt aPnt = BRep_Tool::Pnt( aPoint );
130 // SMESH_HexaBlocks::SMESH_HexaBlocks( SMESH_Mesh* theMesh ):
131 SMESH_HexaBlocks::SMESH_HexaBlocks(SMESH_Mesh& theMesh):
135 _computeVertexOK(false),
136 _computeEdgeOK(false),
137 _computeQuadOK(false),
138 _theMesh(&theMesh), //groups creation
139 _theMeshDS(theMesh.GetMeshDS()) //meshing
144 SMESH_HexaBlocks::~SMESH_HexaBlocks()
149 // --------------------------------------------------------------
151 // --------------------------------------------------------------
153 // --------------------------------------------------------------
155 // --------------------------------------------------------------
156 bool SMESH_HexaBlocks::computeVertex(HEXA_NS::Vertex& vx)
159 ok = computeVertexByAssoc( vx );
161 ok = computeVertexByModel( vx );
164 _computeVertexOK = true;
170 bool SMESH_HexaBlocks::computeVertexByAssoc(HEXA_NS::Vertex& vx)
172 if(MYDEBUG) MESSAGE("computeVertexByAssoc() : : begin <<<<<<");
175 SMDS_MeshNode* newNode = NULL; // new node on mesh
176 double x, y, z; //new node coordinates
178 HEXA_NS::Shape* assoc = vx.getAssociation();
179 if ( assoc == NULL ){
181 MESSAGE("computeVertexByAssoc() : ASSOC not found " << vx.getName ());
187 string strBrep = assoc->getBrep();
188 TopoDS_Shape shape = string2shape( strBrep );
189 ok = shape2coord( shape, x, y, z );
191 if (!ok) throw (SALOME_Exception(LOCALIZED("vertex association : shape2coord() error ")));
192 newNode = _theMeshDS->AddNode(x, y, z);
193 if (_node.count(&vx) >= 1 and MYDEBUG) MESSAGE("_node : ALREADY");
194 _node[&vx] = newNode;//needed in computeEdge()
195 _vertex[newNode] = &vx;
198 MESSAGE("computeVertexByAssoc() : ASSOC found " << vx.getName());
200 MESSAGE("( "<< x <<","<< y <<","<< z <<" )");
203 if(MYDEBUG) MESSAGE("computeVertexByAssoc() : end >>>>>>>>");
207 bool SMESH_HexaBlocks::computeVertexByModel(HEXA_NS::Vertex& vx)
209 if(MYDEBUG) MESSAGE("computeVertexByModel() : : begin <<<<<<");
212 SMDS_MeshNode* newNode = NULL; // new node on mesh
213 double x, y, z; //new node coordinates
216 // std::cout << std::endl;
221 newNode = _theMeshDS->AddNode(x, y, z);
223 if (_node.count(&vx) >= 1 and MYDEBUG) MESSAGE("_node : ALREADY");
224 _node[&vx] = newNode;//needed in computeEdge()
225 _vertex[newNode] = &vx;
227 MESSAGE("computeVertexByModel() :" << vx.getName());
229 MESSAGE("( "<< x <<","<< y <<","<< z <<" )");
232 if(MYDEBUG) MESSAGE("computeVertexByModel() : end >>>>>>>>");
236 // --------------------------------------------------------------
238 // --------------------------------------------------------------
239 bool SMESH_HexaBlocks::computeEdge(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
243 ok = computeEdgeByAssoc( edge, law);
245 ok = computeEdgeByShortestWire( edge, law);
248 ok = computeEdgeByPlanWire( edge, law);
251 ok = computeEdgeByIsoWire( edge, law);
254 ok = computeEdgeBySegment( edge, law);
257 _computeEdgeOK = true;
264 bool SMESH_HexaBlocks::computeEdgeByAssoc( HEXA_NS::Edge& edge, HEXA_NS::Law& law )
266 if(MYDEBUG) MESSAGE("computeEdgeByAssoc(edgeID = "<<edge.getId()<<"): begin <<<<<<");
267 ASSERT( _computeVertexOK );
270 const std::vector <HEXA_NS::Shape*> associations = edge.getAssociations();
271 if ( associations.size() == 0 ){
275 HEXA_NS::Vertex* vx0 = NULL;
276 HEXA_NS::Vertex* vx1 = NULL;
278 // way of discretization
279 if (edge.getWay() == true){
280 vx0 = edge.getVertex(0);
281 vx1 = edge.getVertex(1);
283 vx0 = edge.getVertex(1);
284 vx1 = edge.getVertex(0);
287 SMDS_MeshNode* FIRST_NODE = _node[vx0];
288 SMDS_MeshNode* LAST_NODE = _node[vx1];
292 std::list< BRepAdaptor_Curve* > myCurve;
293 double myCurve_length;
294 std::map< BRepAdaptor_Curve*, double> myCurve_lengths;
295 std::map< BRepAdaptor_Curve*, bool> myCurve_ways;
296 std::map< BRepAdaptor_Curve*, double> myCurve_starts;
297 gp_Pnt myCurve_start( FIRST_NODE->X(), FIRST_NODE->Y(), FIRST_NODE->Z() );
298 gp_Pnt myCurve_end( LAST_NODE->X(), LAST_NODE->Y(), LAST_NODE->Z() );
314 // B) Build nodes and edges on mesh from myCurve
315 SMDS_MeshNode* node_a = NULL;
316 SMDS_MeshNode* node_b = NULL;
317 SMDS_MeshEdge* edge_ab = NULL;
318 SMESHNodes nodesOnEdge;
319 SMESHEdges edgesOnEdge; //backup for group creation
323 nodesOnEdge.push_back(FIRST_NODE);
324 // nodesXxOnEdge.push_back(0.);
325 // _nodeXx[FIRST_NODE] = 0.;
329 double myCurve_start_u = 0.;
330 int nbNodes = law.getNodes(); //law of discretization
331 if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
332 for (int i = 0; i < nbNodes; ++i){
333 u = _Xx(i, law, nbNodes); //u between [0,1]
334 myCurve_u = u*myCurve_length;
337 MESSAGE("myCurve_u -> "<<myCurve_u);
338 MESSAGE("myCurve_length -> "<<myCurve_length);
340 ptOnMyCurve = _getPtOnMyCurve( myCurve_u,
348 node_b = _theMeshDS->AddNode( ptOnMyCurve.X(), ptOnMyCurve.Y(), ptOnMyCurve.Z() );
349 edge_ab = _theMeshDS->AddEdge( node_a, node_b );
350 nodesOnEdge.push_back( node_b );
351 edgesOnEdge.push_back( edge_ab );
352 // nodesXxOnEdge.push_back( u );
353 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
357 edge_ab = _theMeshDS->AddEdge( node_a, LAST_NODE );
358 nodesOnEdge.push_back( LAST_NODE );
359 edgesOnEdge.push_back( edge_ab );
360 // nodesXxOnEdge.push_back( 1. );
361 // _nodeXx[LAST_NODE] = 1.;
362 _nodesOnEdge[&edge] = nodesOnEdge;
363 _edgesOnEdge[&edge] = edgesOnEdge;
367 // _edgeXx[&edge] = nodesXxOnEdge;
369 if(MYDEBUG) MESSAGE("computeEdgeByAssoc() : end >>>>>>>>");
376 bool SMESH_HexaBlocks::computeEdgeByShortestWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
378 if(MYDEBUG) MESSAGE("computeEdgeByShortestWire() not implemented");
382 bool SMESH_HexaBlocks::computeEdgeByPlanWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
384 if(MYDEBUG) MESSAGE("computeEdgeByPlanWire() not implemented");
388 bool SMESH_HexaBlocks::computeEdgeByIsoWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
390 if(MYDEBUG) MESSAGE("computeEdgeByIsoWire() not implemented");
395 bool SMESH_HexaBlocks::computeEdgeBySegment(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
397 if(MYDEBUG) MESSAGE("computeEdgeBySegment() : : begin <<<<<<");
398 ASSERT( _computeVertexOK );
402 HEXA_NS::Vertex* vx0 = NULL;
403 HEXA_NS::Vertex* vx1 = NULL;
405 // way of discretization
406 if (edge.getWay() == true){
407 vx0 = edge.getVertex(0);
408 vx1 = edge.getVertex(1);
410 vx0 = edge.getVertex(1);
411 vx1 = edge.getVertex(0);
415 SMDS_MeshNode* FIRST_NODE = _node[vx0];
416 SMDS_MeshNode* LAST_NODE = _node[vx1];
417 SMDS_MeshNode* node_a = NULL; //new node (to be added)
418 SMDS_MeshNode* node_b = NULL; //new node (to be added)
420 // node and edge creation
421 SMESHNodes nodesOnEdge;
422 SMESHEdges edgesOnEdge;
425 double newNodeX, newNodeY, newNodeZ;
426 SMDS_MeshEdge* newEdge = NULL;
429 nodesOnEdge.push_back(FIRST_NODE);
431 //law of discretization
432 int nbNodes = law.getNodes();
433 if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
434 for (int i = 0; i < nbNodes; ++i){
435 u = _Xx(i, law, nbNodes);
436 newNodeX = FIRST_NODE->X() + u * ( LAST_NODE->X() - FIRST_NODE->X() );
437 newNodeY = FIRST_NODE->Y() + u * ( LAST_NODE->Y() - FIRST_NODE->Y() );
438 newNodeZ = FIRST_NODE->Z() + u * ( LAST_NODE->Z() - FIRST_NODE->Z() );
439 node_b = _theMeshDS->AddNode(newNodeX, newNodeY, newNodeZ);
440 newEdge = _theMeshDS->AddEdge(node_a, node_b);
441 edgesOnEdge.push_back(newEdge);
442 nodesOnEdge.push_back(node_b);
443 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
444 _nodeXx[ node_b ] = u;
445 if(MYDEBUG) MESSAGE("_nodeXx <-"<<u);
448 newEdge = _theMeshDS->AddEdge(node_a, LAST_NODE);
449 nodesOnEdge.push_back(LAST_NODE);
450 edgesOnEdge.push_back(newEdge);
452 _nodesOnEdge[&edge] = nodesOnEdge;
453 _edgesOnEdge[&edge] = edgesOnEdge;
455 if(MYDEBUG) MESSAGE("computeEdgeBySegment() : end >>>>>>>>");
460 // --------------------------------------------------------------
462 // --------------------------------------------------------------
463 std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document* doc )
465 std::map<HEXA_NS::Quad*, bool> quadWays;
466 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
467 std::list<HEXA_NS::Quad*> skinQuad;
468 std::list<HEXA_NS::Quad*> workingQuad;
469 HEXA_NS::Quad* first_q = NULL;
470 HEXA_NS::Quad* q = NULL;
471 HEXA_NS::Edge* e = NULL;
472 HEXA_NS::Vertex *e_0, *e_1 = NULL;
474 // FIRST STEP: eliminate free quad + internal quad
475 int nTotalQuad = doc->countUsedQuad();
476 for (int i=0; i < nTotalQuad; ++i ){
477 q = doc->getUsedQuad(i);
478 switch ( q->getNbrParents() ){ // parent == hexaedron
479 case 0: case 2: quadWays[q] = true; break;
480 case 1: skinQuad.push_back(q); break;
481 default: if ( q->getNbrParents() > 2 ) ASSERT(false);
485 // SECOND STEP: setting edges ways
486 while ( skinQuad.size()>0 ){
487 if(MYDEBUG) MESSAGE("SEARCHING INITIAL QUAD ..." );
488 for ( std::list<HEXA_NS::Quad*>::iterator it = skinQuad.begin(); it != skinQuad.end(); it++ ){
489 _searchInitialQuadWay( *it, e_0, e_1 );
490 if ( e_0 != NULL and e_1 != NULL ){
495 if ( e_0 == NULL and e_1 == NULL ) ASSERT(false);// should never happened,
496 if(MYDEBUG) MESSAGE("INITIAL QUAD FOUND!" );
497 for ( int j=0 ; j < 4 ; ++j ){
499 if ( ((e_0 == e->getVertex(0)) && (e_1 == e->getVertex(1)))
500 || ((e_0 == e->getVertex(1)) && (e_1 == e->getVertex(0))) ){
504 if(MYDEBUG) MESSAGE("INITIAL EDGE WAY FOUND!" );
506 edgeWays[e] = std::make_pair( e_0, e_1 );
507 workingQuad.push_back(q);
509 while ( workingQuad.size() > 0 ){
510 if(MYDEBUG) MESSAGE("COMPUTE QUAD WAY ... ID ="<< q->getId());
511 HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
513 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
514 while ( localEdgeWays.size() != 4 ){
515 HEXA_NS::Edge* e = q->getEdge(i%4);
516 if ( lastVertex == NULL ){
517 if ( edgeWays.count(e) == 1 ){
519 localEdgeWays[e] = std::make_pair( edgeWays[e].first, edgeWays[e].second );
521 localEdgeWays[e] = std::make_pair( edgeWays[e].second, edgeWays[e].first);
523 firstVertex = localEdgeWays[e].first;
524 lastVertex = localEdgeWays[e].second;
527 HEXA_NS::Vertex* e_0 = e->getVertex(0);
528 HEXA_NS::Vertex* e_1 = e->getVertex(1);
530 if ( lastVertex == e_0 ){
531 firstVertex = e_0; lastVertex = e_1;
532 } else if ( lastVertex == e_1 ){
533 firstVertex = e_1; lastVertex = e_0;
534 } else if ( firstVertex == e_0 ) {
535 firstVertex = e_1; lastVertex = e_0;
536 } else if ( firstVertex == e_1 ) {
537 firstVertex = e_0; lastVertex = e_1;
541 localEdgeWays[e] = std::make_pair( firstVertex, lastVertex );
542 if ( edgeWays.count(e) == 0 ){ // keep current value if present otherwise add it
543 edgeWays[e] = localEdgeWays[e];
550 //add new working quad
551 for ( int i=0 ; i < 4; ++i ){
552 HEXA_NS::Quad* next_q = NULL;
553 HEXA_NS::Edge* e = q->getEdge(i);
554 for ( int j=0 ; j < e->getNbrParents(); ++j ){
555 next_q = e->getParent(j);
559 int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
561 int fromWorkingQuad = std::count( workingQuad.begin(), workingQuad.end(), next_q );
562 if ( fromWorkingQuad == 0 ){
563 workingQuad.push_front( next_q );
570 HEXA_NS::Edge* e0 = q->getEdge(0);
571 HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
573 if ( e0_0 == localEdgeWays[ e0 ].first ){
575 } else if ( e0_0 == localEdgeWays[ e0 ].second ){
580 workingQuad.remove( q );
581 skinQuad.remove( q );
582 q = workingQuad.back();
593 // std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document& doc, std::map<HEXA_NS::Quad*, bool> initQuads )
595 // std::map<HEXA_NS::Quad*, bool> quadWays;
596 // // std::map<HEXA_NS::Edge*, bool> edgeWays;
597 // // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
598 // std::map<HEXA_NS::Quad*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > workingQuads;
600 // std::list<HEXA_NS::Quad*> skinQuad;
601 // std::list<HEXA_NS::Quad*> notSkinQuad;
602 // // std::list<HEXA_NS::Quad*> workingQuad;
603 // HEXA_NS::Quad* first_q = NULL;
604 // HEXA_NS::Quad* q = NULL;
605 // HEXA_NS::Edge* e = NULL;
606 // HEXA_NS::Vertex *e_0, *e_1 = NULL;
608 // // FIRST STEP: eliminate free quad + internal quad
609 // int nTotalQuad = doc.countQuad();
610 // for (int i=0; i < nTotalQuad; ++i ){
611 // q = doc.getQuad(i);
612 // switch ( q->getNbrParents() ){ // parent == hexaedron
613 // case 0: case 2: quadWays[q] = true; break;
614 // // case 0: case 2: notSkinQuad.push_back(q); break; //CS_TEST
615 // case 1: skinQuad.push_back(q); break;
616 // default: if ( q->getNbrParents() > 2 ) ASSERT(false);
622 // q = first_q = skinQuad.front();
623 // e = q->getUsedEdge(0);
624 // e_0 = e->getVertex(0);
625 // e_1 = e->getVertex(1);
627 // workingQuads[q] = std::make_pair( e_0, e_1 );
629 // while ( workingQuads.size() > 0 ){
630 // MESSAGE("CURRENT QUAD ------>"<< q->getId());
631 // HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
634 // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
635 // while ( localEdgeWays.size() != 4 ){
636 // HEXA_NS::Edge* e = q->getUsedEdge(i%4);
637 // if ( lastVertex == NULL ){
638 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
639 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
641 // if ( (workingQuads[q].first == e_0 and workingQuads[q].second == e_1)
642 // or (workingQuads[q].first == e_1 and workingQuads[q].second == e_0) ){
643 // if ( q == first_q ){
644 // localEdgeWays[e] = std::make_pair( workingQuads[q].first, workingQuads[q].second );
645 // MESSAGE("FIRST QUAD ");
647 // localEdgeWays[e] = std::make_pair( workingQuads[q].second, workingQuads[q].first);
648 // MESSAGE("NOT FIRST QUAD ");
650 // firstVertex = localEdgeWays[e].first;
651 // lastVertex = localEdgeWays[e].second;
654 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
655 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
656 // if ( lastVertex == e_0 ){
657 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
658 // firstVertex = e_0;
660 // } else if ( lastVertex == e_1 ){
661 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
662 // firstVertex = e_1;
664 // } else if ( firstVertex == e_0 ) {
665 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
666 // firstVertex = e_1;
668 // } else if ( firstVertex == e_1 ) {
669 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
670 // firstVertex = e_0;
680 // //add new working quad
681 // for ( int i=0 ; i < 4; ++i ){
682 // HEXA_NS::Quad* next_q = NULL;
683 // HEXA_NS::Edge* e = q->getUsedEdge(i);
684 // MESSAGE("NB PARENTS ="<< e->getNbrParents() );
685 // for ( int j=0 ; j < e->getNbrParents(); ++j ){
686 // next_q = e->getParent(j);
687 // if ( next_q == q ){
690 // int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
691 // if (fromSkin != 0){
692 // // int fromWorkingQuad = std::count( workingQuads.begin(), workingQuads.end(), next_q );
693 // int fromWorkingQuad = workingQuads.count( next_q );
694 // // MESSAGE("CHECK QUAD:"<< newWorkingQuad->getId());
695 // if ( fromWorkingQuad == 0 ){
696 // // workingQuads.push_front( next_q );
697 // workingQuads[ next_q ] = localEdgeWays[e];
698 // // MESSAGE("EDGE :"<<e->getId()<<"ADD QUAD :"<< newWorkingQuad->getId());
704 // //setting quad way
705 // HEXA_NS::Edge* e0 = q->getUsedEdge(0);
706 // HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
708 // if ( e0_0 == localEdgeWays[ e0 ].first ){
709 // quadWays[q] = true;
710 // } else if ( e0_0 == localEdgeWays[ e0 ].second ){
711 // quadWays[q] = false;
715 // MESSAGE("quadWays ID ="<< q->getId() << ", WAY = " << quadWays[q] );
717 // // workingQuad.remove( q );
718 // workingQuads.erase( q );
719 // skinQuad.remove( q );
720 // *workingQuads.begin();
721 // q = (*workingQuads.begin()).first;
727 bool SMESH_HexaBlocks::computeQuad( HEXA_NS::Quad& quad, bool way )
731 ok = computeQuadByAssoc(quad, way);
733 ok = computeQuadByFindingGeom(quad, way);
736 ok = computeQuadByLinearApproximation(quad, way);
739 _computeQuadOK = true;
745 bool SMESH_HexaBlocks::computeQuadByAssoc( HEXA_NS::Quad& quad, bool way )
747 // int id = quad.getId();
748 // if ( id != 11 ) return false; //7
750 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
751 MESSAGE("quadID = "<<quad.getId());
753 ASSERT( _computeEdgeOK );
756 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
757 SMESHFaces facesOnQuad;
758 SMDS_MeshFace* newFace = NULL;
759 std::vector<double> xx, yy;
761 // Elements for quad computation
762 SMDS_MeshNode *S1, *S2, *S4, *S3;
764 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
765 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
766 if ( initOk == false ){
770 const std::vector <HEXA_NS::Shape*> shapes = quad.getAssociations();
771 if ( shapes.size() == 0 ){
772 if(MYDEBUG) MESSAGE("computeQuadByAssoc() : end >>>>>>>>");
775 TopoDS_Shape shapeOrCompound = _getShapeOrCompound( shapes );
776 // bool quadWay = _computeQuadWay( quad, S1, S2, S3, S4, &shapeOrCompound );
777 // bool quadWay = _computeQuadWay( quad );
780 std::map<SMDS_MeshNode*, gp_Pnt> interpolatedPoints;
781 int iSize = nodesOnQuad.size();
782 int jSize = nodesOnQuad[0].size();
784 S1 = nodesOnQuad[0][0];
785 // S2 = nodesOnQuad[bNodes.size()-1][0];
786 S2 = nodesOnQuad[iSize-1][0];
787 S4 = nodesOnQuad[0][jSize-1];
788 S3 = nodesOnQuad[iSize-1][jSize-1];
791 for (int j = 1; j < jSize; ++j){
792 for (int i = 1; i < iSize; ++i){
793 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
794 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
795 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
796 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
799 double newNodeX, newNodeY, newNodeZ;
800 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
801 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
802 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
803 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
807 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
808 gp_Pnt newPt = gp_Pnt( newNodeX, newNodeY, newNodeZ );//interpolated point
811 if ( interpolatedPoints.count(n1) > 0 ){
812 pt1 = interpolatedPoints[n1];
814 pt1 = gp_Pnt( n1->X(), n1->Y(), n1->Z() );
816 if ( interpolatedPoints.count(n3) > 0 ){
817 pt3 = interpolatedPoints[n3];
819 pt3 = gp_Pnt( n3->X(), n3->Y(), n3->Z() );
821 gp_Vec vec1( newPt, pt1 );
822 gp_Vec vec2( newPt, pt3 );
824 gp_Pnt ptOnShape = _intersect(newPt, vec1, vec2, shapeOrCompound);
825 newNodeX = ptOnShape.X();
826 newNodeY = ptOnShape.Y();
827 newNodeZ = ptOnShape.Z();
828 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
829 nodesOnQuad[i][j] = n4;
830 interpolatedPoints[ n4 ] = newPt;
833 MESSAGE("u parameter is "<<u);
834 MESSAGE("v parameter is "<<v);
835 MESSAGE("point interpolated ("<<newPt.X()<<","<<newPt.Y()<<","<<newPt.Z()<<" )");
836 MESSAGE("point on shape ("<<newNodeX<<","<<newNodeY<<","<<newNodeZ<<" )");
841 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
842 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
843 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
844 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
848 if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
849 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
851 if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
852 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
854 facesOnQuad.push_back(newFace);
857 _quadNodes[ &quad ] = nodesOnQuad;
858 _facesOnQuad[&quad] = facesOnQuad;
860 if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
865 bool SMESH_HexaBlocks::computeQuadByFindingGeom( HEXA_NS::Quad& quad, bool way )
867 if(MYDEBUG) MESSAGE("computeQuadByFindingGeom() not implemented");
871 bool SMESH_HexaBlocks::_computeQuadInit(
873 ArrayOfSMESHNodes& nodesOnQuad,
874 std::vector<double>& xx, std::vector<double>& yy)
876 if(MYDEBUG) MESSAGE("_computeQuadInit() : begin ---------------");
879 SMDS_MeshNode *S1, *S2, *S4, *S3;
880 HEXA_NS::Edge *eh, *eb, *eg, *ed;
881 HEXA_NS::Edge *e1, *e2, *e3, *e4;
882 HEXA_NS::Vertex *e1_0, *e1_1, *e2_0, *e2_1, *e3_0, *e3_1, *e4_0, *e4_1;
884 e1 = quad.getEdge(0);
885 e2 = quad.getEdge(1);
886 e3 = quad.getEdge(2);
887 e4 = quad.getEdge(3);
889 e1_0 = e1->getVertex(0); e1_1 = e1->getVertex(1);
890 e2_0 = e2->getVertex(0); e2_1 = e2->getVertex(1);
891 e3_0 = e3->getVertex(0); e3_1 = e3->getVertex(1);
892 e4_0 = e4->getVertex(0); e4_1 = e4->getVertex(1);
895 S1 = _node[e1_0]; S2 = _node[e1_1];
901 } else if ( e1_0 == e2_1 ){
904 } else if ( e1_0 == e4_0 ){
907 } else if ( e1_0 == e4_1 ){
914 if ( S4 == _node[e3_0] ){
916 } else if ( S4 == _node[e3_1] ){
922 SMESHNodes hNodes = _nodesOnEdge[eh];
923 SMESHNodes bNodes = _nodesOnEdge[eb];
924 SMESHNodes gNodes = _nodesOnEdge[eg];
925 SMESHNodes dNodes = _nodesOnEdge[ed];
926 nodesOnQuad.resize( bNodes.size(), SMESHNodes(gNodes.size(), static_cast<SMDS_MeshNode*>(NULL)) );
930 // int &b_i = i, &h_i = i, &g_j = j, &d_j = j;
931 int *b_i = &i, *h_i = &i, *g_j = &j, *d_j = &j;
932 bool uWay = true, vWay = true;
934 if ( bNodes[0] != S1 ){
937 ASSERT( bNodes[bNodes.size()-1] == S1 );
939 ASSERT( bNodes[0] == S1);
941 if ( hNodes[0] != S4 ){
944 ASSERT( hNodes[0] == S4 );
946 if ( gNodes[0] != S1 ){
950 ASSERT( gNodes[0] == S1 );
952 if ( dNodes[0] != S2 ){
955 ASSERT( dNodes[0] == S2 );
960 for (i = 0, _i = bNodes.size()-1; i < bNodes.size(); ++i, --_i){
961 nodesOnQuad[i][0] = bNodes[*b_i];
962 nodesOnQuad[i][gNodes.size()-1 ] = hNodes[*h_i];
964 u = _nodeXx[ bNodes[*b_i] ];
971 if ( S1 != nodesOnQuad[0][0] ){
972 if(MYDEBUG) MESSAGE("ZZZZZZZZZZZZZZZZ quadID = "<<quad.getId());
974 // ASSERT( S1 == nodesOnQuad[0][0] );
978 for (j = 0, _j = gNodes.size()-1; j < gNodes.size(); ++j, --_j){
979 nodesOnQuad[0][j] = gNodes[*g_j];
980 if ( S1 != nodesOnQuad[0][0] ){
981 if(MYDEBUG) MESSAGE("XXXXXXXXXXXXXXXX quadID = "<<quad.getId());
983 // ASSERT( S1 == nodesOnQuad[0][0] );
984 nodesOnQuad[bNodes.size()-1][j] = dNodes[*d_j];
985 v = _nodeXx[ gNodes[*g_j] ];
994 int iSize = nodesOnQuad.size();
995 int jSize = nodesOnQuad[0].size();
996 ASSERT( iSize = bNodes.size() );
997 ASSERT( jSize = gNodes.size() );
999 // ASSERT( S1 == nodesOnQuad[0][0] );
1000 // ASSERT( S2 == nodesOnQuad[iSize-1][0]);
1001 // ASSERT( S4 == nodesOnQuad[0][jSize-1]);
1002 // ASSERT( S3 == nodesOnQuad[iSize-1][jSize-1]);
1008 bool SMESH_HexaBlocks::computeQuadByLinearApproximation( HEXA_NS::Quad& quad, bool way )
1010 // int id = quad.getId();
1011 // if ( quad.getId() != 66 ) return false; //7, 41
1013 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
1014 MESSAGE("quadID = "<<quad.getId());
1016 ASSERT( _computeEdgeOK );
1019 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
1020 SMESHFaces facesOnQuad;
1021 SMDS_MeshFace* newFace = NULL;
1022 std::vector<double> xx, yy;
1024 // Elements for quad computation
1025 SMDS_MeshNode *S1, *S2, *S4, *S3;
1027 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
1028 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
1029 if ( initOk == false ){
1033 int iSize = nodesOnQuad.size();
1034 int jSize = nodesOnQuad[0].size();
1036 S1 = nodesOnQuad[0][0];
1037 // S2 = nodesOnQuad[bNodes.size()-1][0];
1038 S2 = nodesOnQuad[iSize-1][0];
1039 S4 = nodesOnQuad[0][jSize-1];
1040 S3 = nodesOnQuad[iSize-1][jSize-1];
1042 for (int j = 1; j < jSize; ++j){
1043 for (int i = 1; i < iSize; ++i){
1044 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
1045 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
1046 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
1047 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
1050 double newNodeX, newNodeY, newNodeZ;
1051 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
1052 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
1053 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
1054 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
1058 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
1059 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
1060 nodesOnQuad[i][j] = n4;
1064 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
1065 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
1066 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
1067 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
1071 if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
1072 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
1074 if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
1075 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
1077 facesOnQuad.push_back(newFace);
1080 _quadNodes[ &quad ] = nodesOnQuad;
1081 _facesOnQuad[&quad] = facesOnQuad;
1083 if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
1088 // --------------------------------------------------------------
1090 // --------------------------------------------------------------
1091 bool SMESH_HexaBlocks::computeHexa( HEXA_NS::Document* doc )
1093 if(MYDEBUG) MESSAGE("computeHexa() : : begin <<<<<<");
1096 SMESH_MesherHelper aHelper(*_theMesh);
1097 TopoDS_Shape shape = _theMesh->GetShapeToMesh();
1098 aHelper.SetSubShape( shape );
1099 aHelper.SetElementsOnShape( true );
1101 SMESH_Gen* gen = _theMesh->GetGen();
1102 SMESH_HexaFromSkin_3D algo( gen->GetANewId(), 0, gen, doc );
1103 algo.InitComputeError();
1105 ok = algo.Compute( *_theMesh, &aHelper, _volumesOnHexa, _node );
1107 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D error!!! ");
1110 MESSAGE("SMESH_HexaFromSkin_3D.comment = "<<algo.GetComputeError()->myComment);
1111 MESSAGE("computeHexa() : end >>>>>>>>");
1118 // --------------------------------------------------------------
1119 // Document computing
1120 // --------------------------------------------------------------
1121 bool SMESH_HexaBlocks::computeDoc( HEXA_NS::Document* doc )
1123 if(MYDEBUG) MESSAGE("computeDoc() : : begin <<<<<<");
1126 // A) Vertex computation
1128 int nVertex = doc->countUsedVertex();
1129 HEXA_NS::Vertex* vertex = NULL;
1131 for (int j=0; j <nVertex; ++j ){ //Computing each vertex of the document
1132 vertex = doc->getUsedVertex(j);
1133 ok = computeVertex(*vertex);
1136 // B) Edges computation
1138 HEXA_NS::Propagation* propa = NULL;
1139 HEXA_NS::Law* law = NULL;
1140 HEXA_NS::Edges edges;
1142 nbPropa = doc->countPropagation();
1143 for (int j=0; j < nbPropa; ++j ){//Computing each edge's propagations of the document
1144 propa = doc->getPropagation(j);
1145 edges = propa->getEdges();
1146 law = propa->getLaw();
1149 law = doc->getLaw(0); // default law
1151 for( HEXA_NS::Edges::const_iterator iter = edges.begin();
1152 iter != edges.end();
1154 ok = computeEdge(**iter, *law);
1157 // C) Quad computation
1158 std::map<HEXA_NS::Quad*, bool> quadWays = computeQuadWays(doc);
1159 int nQuad = doc->countUsedQuad();
1160 HEXA_NS::Quad* q = NULL;
1161 for (int j=0; j <nQuad; ++j ){ //Computing each quad of the document
1162 q = doc->getUsedQuad(j);
1163 int id = q->getId();
1164 if ( quadWays.count(q) > 0 )
1165 ok = computeQuad( *q, quadWays[q] );
1167 if(MYDEBUG) MESSAGE("NO QUAD WAY ID = "<<id);
1171 // D) Hexa computation: Calling HexaFromSkin algo
1172 ok = computeHexa(doc);
1174 if(MYDEBUG) MESSAGE("computeDoc() : end >>>>>>>>");
1179 void SMESH_HexaBlocks::buildGroups(HEXA_NS::Document* doc)
1182 MESSAGE("_addGroups() : : begin <<<<<<");
1183 MESSAGE("_addGroups() : : nb. hexas= " << doc->countUsedHexa());
1184 MESSAGE("_addGroups() : : nb. quads= " << doc->countUsedQuad());
1185 MESSAGE("_addGroups() : : nb. edges= " << doc->countUsedEdge());
1186 MESSAGE("_addGroups() : : nb. nodes= " << doc->countUsedVertex());
1188 // Looping on each groups of the document
1189 for ( int i=0; i < doc->countGroup(); i++ ){
1190 _fillGroup( doc->getGroup(i) );
1193 if(MYDEBUG) MESSAGE("_addGroups() : end >>>>>>>>");
1196 // --------------------------------------------------------------
1198 // --------------------------------------------------------------
1199 double SMESH_HexaBlocks::_Xx( double i, HEXA_NS::Law law, double nbNodes) //, double pos0 )
1204 HEXA_NS::KindLaw k = law.getKind();
1205 double coeff = law.getCoefficient();
1207 case HEXA_NS::Uniform:
1208 result = (i+1)/(nbNodes+1);
1209 if(MYDEBUG) MESSAGE( "_Xx():" << " Uniform u("<<i<< ")"<< " = " << result);
1211 case HEXA_NS::Arithmetic:
1212 u0 = 1./(nbNodes + 1.) - (coeff*nbNodes)/2.;
1214 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Arithmetic discretization : check coefficient")));
1218 result = (i + 1.)*u0 + coeff*i*(i+1.)/2.;
1220 if(MYDEBUG) MESSAGE( "_Xx():" << " Arithmetic u("<<i<< ")"<< " = " << result);
1222 case HEXA_NS::Geometric:
1223 u0 = (1.-coeff)/(1.-pow(coeff, nbNodes + 1) ) ;
1225 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Geometric discretization : check coefficient")));
1229 result = u0 * (1.- pow(coeff, i + 1) )/(1.-coeff) ;
1231 if(MYDEBUG) MESSAGE( "_Xx():" << " Geometric u("<<i<< ")"<< " = " << result);
1238 double SMESH_HexaBlocks::_edgeLength(const TopoDS_Edge & E)
1240 if(MYDEBUG) MESSAGE("_edgeLength() : : begin <<<<<<");
1241 double UMin = 0, UMax = 0;
1242 if (BRep_Tool::Degenerated(E))
1245 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
1246 GeomAdaptor_Curve AdaptCurve(C);
1247 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
1248 if(MYDEBUG) MESSAGE("_edgeLength() : end >>>>>>>>");
1253 void SMESH_HexaBlocks::_buildMyCurve(
1254 const std::vector <HEXA_NS::Shape*>& associations, //IN
1255 const gp_Pnt& myCurve_start, //IN
1256 const gp_Pnt& myCurve_end, //IN
1257 std::list< BRepAdaptor_Curve* >& myCurve, //INOUT
1258 double& myCurve_length, //INOUT
1259 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//INOUT
1260 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //INOUT
1261 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //INOUT
1262 HEXA_NS::Edge& edge) // For error diagnostic
1264 if(MYDEBUG) MESSAGE("_buildMyCurve() : : begin <<<<<<");
1265 bool myCurve_way = true;
1266 myCurve_length = 0.;
1267 BRepAdaptor_Curve* thePreviousCurve = NULL;
1268 BRepAdaptor_Curve* theCurve = NULL;
1270 gp_Pnt theCurve_start, theCurve_end;
1271 gp_Pnt thePreviousCurve_start , thePreviousCurve_end;
1273 for ( std::vector <HEXA_NS::Shape*> ::const_iterator assoc = associations.begin();
1274 assoc != associations.end();
1276 string theBrep = (*assoc)->getBrep();
1277 double ass_debut = std::min ((*assoc)->debut, (*assoc)->fin);
1278 double ass_fin = std::max ((*assoc)->debut, (*assoc)->fin);
1279 TopoDS_Shape theShape = string2shape( theBrep );
1280 TopoDS_Edge theEdge = TopoDS::Edge( theShape );
1281 double theCurve_length = _edgeLength( theEdge );
1283 MESSAGE("_edgeLength ->"<<theCurve_length);
1285 if ( theCurve_length > 0 ){
1287 Handle(Geom_Curve) testCurve = BRep_Tool::Curve(theEdge, f, l);
1288 theCurve = new BRepAdaptor_Curve( theEdge );
1290 GCPnts_AbscissaPoint discret_start (*theCurve,
1291 theCurve_length*ass_debut,
1292 theCurve->FirstParameter() );
1293 GCPnts_AbscissaPoint discret_end (*theCurve,
1294 theCurve_length*ass_fin,
1295 theCurve->FirstParameter() );
1296 double u_start = discret_start.Parameter();
1297 double u_end = discret_end.Parameter();
1298 ASSERT( discret_start.IsDone() && discret_end.IsDone() );
1299 theCurve_start = theCurve->Value( u_start);
1300 theCurve_end = theCurve->Value( u_end );
1301 theCurve_length = theCurve_length*( ass_fin - ass_debut);
1304 MESSAGE("testCurve->f ->"<<f);
1305 MESSAGE("testCurve->l ->"<<l);
1306 MESSAGE("testCurve->FirstParameter ->"<<testCurve->FirstParameter());
1307 MESSAGE("testCurve->LastParameter ->"<<testCurve->LastParameter());
1309 MESSAGE("FirstParameter ->"<<theCurve->FirstParameter());
1310 MESSAGE("LastParameter ->"<<theCurve->LastParameter());
1311 MESSAGE("theCurve_length ->"<<theCurve_length);
1312 MESSAGE("(*assoc)->debut ->"<< ass_debut );
1313 MESSAGE("(*assoc)->fin ->"<< ass_fin );
1314 MESSAGE("u_start ->"<<u_start);
1315 MESSAGE("u_end ->"<<u_end);
1316 MESSAGE("myCurve_start( "<<myCurve_start.X()<<", "<<myCurve_start.Y()<<", "<<myCurve_start.Z()<<") ");
1317 MESSAGE("theCurve_start( "<<theCurve_start.X()<<", "<<theCurve_start.Y()<<", "<<theCurve_start.Z()<<") ");
1318 MESSAGE("myCurve_end( "<<myCurve_end.X()<<", "<<myCurve_end.Y()<<", "<<myCurve_end.Z()<<") ");
1319 MESSAGE("theCurve_end( "<<theCurve_end.X()<<", "<<theCurve_end.Y()<<", "<<theCurve_end.Z()<<") ");
1322 if ( thePreviousCurve == NULL ){
1323 // setting myCurve_way and first curve way
1324 if ( myCurve_start.IsEqual(theCurve_start, HEXA_EPSILON) ){
1325 if(MYDEBUG) MESSAGE("myCurve_start.IsEqual(theCurve_start, HEXA_EPSILON)");
1327 myCurve_ways[theCurve] = true;
1328 } else if ( myCurve_start.IsEqual(theCurve_end, HEXA_EPSILON) ){
1329 if(MYDEBUG) MESSAGE("myCurve_start.IsEqual(theCurve_end, HEXA_EPSILON)");
1331 myCurve_ways[theCurve] = false;
1332 } else if ( myCurve_end.IsEqual(theCurve_end, HEXA_EPSILON) ){
1333 if(MYDEBUG) MESSAGE("myCurve_end.IsEqual(theCurve_end, HEXA_EPSILON)");
1334 myCurve_way = false;
1335 myCurve_ways[theCurve] = true;
1336 } else if ( myCurve_end.IsEqual(theCurve_start, HEXA_EPSILON) ){
1337 if(MYDEBUG) MESSAGE("myCurve_end.IsEqual(theCurve_start, HEXA_EPSILON)");
1338 myCurve_way = false;
1339 myCurve_ways[theCurve] = false;
1341 if(MYDEBUG) MESSAGE("SOMETHING WRONG on edge association... Bad script?");
1344 throw (SALOME_Exception(LOCALIZED("Edge association : check association parameters ( first, last ) between HEXA model and CAO")));
1348 // it is not the first or last curve.
1349 // ways are calculated between previous and new one.
1350 if ( thePreviousCurve_end.IsEqual( theCurve_end, HEXA_EPSILON )
1351 or thePreviousCurve_start.IsEqual( theCurve_start, HEXA_EPSILON ) ){
1352 myCurve_ways[theCurve] = !myCurve_ways[thePreviousCurve];// opposite WAY
1353 if(MYDEBUG) MESSAGE("opposite WAY");
1354 } else if ( thePreviousCurve_end.IsEqual( theCurve_start, HEXA_EPSILON )
1355 or thePreviousCurve_start.IsEqual( theCurve_end, HEXA_EPSILON ) ){
1356 myCurve_ways[theCurve] = myCurve_ways[thePreviousCurve];// same WAY
1357 if(MYDEBUG) MESSAGE("same WAY");
1359 if(MYDEBUG) MESSAGE("SOMETHING WRONG on edge association... bad script?");
1361 throw (SALOME_Exception(LOCALIZED("Edge association : Check association parameters ( first, last ) between HEXA model and CAO")));
1365 myCurve_starts[theCurve] = u_start;
1366 myCurve_lengths[theCurve] = theCurve_length;
1367 myCurve_length += theCurve_length;
1368 myCurve.push_back( theCurve );
1370 thePreviousCurve = theCurve;
1371 thePreviousCurve_start = theCurve_start;
1372 thePreviousCurve_end = theCurve_end;
1374 }//if ( theCurveLength > 0 ){
1379 if ( myCurve_way == false ){
1380 std::list< BRepAdaptor_Curve* > tmp( myCurve.size() );
1381 std::copy( myCurve.rbegin(), myCurve.rend(), tmp.begin() );
1386 MESSAGE("myCurve_way was :"<<myCurve_way);
1387 MESSAGE("_buildMyCurve() : end >>>>>>>>");
1394 gp_Pnt SMESH_HexaBlocks::_getPtOnMyCurve(
1395 const double& myCurve_u, //IN
1396 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //IN
1397 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//IN
1398 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //IN
1399 std::list< BRepAdaptor_Curve* >& myCurve, //INOUT
1400 double& myCurve_start ) //INOUT
1401 // std::map< BRepAdaptor_Curve*, double>& myCurve_firsts,
1402 // std::map< BRepAdaptor_Curve*, double>& myCurve_lasts,
1404 if(MYDEBUG) MESSAGE("_getPtOnMyCurve() : : begin <<<<<<");
1407 // looking for curve which contains parameter myCurve_u
1408 BRepAdaptor_Curve* curve = myCurve.front();
1409 double curve_start = myCurve_start;
1410 double curve_end = curve_start + myCurve_lengths[curve];
1412 GCPnts_AbscissaPoint discret;
1415 MESSAGE("looking for curve: c = "<<myCurve_u);
1416 MESSAGE("looking for curve: curve_start = "<<curve_start);
1417 MESSAGE("looking for curve: curve_end = "<<curve_end);
1418 MESSAGE("looking for curve: curve_lenght = "<<myCurve_lengths[curve]);
1419 MESSAGE("looking for curve: curve.size _lenght= "<<myCurve.size());
1421 while ( not ( (myCurve_u >= curve_start) and (myCurve_u <= curve_end) ) ) {
1422 // if (myCurve.size() == 0 )
1424 // PutData (myCurve.size());
1425 // PutData (myCurve_u);
1428 ASSERT( myCurve.size() != 0 );
1429 myCurve.pop_front();
1430 curve = myCurve.front();
1431 curve_start = curve_end;
1432 curve_end = curve_start + myCurve_lengths[curve];
1434 MESSAGE("go next curve: curve_lenght = "<<myCurve_lengths[curve]);
1435 MESSAGE("go next curve: curve_start = "<<curve_start);
1436 MESSAGE("go next curve: curve_end = "<<curve_end);
1437 MESSAGE("go next curve: myCurve_u = "<<myCurve_u);
1440 myCurve_start = curve_start;
1443 if ( myCurve_ways[curve] ){
1444 // curve_u = myCurve_firsts[curve] + (myCurve_u - curve_start);
1445 // discret = GCPnts_AbscissaPoint( *curve, (myCurve_u - curve_start), curve->FirstParameter() );
1446 discret = GCPnts_AbscissaPoint( *curve, (myCurve_u - curve_start), myCurve_starts[curve] );
1448 // discret = GCPnts_AbscissaPoint( *curve, myCurve_lengths[curve]- (myCurve_u - curve_start), curve->FirstParameter() );
1449 discret = GCPnts_AbscissaPoint( *curve, myCurve_lengths[curve]- (myCurve_u - curve_start), myCurve_starts[curve] );
1451 ASSERT(discret.IsDone());
1452 curve_u = discret.Parameter();
1453 ptOnMyCurve = curve->Value( curve_u );
1456 MESSAGE("curve found!");
1457 MESSAGE("curve_u = "<< curve_u);
1458 MESSAGE("curve way = "<< myCurve_ways[curve]);
1459 MESSAGE("_getPtOnMyCurve() : end >>>>>>>>");
1470 void SMESH_HexaBlocks::_nodeInterpolationUV(double u, double v,
1471 SMDS_MeshNode* Pg, SMDS_MeshNode* Pd, SMDS_MeshNode* Ph, SMDS_MeshNode* Pb,
1472 SMDS_MeshNode* S1, SMDS_MeshNode* S2, SMDS_MeshNode* S3, SMDS_MeshNode* S4,
1473 double& xOut, double& yOut, double& zOut )
1476 MESSAGE("_nodeInterpolationUV() IN:");
1477 MESSAGE("u ( "<< u <<" )");
1478 MESSAGE("v ( "<< v <<" )");
1480 MESSAGE("S1 (" << S1->X() << "," << S1->Y() << "," << S1->Z() << ")");
1481 MESSAGE("S2 (" << S2->X() << "," << S2->Y() << "," << S2->Z() << ")");
1482 MESSAGE("S4 (" << S4->X() << "," << S4->Y() << "," << S4->Z() << ")");
1483 MESSAGE("S3 (" << S3->X() << "," << S3->Y() << "," << S3->Z() << ")");
1485 MESSAGE("Pg (" << Pg->X() << "," << Pg->Y() << "," << Pg->Z() << ")");
1486 MESSAGE("Pd (" << Pd->X() << "," << Pd->Y() << "," << Pd->Z() << ")");
1487 MESSAGE("Ph (" << Ph->X() << "," << Ph->Y() << "," << Ph->Z() << ")");
1488 MESSAGE("Pb (" << Pb->X() << "," << Pb->Y() << "," << Pb->Z() << ")");
1491 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();
1492 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();
1493 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();
1496 MESSAGE("_nodeInterpolationUV() OUT("<<xOut<<","<<yOut<<","<<zOut<<" )");
1501 TopoDS_Shape SMESH_HexaBlocks::_getShapeOrCompound( const std::vector<HEXA_NS::Shape*>& shapesIn)
1503 ASSERT( shapesIn.size()!=0 );
1505 if (shapesIn.size() == 1) {
1506 HEXA_NS::Shape* assoc = shapesIn.front();
1507 string strBrep = assoc->getBrep();
1508 return string2shape( strBrep );
1510 TopoDS_Compound aCompound;
1511 BRep_Builder aBuilder;
1512 aBuilder.MakeCompound( aCompound );
1514 for ( std::vector <HEXA_NS::Shape*> ::const_iterator assoc = shapesIn.begin();
1515 assoc != shapesIn.end();
1517 string strBrep = (*assoc)->getBrep();
1518 TopoDS_Shape shape = string2shape( strBrep );
1519 aBuilder.Add( aCompound, shape );
1526 // ================================================== carre
1527 inline double carre (double val)
1531 // ================================================== dist2
1532 inline double dist2 (const gp_Pnt& pt1, const gp_Pnt& pt2)
1534 double dist = carre (pt2.X()-pt1.X()) + carre (pt2.Y()-pt1.Y())
1535 + carre (pt2.Z()-pt1.Z());
1538 // ================================================== _intersect
1539 gp_Pnt SMESH_HexaBlocks::_intersect( const gp_Pnt& Pt,
1540 const gp_Vec& u, const gp_Vec& v,
1541 const TopoDS_Shape& shape,
1546 gp_Vec normale = u^v;
1547 gp_Dir dir(normale);
1548 gp_Lin li( Pt, dir );
1550 Standard_Real s = -Precision::Infinite();
1551 Standard_Real e = +Precision::Infinite();
1553 IntCurvesFace_ShapeIntersector inter;
1554 inter.Load(shape, tol);
1555 // inter.Load(S, tol);
1556 inter.Perform(li, s, e);//inter.PerformNearest(li, s, e);
1558 /************************************************************** Abu 2011-11-04 */
1559 /// if ( inter.IsDone() && (inter.NbPnt()==1) ) {
1560 if ( inter.IsDone() )
1562 result = inter.Pnt(1);//first
1563 int nbrpts = inter.NbPnt();
1566 double d0 = dist2 (result, Pt);
1567 for (int i=2; i <= inter.NbPnt(); ++i )
1569 double d1 = dist2 (Pt, inter.Pnt(i));
1573 result = inter.Pnt (i);
1577 /************************************************************** Abu 2011-11-04 (fin) */
1579 MESSAGE("_intersect() : OK");
1580 for ( int i=1; i <= inter.NbPnt(); ++i ){
1581 gp_Pnt tmp = inter.Pnt(i);
1582 MESSAGE("_intersect() : pnt("<<i<<") = ("<<tmp.X()<<","<<tmp.Y()<<","<<tmp.Z()<<" )");
1587 if(MYDEBUG) MESSAGE("_intersect() : KO");
1596 // parameters q : IN, v0: INOUT, v1: INOUT
1597 void SMESH_HexaBlocks::_searchInitialQuadWay( HEXA_NS::Quad* q, HEXA_NS::Vertex*& v0, HEXA_NS::Vertex*& v1 )
1599 if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : begin");
1600 v0 = NULL; v1 = NULL;
1601 if ( q->getNbrParents() != 1 ) return; // q must be a skin quad
1603 HEXA_NS::Vertex* qA = q->getVertex(0);
1604 HEXA_NS::Vertex* qB = q->getVertex(1);
1605 HEXA_NS::Vertex* qC = q->getVertex(2);
1606 HEXA_NS::Vertex* qD = q->getVertex(3);
1608 // searching for vertex on opposed quad
1609 HEXA_NS::Vertex *qAA = NULL, *qBB = NULL, *qCC = NULL, *qDD = NULL;
1610 HEXA_NS::Hexa* h = q->getParent(0);
1611 for( int i=0; i < h->countEdge(); ++i ){
1612 HEXA_NS::Edge* e = h->getEdge(i);
1613 HEXA_NS::Vertex* e0 = e->getVertex(0);
1614 HEXA_NS::Vertex* e1 = e->getVertex(1);
1616 if ( e0 == qA and e1 != qB and e1 != qC and e1 != qD ){
1618 } else if ( e1 == qA and e0 != qB and e0 != qC and e0 != qD ){
1620 } else if ( e0 == qB and e1 != qA and e1 != qC and e1 != qD ){
1622 } else if ( e1 == qB and e0 != qA and e0 != qC and e0 != qD ){
1624 } else if ( e0 == qC and e1 != qA and e1 != qB and e1 != qD ){
1626 } else if ( e1 == qC and e0 != qA and e0 != qB and e0 != qD ){
1628 } else if ( e0 == qD and e1 != qA and e1 != qB and e1 != qC ){
1630 } else if ( e1 == qD and e0 != qA and e0 != qB and e0 != qC ){
1635 // working on final value ( point on CAO ), not on model
1636 SMDS_MeshNode *nA = _node[qA], *nAA = _node[qAA];
1637 SMDS_MeshNode *nB = _node[qB], *nBB = _node[qBB];
1638 SMDS_MeshNode *nC = _node[qC], *nCC = _node[qCC];
1639 SMDS_MeshNode *nD = _node[qD], *nDD = _node[qDD];
1641 gp_Pnt pA( nA->X(), nA->Y(), nA->Z() );
1642 gp_Pnt pB( nB->X(), nB->Y(), nB->Z() );
1643 gp_Pnt pC( nC->X(), nC->Y(), nC->Z() );
1644 gp_Pnt pD( nD->X(), nD->Y(), nD->Z() );
1646 gp_Pnt pAA( nAA->X(), nAA->Y(), nAA->Z() );
1647 gp_Pnt pBB( nBB->X(), nBB->Y(), nBB->Z() );
1648 gp_Pnt pCC( nCC->X(), nCC->Y(), nCC->Z() );
1649 gp_Pnt pDD( nDD->X(), nDD->Y(), nDD->Z() );
1651 gp_Vec AB( pA, pB );
1652 gp_Vec AC( pA, pC );
1653 gp_Vec normP = AB^AC;
1654 gp_Dir dirP( normP );
1656 // building plane for point projection
1657 gp_Pln plnP( gp_Pnt(nA->X(), nA->Y(), nA->Z()), dirP);
1658 TopoDS_Shape sPlnP = BRepBuilderAPI_MakeFace(plnP).Face();
1660 // PAAA is the result of PAA projection
1661 gp_Pnt pAAA = _intersect( pAA, AB, AC, sPlnP );
1662 gp_Pnt pBBB = _intersect( pBB, AB, AC, sPlnP );
1663 gp_Pnt pCCC = _intersect( pCC, AB, AC, sPlnP );
1664 gp_Pnt pDDD = _intersect( pDD, AB, AC, sPlnP );
1666 gp_Dir AA( gp_Vec(pAA, pAAA) );
1667 gp_Dir BB( gp_Vec(pBB, pBBB) );
1668 gp_Dir CC( gp_Vec(pCC, pCCC) );
1669 gp_Dir DD( gp_Vec(pDD, pDDD) );
1671 // eventually, we are able to know if the input quad is a good client!
1672 // exit the fonction otherwise
1673 if ( AA.IsOpposite(BB, HEXA_QUAD_WAY) ) return;
1674 if ( BB.IsOpposite(CC, HEXA_QUAD_WAY) ) return;
1675 if ( CC.IsOpposite(DD, HEXA_QUAD_WAY) ) return;
1677 // ok, give the input quad the good orientation by
1679 if ( !dirP.IsOpposite(AA, HEXA_QUAD_WAY) ) { //OK
1685 if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : end");
1688 SMESH_Group* SMESH_HexaBlocks::_createGroup(HEXA_NS::Group& grHex)
1690 if(MYDEBUG) MESSAGE("_createGroup() : : begin <<<<<<");
1692 std::string aGrName = grHex.getName();
1693 HEXA_NS::EnumGroup grHexKind = grHex.getKind();
1695 if(MYDEBUG) MESSAGE("aGrName"<<aGrName);
1697 SMDSAbs_ElementType aGrType;
1698 switch ( grHexKind ){
1699 case HEXA_NS::HexaCell : aGrType = SMDSAbs_Volume; break;
1700 case HEXA_NS::QuadCell : aGrType = SMDSAbs_Face ; break;
1701 case HEXA_NS::EdgeCell : aGrType = SMDSAbs_Edge ; break;
1702 case HEXA_NS::HexaNode : aGrType = SMDSAbs_Node ; break;
1703 case HEXA_NS::QuadNode : aGrType = SMDSAbs_Node ; break;
1704 case HEXA_NS::EdgeNode : aGrType = SMDSAbs_Node ; break;
1705 case HEXA_NS::VertexNode : aGrType = SMDSAbs_Node ; break;
1706 default : ASSERT(false);
1710 SMESH_Group* aGr = _theMesh->AddGroup(aGrType, aGrName.c_str(), aId);
1712 if(MYDEBUG) MESSAGE("_createGroup() : end >>>>>>>>");
1716 void SMESH_HexaBlocks::_fillGroup(HEXA_NS::Group* grHex )
1718 if(MYDEBUG) MESSAGE("_fillGroup() : : begin <<<<<<");
1720 SMESH_Group* aGr = _createGroup( *grHex );
1721 HEXA_NS::EltBase* grHexElt = NULL;
1722 HEXA_NS::EnumGroup grHexKind = grHex->getKind();
1723 int grHexNbElt = grHex->countElement();
1725 if(MYDEBUG) MESSAGE("_fillGroup() : kind = " << grHexKind);
1726 if(MYDEBUG) MESSAGE("_fillGroup() : count= " << grHexNbElt);
1728 // A)Looking for elements ID
1729 std::vector<const SMDS_MeshElement*> aGrEltIDs;
1731 for ( int n=0; n<grHexNbElt; ++n ){
1732 grHexElt = grHex->getElement(n);
1734 switch ( grHexKind ){
1735 case HEXA_NS::HexaCell:
1737 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1738 // HEXA_NS::Hexa* h = dynamic_cast<HEXA_NS::Hexa*>(grHexElt);
1740 if ( _volumesOnHexa.count(h)>0 ){
1741 SMESHVolumes volumes = _volumesOnHexa[h];
1742 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1743 aGrEltIDs.push_back(*aVolume);
1746 if(MYDEBUG) MESSAGE("GROUP OF VOLUME: volume for hexa (id = "<<h->getId()<<") not found");
1750 case HEXA_NS::QuadCell:
1752 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1753 // HEXA_NS::Quad* q = dynamic_cast<HEXA_NS::Quad*>(grHexElt);
1755 if ( _facesOnQuad.count(q)>0 ){
1756 SMESHFaces faces = _facesOnQuad[q];
1757 for ( SMESHFaces::iterator aFace = faces.begin(); aFace != faces.end(); ++aFace ){
1758 aGrEltIDs.push_back(*aFace);
1761 if(MYDEBUG) MESSAGE("GROUP OF FACE: face for quad (id = "<<q->getId()<<") not found");
1765 case HEXA_NS::EdgeCell:
1767 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1768 // HEXA_NS::Edge* e = dynamic_cast<HEXA_NS::Edge*>(grHexElt);
1770 if ( _edgesOnEdge.count(e)>0 ){
1771 SMESHEdges edges = _edgesOnEdge[e];
1772 for ( SMESHEdges::iterator anEdge = edges.begin(); anEdge != edges.end(); ++anEdge ){
1773 aGrEltIDs.push_back(*anEdge);
1776 if(MYDEBUG) MESSAGE("GROUP OF Edge: edge for edge (id = "<<e->getId()<<") not found");
1780 case HEXA_NS::HexaNode:
1782 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1783 // HEXA_NS::Hexa* h = dynamic_cast<HEXA_NS::Hexa*>(grHexElt);
1785 if ( _volumesOnHexa.count(h)>0 ){
1786 SMESHVolumes volumes = _volumesOnHexa[h];
1787 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1788 SMDS_ElemIteratorPtr aNodeIter = (*aVolume)->nodesIterator();
1789 while( aNodeIter->more() ){
1790 const SMDS_MeshNode* aNode =
1791 dynamic_cast<const SMDS_MeshNode*>( aNodeIter->next() );
1793 aGrEltIDs.push_back(aNode);
1798 if(MYDEBUG) MESSAGE("GROUP OF HEXA NODES: nodes on hexa (id = "<<h->getId()<<") not found");
1802 case HEXA_NS::QuadNode:
1804 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1805 // HEXA_NS::Quad* q = dynamic_cast<HEXA_NS::Quad*>(grHexElt);
1807 if ( _quadNodes.count(q)>0 ){
1808 ArrayOfSMESHNodes nodesOnQuad = _quadNodes[q];
1809 for ( ArrayOfSMESHNodes::iterator nodes = nodesOnQuad.begin(); nodes != nodesOnQuad.end(); ++nodes){
1810 for ( SMESHNodes::iterator aNode = nodes->begin(); aNode != nodes->end(); ++aNode){
1811 aGrEltIDs.push_back(*aNode);
1815 if(MYDEBUG) MESSAGE("GROUP OF QUAD NODES: nodes on quad (id = "<<q->getId()<<") not found");
1819 case HEXA_NS::EdgeNode:
1821 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1822 // HEXA_NS::Edge* e = dynamic_cast<HEXA_NS::Edge*>(grHexElt);
1824 if ( _nodesOnEdge.count(e)>0 ){
1825 SMESHNodes nodes = _nodesOnEdge[e];
1826 for ( SMESHNodes::iterator aNode = nodes.begin(); aNode != nodes.end(); ++aNode){
1827 aGrEltIDs.push_back(*aNode);
1830 if(MYDEBUG) MESSAGE("GROUP OF EDGE NODES: nodes on edge (id = "<<e->getId()<<") not found");
1834 case HEXA_NS::VertexNode:
1836 HEXA_NS::Vertex* v = reinterpret_cast<HEXA_NS::Vertex*>(grHexElt);
1837 // HEXA_NS::Vertex* v = dynamic_cast<HEXA_NS::Vertex*>(grHexElt);
1839 if ( _node.count(v)>0 ){
1840 aGrEltIDs.push_back(_node[v]);
1842 if(MYDEBUG) MESSAGE("GROUP OF VERTEX NODES: nodes for vertex (id = "<<v->getId()<<") not found");
1846 default : ASSERT(false);
1850 // B)Filling the group on SMESH
1851 SMESHDS_Group* aGroupDS = dynamic_cast<SMESHDS_Group*>( aGr->GetGroupDS() );
1853 for ( int i=0; i < aGrEltIDs.size(); i++ ) {
1854 aGroupDS->SMDSGroup().Add( aGrEltIDs[i] );
1857 if(MYDEBUG) MESSAGE("_fillGroup() : end >>>>>>>>");
1864 // not used, for backup purpose only:
1865 void SMESH_HexaBlocks::_getCurve( const std::vector<HEXA_NS::Shape*>& shapesIn,
1866 Handle_Geom_Curve& curveOut, double& curveFirstOut, double& curveLastOut )
1868 // std::cout<<"------------------- _getCurve ------------ "<<std::endl;
1869 GeomConvert_CompCurveToBSplineCurve* gen = NULL;
1871 double curvesLenght = 0.;
1872 double curvesFirst = shapesIn.front()->debut;
1873 double curvesLast = shapesIn.back()->fin;
1875 for ( std::vector <HEXA_NS::Shape*> ::const_iterator assoc = shapesIn.begin();
1876 assoc != shapesIn.end();
1878 string strBrep = (*assoc)->getBrep();
1879 TopoDS_Shape shape = string2shape( strBrep );
1880 TopoDS_Edge Edge = TopoDS::Edge(shape);
1882 Handle(Geom_Curve) curve = BRep_Tool::Curve(Edge, f, l);
1883 curvesLenght += l-f;
1884 Handle(Geom_BoundedCurve) bCurve = Handle(Geom_BoundedCurve)::DownCast(curve);
1886 gen = new GeomConvert_CompCurveToBSplineCurve(bCurve);
1888 bool bb=gen->Add(bCurve, Precision::Confusion(), Standard_True, Standard_False, 1);
1892 curveFirstOut = curvesFirst/curvesLenght;
1893 curveLastOut = curvesLenght - (1.-curvesLast)/curvesLenght;
1894 curveOut = gen->BSplineCurve();
1896 std::cout<<"curvesFirst -> "<<curvesFirst<<std::endl;
1897 std::cout<<"curvesLast -> "<<curvesLast<<std::endl;
1898 std::cout<<"curvesLenght -> "<<curvesLenght<<std::endl;
1899 std::cout<<"curveFirstOut -> "<<curveFirstOut<<std::endl;
1900 std::cout<<"curveLastOut -> "<<curveLastOut<<std::endl;
1908 // bool SMESH_HexaBlocks::_areSame(double a, double b)
1910 // return fabs(a - b) < HEXA_EPSILON;
1912 // // MESSAGE("Angular() :" << dir2.IsOpposite(dir1, Precision::Angular()));
1913 // // ASSERT( dir2.IsParallel(dir1, HEXA_QUAD_WAY) );
1914 // // bool test2 = norm2.IsOpposite(norm1, HEXA_QUAD_WAY2) == norm3.IsOpposite(norm1, HEXA_QUAD_WAY2);
1915 // // way = norm1.IsOpposite(norm3.Reversed(), HEXA_QUAD_WAY2);
1916 // gp_Pnt p( n->X(), n->Y(), n->Z() );
1917 // gp_Pnt ptOnPlane;
1918 // gp_Pnt ptOnSurface;
1919 // gp_Pnt ptOnPlaneOrSurface;
1920 // // gp_Vec norm2(p1, p);
1921 // TopoDS_Shape planeOrSurface;
1924 // gp_Pln pln(p1, dir1);
1925 // TopoDS_Shape shapePln = BRepBuilderAPI_MakeFace(pln).Face();
1926 // ptOnPlane = _intersect( p, a1, b1, shapePln );
1927 // ptOnPlaneOrSurface = ptOnPlane;
1930 // // if ( assoc != NULL ){
1931 // // MESSAGE("_computeQuadWay with assoc");
1932 // for( int i=0; i < h->countEdge(); ++i ){
1933 // HEXA_NS::Edge* e = h->getUsedEdge(i);
1934 // if ( e->definedBy(v1,v2) ){
1935 // const std::vector <HEXA_NS::Shape*> assocs = e->getAssociations();
1936 // if ( assocs.size() != 0 ){
1937 // HEXA_NS::Shape* assoc = assocs[0]; //CS_TODO
1938 // string theBrep = assoc->getBrep();
1939 // TopoDS_Shape theShape = string2shape( theBrep );
1940 // ptOnSurface = _intersect( p, a1, b1, theShape );
1941 // if ( !ptOnSurface.IsEqual(p, HEXA_EPSILON) ){
1942 // ptOnPlaneOrSurface = ptOnSurface;