1 // Copyright (C) 2009-2013 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 "HexGroup.hxx"
65 #include "hexa_base.hxx"
66 #include "HexAssoEdge.hxx"
68 // HEXABLOCKPLUGIN includes
69 #include "HEXABLOCKPlugin_mesh.hxx"
70 #include "HEXABLOCKPlugin_FromSkin_3D.hxx"
73 #include "Basics_Utils.hxx"
74 #include "utilities.h"
85 #define EXCEPTION(TYPE, MSG) {\
86 std::ostringstream aStream;\
87 aStream<<__FILE__<<"["<<__LINE__<<"]::"<<MSG;\
88 throw TYPE(aStream.str());\
93 static int MYDEBUG = HEXA_NS::on_debug ();
95 static int MYDEBUG = 0;
99 static double HEXA_EPS = 1.0e-3; //1E-3;
100 static double HEXA_QUAD_WAY = M_PI/4.; //3.*PI/8.;
102 // ============================================================ string2shape
103 TopoDS_Shape string2shape( const string& brep )
106 istringstream streamBrep(brep);
107 BRep_Builder aBuilder;
108 BRepTools::Read(shape, streamBrep, aBuilder);
112 // ============================================================ shape2coord
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 );
127 // ============================================================= edge_length
128 double edge_length(const TopoDS_Edge & E)
130 double UMin = 0, UMax = 0;
131 if (BRep_Tool::Degenerated(E))
134 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
135 GeomAdaptor_Curve AdaptCurve(C);
136 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
139 // =============================================================== make_curve
140 BRepAdaptor_Curve* make_curve (gp_Pnt& pstart, gp_Pnt& pend,
141 double& length, double& u_start,
142 HEXA_NS::Edge& edge, int nro)
144 Hex::AssoEdge* asso = edge.getAssociation (nro);
145 BRepAdaptor_Curve* the_curve = asso->getCurve();
146 length = asso->length ();
150 const double* point1 = asso->getOrigin ();
151 const double* point2 = asso->getExtrem ();
152 pstart = gp_Pnt (point1[0], point1[1], point1[2]);
153 pend = gp_Pnt (point2[0], point2[1], point2[2]);
154 u_start = asso->getUstart();
157 // ============================================================== Constructeur
158 // SMESH_HexaBlocks::SMESH_HexaBlocks( SMESH_Mesh* theMesh ):
159 SMESH_HexaBlocks::SMESH_HexaBlocks(SMESH_Mesh& theMesh):
163 _computeVertexOK(false),
164 _computeEdgeOK(false),
165 _computeQuadOK(false),
166 _theMesh(&theMesh), //groups creation
167 _theMeshDS(theMesh.GetMeshDS()) //meshing
172 SMESH_HexaBlocks::~SMESH_HexaBlocks()
177 // --------------------------------------------------------------
179 // --------------------------------------------------------------
181 // --------------------------------------------------------------
183 // --------------------------------------------------------------
184 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
185 // ============================================================== computeVertex
186 bool SMESH_HexaBlocks::computeVertex(HEXA_NS::Vertex& vx)
189 vx.getAssoCoord (px, py, pz);
191 SMDS_MeshNode* new_node = _theMeshDS->AddNode (px, py, pz);
192 _vertex [new_node] = &vx;
193 _node [&vx] = new_node; //needed in computeEdge()
194 _computeVertexOK = true;
197 // --------------------------------------------------------------
199 // --------------------------------------------------------------
200 bool SMESH_HexaBlocks::computeEdge(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
204 ok = computeEdgeByAssoc( edge, law);
206 ok = computeEdgeByShortestWire( edge, law);
209 ok = computeEdgeByPlanWire( edge, law);
212 ok = computeEdgeByIsoWire( edge, law);
215 ok = computeEdgeBySegment( edge, law);
218 _computeEdgeOK = true;
225 bool SMESH_HexaBlocks::computeEdgeByAssoc( HEXA_NS::Edge& edge, HEXA_NS::Law& law )
227 if(MYDEBUG) MESSAGE("computeEdgeByAssoc(edgeID = "<<edge.getId()<<"): begin <<<<<<");
228 ASSERT( _computeVertexOK );
231 if (NOT edge.isAssociated())
235 HEXA_NS::Vertex* vx0 = NULL;
236 HEXA_NS::Vertex* vx1 = NULL;
238 // way of discretization
239 if (edge.getWay() == true){
240 vx0 = edge.getVertex(0);
241 vx1 = edge.getVertex(1);
243 vx0 = edge.getVertex(1);
244 vx1 = edge.getVertex(0);
247 SMDS_MeshNode* FIRST_NODE = _node[vx0];
248 SMDS_MeshNode* LAST_NODE = _node[vx1];
251 // A) Build myCurve_list
252 std::list< BRepAdaptor_Curve* > myCurve_list;
253 double myCurve_tot_len;
254 std::map< BRepAdaptor_Curve*, double> myCurve_lengths;
255 std::map< BRepAdaptor_Curve*, bool> myCurve_ways;
256 std::map< BRepAdaptor_Curve*, double> myCurve_starts;
258 gp_Pnt myCurve_pt_start(FIRST_NODE->X(), FIRST_NODE->Y(), FIRST_NODE->Z());
259 gp_Pnt myCurve_pt_end (LAST_NODE->X(), LAST_NODE->Y(), LAST_NODE->Z());
273 // B) Build nodes and edges on mesh from myCurve_list
274 SMDS_MeshNode* node_a = NULL;
275 SMDS_MeshNode* node_b = NULL;
276 SMDS_MeshEdge* edge_ab = NULL;
277 SMESHNodes nodesOnEdge;
278 SMESHEdges edgesOnEdge; //backup for group creation
282 nodesOnEdge.push_back(FIRST_NODE);
283 // nodesXxOnEdge.push_back(0.);
284 // _nodeXx[FIRST_NODE] = 0.;
288 double myCurve_start_u = 0.;
289 int nbNodes = law.getNodes(); //law of discretization
290 if (myCurve_list.size()==0)
292 PutData (edge.getName());
295 if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
296 for (int i = 0; i < nbNodes; ++i){
297 u = _Xx(i, law, nbNodes); //u between [0,1]
298 myCurve_u = u*myCurve_tot_len;
301 MESSAGE("myCurve_u -> "<<myCurve_u);
302 MESSAGE("myCurve_tot_len -> "<<myCurve_tot_len);
304 ptOnMyCurve = _getPtOnMyCurve( myCurve_u,
312 node_b = _theMeshDS->AddNode( ptOnMyCurve.X(), ptOnMyCurve.Y(), ptOnMyCurve.Z() );
313 edge_ab = _theMeshDS->AddEdge( node_a, node_b );
314 nodesOnEdge.push_back( node_b );
315 edgesOnEdge.push_back( edge_ab );
316 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
320 edge_ab = _theMeshDS->AddEdge( node_a, LAST_NODE );
321 nodesOnEdge.push_back( LAST_NODE );
322 edgesOnEdge.push_back( edge_ab );
323 _nodesOnEdge[&edge] = nodesOnEdge;
324 _edgesOnEdge[&edge] = edgesOnEdge;
327 if(MYDEBUG) MESSAGE("computeEdgeByAssoc() : end >>>>>>>>");
334 bool SMESH_HexaBlocks::computeEdgeByShortestWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
336 if(MYDEBUG) MESSAGE("computeEdgeByShortestWire() not implemented");
340 bool SMESH_HexaBlocks::computeEdgeByPlanWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
342 if(MYDEBUG) MESSAGE("computeEdgeByPlanWire() not implemented");
346 bool SMESH_HexaBlocks::computeEdgeByIsoWire( HEXA_NS::Edge& edge, HEXA_NS::Law& law)
348 if(MYDEBUG) MESSAGE("computeEdgeByIsoWire() not implemented");
353 bool SMESH_HexaBlocks::computeEdgeBySegment(HEXA_NS::Edge& edge, HEXA_NS::Law& law)
355 if(MYDEBUG) MESSAGE("computeEdgeBySegment() : : begin <<<<<<");
356 ASSERT( _computeVertexOK );
360 HEXA_NS::Vertex* vx0 = NULL;
361 HEXA_NS::Vertex* vx1 = NULL;
363 // way of discretization
364 if (edge.getWay() == true){
365 vx0 = edge.getVertex(0);
366 vx1 = edge.getVertex(1);
368 vx0 = edge.getVertex(1);
369 vx1 = edge.getVertex(0);
373 SMDS_MeshNode* FIRST_NODE = _node[vx0];
374 SMDS_MeshNode* LAST_NODE = _node[vx1];
375 SMDS_MeshNode* node_a = NULL; //new node (to be added)
376 SMDS_MeshNode* node_b = NULL; //new node (to be added)
378 // node and edge creation
379 SMESHNodes nodesOnEdge;
380 SMESHEdges edgesOnEdge;
383 double newNodeX, newNodeY, newNodeZ;
384 SMDS_MeshEdge* newEdge = NULL;
387 nodesOnEdge.push_back(FIRST_NODE);
389 //law of discretization
390 int nbNodes = law.getNodes();
391 if (MYDEBUG) MESSAGE("nbNodes -> "<<nbNodes);
392 for (int i = 0; i < nbNodes; ++i){
393 u = _Xx(i, law, nbNodes);
394 newNodeX = FIRST_NODE->X() + u * ( LAST_NODE->X() - FIRST_NODE->X() );
395 newNodeY = FIRST_NODE->Y() + u * ( LAST_NODE->Y() - FIRST_NODE->Y() );
396 newNodeZ = FIRST_NODE->Z() + u * ( LAST_NODE->Z() - FIRST_NODE->Z() );
397 node_b = _theMeshDS->AddNode(newNodeX, newNodeY, newNodeZ);
398 newEdge = _theMeshDS->AddEdge(node_a, node_b);
399 edgesOnEdge.push_back(newEdge);
400 nodesOnEdge.push_back(node_b);
401 if (_nodeXx.count(node_b) >= 1 ) ASSERT(false);
402 _nodeXx[ node_b ] = u;
403 if(MYDEBUG) MESSAGE("_nodeXx <-"<<u);
406 newEdge = _theMeshDS->AddEdge(node_a, LAST_NODE);
407 nodesOnEdge.push_back(LAST_NODE);
408 edgesOnEdge.push_back(newEdge);
410 _nodesOnEdge[&edge] = nodesOnEdge;
411 _edgesOnEdge[&edge] = edgesOnEdge;
413 if(MYDEBUG) MESSAGE("computeEdgeBySegment() : end >>>>>>>>");
418 // --------------------------------------------------------------
420 // --------------------------------------------------------------
421 std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document* doc )
423 std::map<HEXA_NS::Quad*, bool> quadWays;
424 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
425 std::list<HEXA_NS::Quad*> skinQuad;
426 std::list<HEXA_NS::Quad*> workingQuad;
427 HEXA_NS::Quad* first_q = NULL;
428 HEXA_NS::Quad* q = NULL;
429 HEXA_NS::Edge* e = NULL;
430 HEXA_NS::Vertex *e_0, *e_1 = NULL;
432 // FIRST STEP: eliminate free quad + internal quad
433 int nTotalQuad = doc->countUsedQuad();
434 for (int i=0; i < nTotalQuad; ++i ){
435 q = doc->getUsedQuad(i);
436 switch ( q->getNbrParents() ){ // parent == hexaedron
437 case 0: case 2: quadWays[q] = true; break;
438 case 1: skinQuad.push_back(q); break;
439 default: if ( q->getNbrParents() > 2 ) ASSERT(false);
443 // SECOND STEP: setting edges ways
444 while ( skinQuad.size()>0 ){
445 if(MYDEBUG) MESSAGE("SEARCHING INITIAL QUAD ..." );
446 for ( std::list<HEXA_NS::Quad*>::iterator it = skinQuad.begin(); it != skinQuad.end(); it++ ){
447 _searchInitialQuadWay( *it, e_0, e_1 );
448 if ( e_0 != NULL and e_1 != NULL ){
453 if ( e_0 == NULL and e_1 == NULL ) ASSERT(false);// should never happened,
454 if(MYDEBUG) MESSAGE("INITIAL QUAD FOUND!" );
455 for ( int j=0 ; j < 4 ; ++j ){
457 if ( ((e_0 == e->getVertex(0)) && (e_1 == e->getVertex(1)))
458 || ((e_0 == e->getVertex(1)) && (e_1 == e->getVertex(0))) ){
462 if(MYDEBUG) MESSAGE("INITIAL EDGE WAY FOUND!" );
464 edgeWays[e] = std::make_pair( e_0, e_1 );
465 workingQuad.push_back(q);
467 while ( workingQuad.size() > 0 ){
468 if(MYDEBUG) MESSAGE("COMPUTE QUAD WAY ... ID ="<< q->getId());
469 HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
471 std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
472 while ( localEdgeWays.size() != 4 ){
473 HEXA_NS::Edge* e = q->getEdge(i%4);
474 if ( lastVertex == NULL ){
475 if ( edgeWays.count(e) == 1 ){
477 localEdgeWays[e] = std::make_pair( edgeWays[e].first, edgeWays[e].second );
479 localEdgeWays[e] = std::make_pair( edgeWays[e].second, edgeWays[e].first);
481 firstVertex = localEdgeWays[e].first;
482 lastVertex = localEdgeWays[e].second;
485 HEXA_NS::Vertex* e_0 = e->getVertex(0);
486 HEXA_NS::Vertex* e_1 = e->getVertex(1);
488 if ( lastVertex == e_0 ){
489 firstVertex = e_0; lastVertex = e_1;
490 } else if ( lastVertex == e_1 ){
491 firstVertex = e_1; lastVertex = e_0;
492 } else if ( firstVertex == e_0 ) {
493 firstVertex = e_1; lastVertex = e_0;
494 } else if ( firstVertex == e_1 ) {
495 firstVertex = e_0; lastVertex = e_1;
499 localEdgeWays[e] = std::make_pair( firstVertex, lastVertex );
500 if ( edgeWays.count(e) == 0 ){ // keep current value if present otherwise add it
501 edgeWays[e] = localEdgeWays[e];
508 //add new working quad
509 for ( int i=0 ; i < 4; ++i ){
510 HEXA_NS::Quad* next_q = NULL;
511 HEXA_NS::Edge* e = q->getEdge(i);
512 for ( int j=0 ; j < e->getNbrParents(); ++j ){
513 next_q = e->getParent(j);
517 int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
519 int fromWorkingQuad = std::count( workingQuad.begin(), workingQuad.end(), next_q );
520 if ( fromWorkingQuad == 0 ){
521 workingQuad.push_front( next_q );
528 HEXA_NS::Edge* e0 = q->getEdge(0);
529 HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
531 if ( e0_0 == localEdgeWays[ e0 ].first ){
533 } else if ( e0_0 == localEdgeWays[ e0 ].second ){
538 workingQuad.remove( q );
539 skinQuad.remove( q );
540 q = workingQuad.back();
551 // std::map<HEXA_NS::Quad*, bool> SMESH_HexaBlocks::computeQuadWays( HEXA_NS::Document& doc, std::map<HEXA_NS::Quad*, bool> initQuads )
553 // std::map<HEXA_NS::Quad*, bool> quadWays;
554 // // std::map<HEXA_NS::Edge*, bool> edgeWays;
555 // // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > edgeWays;
556 // std::map<HEXA_NS::Quad*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > workingQuads;
558 // std::list<HEXA_NS::Quad*> skinQuad;
559 // std::list<HEXA_NS::Quad*> notSkinQuad;
560 // // std::list<HEXA_NS::Quad*> workingQuad;
561 // HEXA_NS::Quad* first_q = NULL;
562 // HEXA_NS::Quad* q = NULL;
563 // HEXA_NS::Edge* e = NULL;
564 // HEXA_NS::Vertex *e_0, *e_1 = NULL;
566 // // FIRST STEP: eliminate free quad + internal quad
567 // int nTotalQuad = doc.countQuad();
568 // for (int i=0; i < nTotalQuad; ++i ){
569 // q = doc.getQuad(i);
570 // switch ( q->getNbrParents() ){ // parent == hexaedron
571 // case 0: case 2: quadWays[q] = true; break;
572 // // case 0: case 2: notSkinQuad.push_back(q); break; //CS_TEST
573 // case 1: skinQuad.push_back(q); break;
574 // default: if ( q->getNbrParents() > 2 ) ASSERT(false);
580 // q = first_q = skinQuad.front();
581 // e = q->getUsedEdge(0);
582 // e_0 = e->getVertex(0);
583 // e_1 = e->getVertex(1);
585 // workingQuads[q] = std::make_pair( e_0, e_1 );
587 // while ( workingQuads.size() > 0 ){
588 // MESSAGE("CURRENT QUAD ------>"<< q->getId());
589 // HEXA_NS::Vertex *lastVertex=NULL, *firstVertex = NULL;
592 // std::map<HEXA_NS::Edge*, std::pair<HEXA_NS::Vertex*, HEXA_NS::Vertex*> > localEdgeWays;
593 // while ( localEdgeWays.size() != 4 ){
594 // HEXA_NS::Edge* e = q->getUsedEdge(i%4);
595 // if ( lastVertex == NULL ){
596 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
597 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
599 // if ( (workingQuads[q].first == e_0 and workingQuads[q].second == e_1)
600 // or (workingQuads[q].first == e_1 and workingQuads[q].second == e_0) ){
601 // if ( q == first_q ){
602 // localEdgeWays[e] = std::make_pair( workingQuads[q].first, workingQuads[q].second );
603 // MESSAGE("FIRST QUAD ");
605 // localEdgeWays[e] = std::make_pair( workingQuads[q].second, workingQuads[q].first);
606 // MESSAGE("NOT FIRST QUAD ");
608 // firstVertex = localEdgeWays[e].first;
609 // lastVertex = localEdgeWays[e].second;
612 // HEXA_NS::Vertex* e_0 = e->getVertex(0);
613 // HEXA_NS::Vertex* e_1 = e->getVertex(1);
614 // if ( lastVertex == e_0 ){
615 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
616 // firstVertex = e_0;
618 // } else if ( lastVertex == e_1 ){
619 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
620 // firstVertex = e_1;
622 // } else if ( firstVertex == e_0 ) {
623 // localEdgeWays[e] = std::make_pair( e_1, e_0 );
624 // firstVertex = e_1;
626 // } else if ( firstVertex == e_1 ) {
627 // localEdgeWays[e] = std::make_pair( e_0, e_1 );
628 // firstVertex = e_0;
638 // //add new working quad
639 // for ( int i=0 ; i < 4; ++i ){
640 // HEXA_NS::Quad* next_q = NULL;
641 // HEXA_NS::Edge* e = q->getUsedEdge(i);
642 // MESSAGE("NB PARENTS ="<< e->getNbrParents() );
643 // for ( int j=0 ; j < e->getNbrParents(); ++j ){
644 // next_q = e->getParent(j);
645 // if ( next_q == q ){
648 // int fromSkin = std::count( skinQuad.begin(), skinQuad.end(), next_q );
649 // if (fromSkin != 0){
650 // // int fromWorkingQuad = std::count( workingQuads.begin(), workingQuads.end(), next_q );
651 // int fromWorkingQuad = workingQuads.count( next_q );
652 // // MESSAGE("CHECK QUAD:"<< newWorkingQuad->getId());
653 // if ( fromWorkingQuad == 0 ){
654 // // workingQuads.push_front( next_q );
655 // workingQuads[ next_q ] = localEdgeWays[e];
656 // // MESSAGE("EDGE :"<<e->getId()<<"ADD QUAD :"<< newWorkingQuad->getId());
662 // //setting quad way
663 // HEXA_NS::Edge* e0 = q->getUsedEdge(0);
664 // HEXA_NS::Vertex* e0_0 = e0->getVertex(0);
666 // if ( e0_0 == localEdgeWays[ e0 ].first ){
667 // quadWays[q] = true;
668 // } else if ( e0_0 == localEdgeWays[ e0 ].second ){
669 // quadWays[q] = false;
673 // MESSAGE("quadWays ID ="<< q->getId() << ", WAY = " << quadWays[q] );
675 // // workingQuad.remove( q );
676 // workingQuads.erase( q );
677 // skinQuad.remove( q );
678 // *workingQuads.begin();
679 // q = (*workingQuads.begin()).first;
685 bool SMESH_HexaBlocks::computeQuad( HEXA_NS::Quad& quad, bool way )
689 ok = computeQuadByAssoc(quad, way);
691 ok = computeQuadByFindingGeom(quad, way);
694 ok = computeQuadByLinearApproximation(quad, way);
697 _computeQuadOK = true;
703 bool SMESH_HexaBlocks::computeQuadByAssoc( HEXA_NS::Quad& quad, bool way )
705 // int id = quad.getId();
706 // if ( id != 11 ) return false; //7
708 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
709 MESSAGE("quadID = "<<quad.getId());
711 ASSERT( _computeEdgeOK );
714 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
715 SMESHFaces facesOnQuad;
716 SMDS_MeshFace* newFace = NULL;
717 std::vector<double> xx, yy;
719 // Elements for quad computation
720 SMDS_MeshNode *S1, *S2, *S4, *S3;
722 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
723 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
724 if ( initOk == false ){
728 int nbass = quad.countAssociation ();
731 if(MYDEBUG) MESSAGE("computeQuadByAssoc() : end >>>>>>>>");
735 TopoDS_Shape shapeOrCompound = getFaceShapes ( quad );
738 std::map<SMDS_MeshNode*, gp_Pnt> interpolatedPoints;
739 int iSize = nodesOnQuad.size();
740 int jSize = nodesOnQuad[0].size();
742 S1 = nodesOnQuad[0][0];
743 S2 = nodesOnQuad[iSize-1][0];
744 S4 = nodesOnQuad[0][jSize-1];
745 S3 = nodesOnQuad[iSize-1][jSize-1];
748 for (int j = 1; j < jSize; ++j){
749 for (int i = 1; i < iSize; ++i){
750 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
751 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
752 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
753 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
756 double newNodeX, newNodeY, newNodeZ;
757 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
758 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
759 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
760 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
764 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
765 gp_Pnt newPt = gp_Pnt( newNodeX, newNodeY, newNodeZ );//interpolated point
768 if ( interpolatedPoints.count(n1) > 0 ){
769 pt1 = interpolatedPoints[n1];
771 pt1 = gp_Pnt( n1->X(), n1->Y(), n1->Z() );
773 if ( interpolatedPoints.count(n3) > 0 ){
774 pt3 = interpolatedPoints[n3];
776 pt3 = gp_Pnt( n3->X(), n3->Y(), n3->Z() );
778 gp_Vec vec1( newPt, pt1 );
779 gp_Vec vec2( newPt, pt3 );
781 gp_Pnt ptOnShape = _intersect(newPt, vec1, vec2, shapeOrCompound);
782 newNodeX = ptOnShape.X();
783 newNodeY = ptOnShape.Y();
784 newNodeZ = ptOnShape.Z();
785 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
786 nodesOnQuad[i][j] = n4;
787 interpolatedPoints[ n4 ] = newPt;
790 MESSAGE("u parameter is "<<u);
791 MESSAGE("v parameter is "<<v);
792 MESSAGE("point interpolated ("<<newPt.X()<<","<<newPt.Y()<<","<<newPt.Z()<<" )");
793 MESSAGE("point on shape ("<<newNodeX<<","<<newNodeY<<","<<newNodeZ<<" )");
798 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
799 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
800 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
801 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
805 if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
806 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
808 if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
809 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
811 facesOnQuad.push_back(newFace);
814 _quadNodes[ &quad ] = nodesOnQuad;
815 _facesOnQuad[&quad] = facesOnQuad;
817 if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
822 bool SMESH_HexaBlocks::computeQuadByFindingGeom( HEXA_NS::Quad& quad, bool way )
824 if(MYDEBUG) MESSAGE("computeQuadByFindingGeom() not implemented");
828 bool SMESH_HexaBlocks::_computeQuadInit(
830 ArrayOfSMESHNodes& nodesOnQuad,
831 std::vector<double>& xx, std::vector<double>& yy)
833 if(MYDEBUG) MESSAGE("_computeQuadInit() : begin ---------------");
836 SMDS_MeshNode *S1, *S2, *S4, *S3;
837 HEXA_NS::Edge *eh, *eb, *eg, *ed;
838 HEXA_NS::Edge *e1, *e2, *e3, *e4;
839 HEXA_NS::Vertex *e1_0, *e1_1, *e2_0, *e2_1, *e3_0, *e3_1, *e4_0, *e4_1;
841 e1 = quad.getEdge(0);
842 e2 = quad.getEdge(1);
843 e3 = quad.getEdge(2);
844 e4 = quad.getEdge(3);
846 e1_0 = e1->getVertex(0); e1_1 = e1->getVertex(1);
847 e2_0 = e2->getVertex(0); e2_1 = e2->getVertex(1);
848 e3_0 = e3->getVertex(0); e3_1 = e3->getVertex(1);
849 e4_0 = e4->getVertex(0); e4_1 = e4->getVertex(1);
852 S1 = _node[e1_0]; S2 = _node[e1_1];
858 } else if ( e1_0 == e2_1 ){
861 } else if ( e1_0 == e4_0 ){
864 } else if ( e1_0 == e4_1 ){
871 if ( S4 == _node[e3_0] ){
873 } else if ( S4 == _node[e3_1] ){
879 SMESHNodes hNodes = _nodesOnEdge[eh];
880 SMESHNodes bNodes = _nodesOnEdge[eb];
881 SMESHNodes gNodes = _nodesOnEdge[eg];
882 SMESHNodes dNodes = _nodesOnEdge[ed];
883 nodesOnQuad.resize( bNodes.size(), SMESHNodes(gNodes.size(), static_cast<SMDS_MeshNode*>(NULL)) );
887 // int &b_i = i, &h_i = i, &g_j = j, &d_j = j;
888 int *b_i = &i, *h_i = &i, *g_j = &j, *d_j = &j;
889 bool uWay = true, vWay = true;
891 if ( bNodes[0] != S1 ){
894 ASSERT( bNodes[bNodes.size()-1] == S1 );
896 ASSERT( bNodes[0] == S1);
898 if ( hNodes[0] != S4 ){
901 ASSERT( hNodes[0] == S4 );
903 if ( gNodes[0] != S1 ){
907 ASSERT( gNodes[0] == S1 );
909 if ( dNodes[0] != S2 ){
912 ASSERT( dNodes[0] == S2 );
917 for (i = 0, _i = bNodes.size()-1; i < bNodes.size(); ++i, --_i){
918 nodesOnQuad[i][0] = bNodes[*b_i];
919 nodesOnQuad[i][gNodes.size()-1 ] = hNodes[*h_i];
921 u = _nodeXx[ bNodes[*b_i] ];
928 if ( S1 != nodesOnQuad[0][0] ){
929 if(MYDEBUG) MESSAGE("ZZZZZZZZZZZZZZZZ quadID = "<<quad.getId());
931 // ASSERT( S1 == nodesOnQuad[0][0] );
935 for (j = 0, _j = gNodes.size()-1; j < gNodes.size(); ++j, --_j){
936 nodesOnQuad[0][j] = gNodes[*g_j];
937 if ( S1 != nodesOnQuad[0][0] ){
938 if(MYDEBUG) MESSAGE("XXXXXXXXXXXXXXXX quadID = "<<quad.getId());
940 // ASSERT( S1 == nodesOnQuad[0][0] );
941 nodesOnQuad[bNodes.size()-1][j] = dNodes[*d_j];
942 v = _nodeXx[ gNodes[*g_j] ];
951 int iSize = nodesOnQuad.size();
952 int jSize = nodesOnQuad[0].size();
953 ASSERT( iSize = bNodes.size() );
954 ASSERT( jSize = gNodes.size() );
956 // ASSERT( S1 == nodesOnQuad[0][0] );
957 // ASSERT( S2 == nodesOnQuad[iSize-1][0]);
958 // ASSERT( S4 == nodesOnQuad[0][jSize-1]);
959 // ASSERT( S3 == nodesOnQuad[iSize-1][jSize-1]);
965 bool SMESH_HexaBlocks::computeQuadByLinearApproximation( HEXA_NS::Quad& quad, bool way )
967 // int id = quad.getId();
968 // if ( quad.getId() != 66 ) return false; //7, 41
970 MESSAGE("computeQuadByLinearApproximation() : : begin <<<<<<");
971 MESSAGE("quadID = "<<quad.getId());
973 ASSERT( _computeEdgeOK );
976 ArrayOfSMESHNodes nodesOnQuad; // nodes in this quad ( to be added on the mesh )
977 SMESHFaces facesOnQuad;
978 SMDS_MeshFace* newFace = NULL;
979 std::vector<double> xx, yy;
981 // Elements for quad computation
982 SMDS_MeshNode *S1, *S2, *S4, *S3;
984 // bool initOk = _computeQuadInit( quad, eh, eb, eg, ed, S1, S2, S3, S4 );
985 bool initOk = _computeQuadInit( quad, nodesOnQuad, xx, yy );
986 if ( initOk == false ){
990 int iSize = nodesOnQuad.size();
991 int jSize = nodesOnQuad[0].size();
993 S1 = nodesOnQuad[0][0];
994 // S2 = nodesOnQuad[bNodes.size()-1][0];
995 S2 = nodesOnQuad[iSize-1][0];
996 S4 = nodesOnQuad[0][jSize-1];
997 S3 = nodesOnQuad[iSize-1][jSize-1];
999 for (int j = 1; j < jSize; ++j){
1000 for (int i = 1; i < iSize; ++i){
1001 SMDS_MeshNode* n1 = nodesOnQuad[i-1][j];
1002 SMDS_MeshNode* n2 = nodesOnQuad[i-1][j-1];
1003 SMDS_MeshNode* n3 = nodesOnQuad[i][j-1];
1004 SMDS_MeshNode* n4 = nodesOnQuad[i][j];
1007 double newNodeX, newNodeY, newNodeZ;
1008 SMDS_MeshNode* Ph = nodesOnQuad[i][jSize-1]; //dNodes[h_i];
1009 SMDS_MeshNode* Pb = nodesOnQuad[i][0]; //bNodes[b_i];
1010 SMDS_MeshNode* Pg = nodesOnQuad[0][j]; //gNodes[g_j];
1011 SMDS_MeshNode* Pd = nodesOnQuad[iSize-1][j]; //dNodes[d_j];
1015 _nodeInterpolationUV(u, v, Pg, Pd, Ph, Pb, S1, S2, S3, S4, newNodeX, newNodeY, newNodeZ);
1016 n4 = _theMeshDS->AddNode( newNodeX, newNodeY, newNodeZ );
1017 nodesOnQuad[i][j] = n4;
1021 MESSAGE("n1 (" << n1->X() << "," << n1->Y() << "," << n1->Z() << ")");
1022 MESSAGE("n2 (" << n2->X() << "," << n2->Y() << "," << n2->Z() << ")");
1023 MESSAGE("n4 (" << n4->X() << "," << n4->Y() << "," << n4->Z() << ")");
1024 MESSAGE("n3 (" << n3->X() << "," << n3->Y() << "," << n3->Z() << ")");
1028 if (MYDEBUG) MESSAGE("AddFace( n1, n2, n3, n4 )");
1029 newFace = _theMeshDS->AddFace( n1, n2, n3, n4 );
1031 if (MYDEBUG) MESSAGE("AddFace( n4, n3, n2, n1 )");
1032 newFace = _theMeshDS->AddFace( n4, n3, n2, n1 );
1034 facesOnQuad.push_back(newFace);
1037 _quadNodes[ &quad ] = nodesOnQuad;
1038 _facesOnQuad[&quad] = facesOnQuad;
1040 if(MYDEBUG) MESSAGE("computeQuadByLinearApproximation() : end >>>>>>>>");
1045 // --------------------------------------------------------------
1047 // --------------------------------------------------------------
1048 bool SMESH_HexaBlocks::computeHexa( HEXA_NS::Document* doc )
1050 if(MYDEBUG) MESSAGE("computeHexa() : : begin <<<<<<");
1053 SMESH_MesherHelper aHelper(*_theMesh);
1054 TopoDS_Shape shape = _theMesh->GetShapeToMesh();
1055 aHelper.SetSubShape( shape );
1056 aHelper.SetElementsOnShape( true );
1058 SMESH_Gen* gen = _theMesh->GetGen();
1059 SMESH_HexaFromSkin_3D algo( gen->GetANewId(), 0, gen, doc );
1060 algo.InitComputeError();
1062 ok = algo.Compute( *_theMesh, &aHelper, _volumesOnHexa, _node );
1064 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D error!!! ");
1067 MESSAGE("SMESH_HexaFromSkin_3D.comment = "<<algo.GetComputeError()->myComment);
1068 MESSAGE("computeHexa() : end >>>>>>>>");
1075 // --------------------------------------------------------------
1076 // Document computing
1077 // --------------------------------------------------------------
1078 bool SMESH_HexaBlocks::computeDoc( HEXA_NS::Document* doc )
1080 if(MYDEBUG) MESSAGE("computeDoc() : : begin <<<<<<");
1083 // A) Vertex computation
1086 int nVertex = doc->countUsedVertex();
1087 HEXA_NS::Vertex* vertex = NULL;
1089 for (int j=0; j <nVertex; ++j ){ //Computing each vertex of the document
1090 vertex = doc->getUsedVertex(j);
1091 ok = computeVertex(*vertex);
1094 // B) Edges computation
1096 HEXA_NS::Propagation* propa = NULL;
1097 HEXA_NS::Law* law = NULL;
1098 HEXA_NS::Edges edges;
1100 nbPropa = doc->countPropagation();
1101 for (int j=0; j < nbPropa; ++j ){//Computing each edge's propagations of the document
1102 propa = doc->getPropagation(j);
1103 edges = propa->getEdges();
1104 law = propa->getLaw();
1107 law = doc->getLaw(0); // default law
1109 for( HEXA_NS::Edges::const_iterator iter = edges.begin();
1110 iter != edges.end();
1112 ok = computeEdge(**iter, *law);
1115 // C) Quad computation
1116 std::map<HEXA_NS::Quad*, bool> quadWays = computeQuadWays(doc);
1117 int nQuad = doc->countUsedQuad();
1118 HEXA_NS::Quad* q = NULL;
1119 for (int j=0; j <nQuad; ++j ){ //Computing each quad of the document
1120 q = doc->getUsedQuad(j);
1121 int id = q->getId();
1122 if ( quadWays.count(q) > 0 )
1123 ok = computeQuad( *q, quadWays[q] );
1125 if(MYDEBUG) MESSAGE("NO QUAD WAY ID = "<<id);
1129 // D) Hexa computation: Calling HexaFromSkin algo
1130 ok = computeHexa(doc);
1132 if(MYDEBUG) MESSAGE("computeDoc() : end >>>>>>>>");
1138 void SMESH_HexaBlocks::buildGroups(HEXA_NS::Document* doc)
1141 MESSAGE("_addGroups() : : begin <<<<<<");
1142 MESSAGE("_addGroups() : : nb. hexas= " << doc->countUsedHexa());
1143 MESSAGE("_addGroups() : : nb. quads= " << doc->countUsedQuad());
1144 MESSAGE("_addGroups() : : nb. edges= " << doc->countUsedEdge());
1145 MESSAGE("_addGroups() : : nb. nodes= " << doc->countUsedVertex());
1147 // Looping on each groups of the document
1148 for ( int i=0; i < doc->countGroup(); i++ ){
1149 _fillGroup( doc->getGroup(i) );
1152 if(MYDEBUG) MESSAGE("_addGroups() : end >>>>>>>>");
1155 // --------------------------------------------------------------
1157 // --------------------------------------------------------------
1158 double SMESH_HexaBlocks::_Xx( double i, HEXA_NS::Law law, double nbNodes) //, double pos0 )
1163 HEXA_NS::KindLaw k = law.getKind();
1164 double coeff = law.getCoefficient();
1166 case HEXA_NS::Uniform:
1167 result = (i+1)/(nbNodes+1);
1168 if(MYDEBUG) MESSAGE( "_Xx():" << " Uniform u("<<i<< ")"<< " = " << result);
1170 case HEXA_NS::Arithmetic:
1171 u0 = 1./(nbNodes + 1.) - (coeff*nbNodes)/2.;
1173 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Arithmetic discretization : check coefficient")));
1177 result = (i + 1.)*u0 + coeff*i*(i+1.)/2.;
1179 if(MYDEBUG) MESSAGE( "_Xx():" << " Arithmetic u("<<i<< ")"<< " = " << result);
1181 case HEXA_NS::Geometric:
1182 u0 = (1.-coeff)/(1.-pow(coeff, nbNodes + 1) ) ;
1184 if ( u0 <= 0 ) throw (SALOME_Exception(LOCALIZED("Geometric discretization : check coefficient")));
1188 result = u0 * (1.- pow(coeff, i + 1) )/(1.-coeff) ;
1190 if(MYDEBUG) MESSAGE( "_Xx():" << " Geometric u("<<i<< ")"<< " = " << result);
1197 // ============================================================= _edgeLength
1198 double SMESH_HexaBlocks::_edgeLength(const TopoDS_Edge & E)
1200 if(MYDEBUG) MESSAGE("_edgeLength() : : begin <<<<<<");
1201 double UMin = 0, UMax = 0;
1202 if (BRep_Tool::Degenerated(E))
1205 Handle(Geom_Curve) C = BRep_Tool::Curve(E, L, UMin, UMax);
1206 GeomAdaptor_Curve AdaptCurve(C);
1207 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, UMin, UMax);
1208 if(MYDEBUG) MESSAGE("_edgeLength() : end >>>>>>>>");
1211 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1212 // ============================================================== _buildMyCurve
1213 // === construire ma courbe a moi
1214 void SMESH_HexaBlocks::_buildMyCurve(
1215 const gp_Pnt& myCurve_pt_start, //IN
1216 const gp_Pnt& myCurve_pt_end, //IN
1217 std::list< BRepAdaptor_Curve* >& myCurve_list, //INOUT
1218 double& myCurve_tot_len, //INOUT
1219 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//INOUT
1220 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //INOUT
1221 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //INOUT
1222 HEXA_NS::Edge& edge) // For error diagnostic
1224 if(MYDEBUG) MESSAGE("_buildMyCurve() : : begin <<<<<<");
1225 bool current_way = true;
1226 myCurve_tot_len = 0.;
1227 BRepAdaptor_Curve* thePreviousCurve = NULL;
1228 BRepAdaptor_Curve* theCurve = NULL;
1230 gp_Pnt curv_start, curv_end;
1231 gp_Pnt p_curv_start , p_curv_end;
1233 double curv_length = 0;
1235 int nbass = edge.countAssociation();
1236 for (int nro = 0 ; nro < nbass ; ++nro)
1238 theCurve = make_curve (curv_start, curv_end, curv_length, u_start,
1240 if (theCurve != NULL)
1243 if ( thePreviousCurve == NULL )
1245 // setting current_way and first curve way
1246 if ( myCurve_pt_start.IsEqual(curv_start, HEXA_EPS) )
1251 else if ( myCurve_pt_start.IsEqual(curv_end, HEXA_EPS) )
1256 else if ( myCurve_pt_end.IsEqual(curv_end, HEXA_EPS) )
1258 current_way = false;
1261 else if ( myCurve_pt_end.IsEqual(curv_start, HEXA_EPS) )
1263 current_way = false;
1269 MESSAGE("SOMETHING WRONG on edge association... Bad script?");
1271 throw (SALOME_Exception (LOCALIZED("Edge association : check association parameters ( first, last ) between HEXA model and CAO")));
1275 // it is not the first or last curve.
1276 // ways are calculated between previous and new one.
1279 if ( p_curv_end.IsEqual( curv_end, HEXA_EPS )
1280 || p_curv_start.IsEqual( curv_start, HEXA_EPS ) )
1282 sens = NOT myCurve_ways[thePreviousCurve];// opposite WAY
1284 else if ( p_curv_end.IsEqual( curv_start, HEXA_EPS )
1285 || p_curv_start.IsEqual( curv_end, HEXA_EPS ) )
1287 sens = myCurve_ways[thePreviousCurve];// same WAY
1292 MESSAGE("SOMETHING WRONG on edge association... bad script?");
1295 throw (SALOME_Exception(LOCALIZED("Edge association : Check association parameters ( first, last ) between HEXA model and CAO")));
1299 myCurve_list.push_back( theCurve );
1300 myCurve_ways [theCurve] = sens;
1301 myCurve_starts [theCurve] = u_start;
1302 myCurve_lengths[theCurve] = curv_length;
1303 myCurve_tot_len += curv_length;
1305 thePreviousCurve = theCurve;
1306 p_curv_start = curv_start;
1307 p_curv_end = curv_end;
1308 }//if ( theCurveLength > 0 )
1312 if ( NOT current_way ){
1313 std::list< BRepAdaptor_Curve* > tmp( myCurve_list.size() );
1314 std::copy( myCurve_list.rbegin(), myCurve_list.rend(), tmp.begin() );
1319 MESSAGE("current_way was :" << current_way);
1320 MESSAGE("_buildMyCurve() : end >>>>>>>>");
1327 gp_Pnt SMESH_HexaBlocks::_getPtOnMyCurve(
1328 const double& myCurve_u, //IN
1329 std::map< BRepAdaptor_Curve*, bool>& myCurve_ways, //IN
1330 std::map< BRepAdaptor_Curve*, double>& myCurve_lengths,//IN
1331 std::map< BRepAdaptor_Curve*, double>& myCurve_starts, //IN
1332 std::list< BRepAdaptor_Curve* >& myCurve_list, //INOUT
1333 double& myCurve_start ) //INOUT
1334 // std::map< BRepAdaptor_Curve*, double>& myCurve_firsts,
1335 // std::map< BRepAdaptor_Curve*, double>& myCurve_lasts,
1337 if(MYDEBUG) MESSAGE("_getPtOnMyCurve() : : begin <<<<<<");
1340 // looking for curve which contains parameter myCurve_u
1341 BRepAdaptor_Curve* curve = myCurve_list.front();
1342 double curve_start = myCurve_start;
1343 double curve_end = curve_start + myCurve_lengths[curve];
1345 GCPnts_AbscissaPoint discret;
1348 MESSAGE("looking for curve = "<<(long) curve);
1349 MESSAGE("looking for curve: curve_u = "<<myCurve_u);
1350 MESSAGE("looking for curve: curve_start = "<<curve_start);
1351 MESSAGE("looking for curve: curve_end = "<<curve_end);
1352 MESSAGE("looking for curve: curve_lenght = "<<myCurve_lengths[curve]);
1353 MESSAGE("looking for curve: curve.size _lenght= "<<myCurve_list.size());
1355 while ( not ( (myCurve_u >= curve_start) and (myCurve_u <= curve_end) ) ) {
1357 ASSERT( myCurve_list.size() != 0 );
1358 myCurve_list.pop_front();
1359 curve = myCurve_list.front();
1360 curve_start = curve_end;
1361 curve_end = curve_start + myCurve_lengths[curve];
1363 MESSAGE("go next curve: curve_lenght = "<<myCurve_lengths[curve]);
1364 MESSAGE("go next curve: curve_start = "<<curve_start);
1365 MESSAGE("go next curve: curve_end = "<<curve_end);
1366 MESSAGE("go next curve: myCurve_u = "<<myCurve_u);
1369 myCurve_start = curve_start;
1372 if ( myCurve_ways[curve] ){
1373 discret = GCPnts_AbscissaPoint( *curve, (myCurve_u - curve_start), myCurve_starts[curve] );
1375 discret = GCPnts_AbscissaPoint( *curve, myCurve_lengths[curve]- (myCurve_u - curve_start), myCurve_starts[curve] );
1377 // PutData (discret);
1378 ASSERT(discret.IsDone());
1379 curve_u = discret.Parameter();
1380 ptOnMyCurve = curve->Value( curve_u );
1383 MESSAGE("curve found!");
1384 MESSAGE("curve_u = "<< curve_u);
1385 MESSAGE("curve way = "<< myCurve_ways[curve]);
1386 MESSAGE("_getPtOnMyCurve() : end >>>>>>>>");
1394 void SMESH_HexaBlocks::_nodeInterpolationUV(double u, double v,
1395 SMDS_MeshNode* Pg, SMDS_MeshNode* Pd, SMDS_MeshNode* Ph, SMDS_MeshNode* Pb,
1396 SMDS_MeshNode* S1, SMDS_MeshNode* S2, SMDS_MeshNode* S3, SMDS_MeshNode* S4,
1397 double& xOut, double& yOut, double& zOut )
1400 MESSAGE("_nodeInterpolationUV() IN:");
1401 MESSAGE("u ( "<< u <<" )");
1402 MESSAGE("v ( "<< v <<" )");
1404 MESSAGE("S1 (" << S1->X() << "," << S1->Y() << "," << S1->Z() << ")");
1405 MESSAGE("S2 (" << S2->X() << "," << S2->Y() << "," << S2->Z() << ")");
1406 MESSAGE("S4 (" << S4->X() << "," << S4->Y() << "," << S4->Z() << ")");
1407 MESSAGE("S3 (" << S3->X() << "," << S3->Y() << "," << S3->Z() << ")");
1409 MESSAGE("Pg (" << Pg->X() << "," << Pg->Y() << "," << Pg->Z() << ")");
1410 MESSAGE("Pd (" << Pd->X() << "," << Pd->Y() << "," << Pd->Z() << ")");
1411 MESSAGE("Ph (" << Ph->X() << "," << Ph->Y() << "," << Ph->Z() << ")");
1412 MESSAGE("Pb (" << Pb->X() << "," << Pb->Y() << "," << Pb->Z() << ")");
1415 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();
1416 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();
1417 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();
1420 MESSAGE("_nodeInterpolationUV() OUT("<<xOut<<","<<yOut<<","<<zOut<<" )");
1424 // =========================================================== getFaceShapes
1425 TopoDS_Shape SMESH_HexaBlocks::getFaceShapes (Hex::Quad& quad)
1427 int nbass = quad.countAssociation ();
1428 Hex::FaceShape* face = quad.getAssociation (0);
1430 return face->getShape ();
1432 TopoDS_Compound compound;
1433 BRep_Builder builder;
1434 builder.MakeCompound (compound);
1436 for (int nro=0 ; nro < nbass ; nro++)
1438 face = quad.getAssociation (nro);
1439 TopoDS_Shape shape = face->getShape ();
1440 builder.Add (compound, shape);
1447 // ================================================== carre
1448 inline double carre (double val)
1452 // ================================================== dist2
1453 inline double dist2 (const gp_Pnt& pt1, const gp_Pnt& pt2)
1455 double dist = carre (pt2.X()-pt1.X()) + carre (pt2.Y()-pt1.Y())
1456 + carre (pt2.Z()-pt1.Z());
1459 // ================================================== _intersect
1460 gp_Pnt SMESH_HexaBlocks::_intersect( const gp_Pnt& Pt,
1461 const gp_Vec& u, const gp_Vec& v,
1462 const TopoDS_Shape& shape,
1467 gp_Vec normale = u^v;
1468 gp_Dir dir(normale);
1469 gp_Lin li( Pt, dir );
1471 Standard_Real s = -Precision::Infinite();
1472 Standard_Real e = +Precision::Infinite();
1474 IntCurvesFace_ShapeIntersector inter;
1475 inter.Load(shape, tol);
1476 // inter.Load(S, tol);
1477 inter.Perform(li, s, e);//inter.PerformNearest(li, s, e);
1479 /*********************************************** Abu 2011-11-04 */
1480 if ( inter.IsDone() )
1482 result = inter.Pnt(1);//first
1483 int nbrpts = inter.NbPnt();
1486 double d0 = dist2 (result, Pt);
1487 for (int i=2; i <= inter.NbPnt(); ++i )
1489 double d1 = dist2 (Pt, inter.Pnt(i));
1493 result = inter.Pnt (i);
1497 /********************************************** Abu 2011-11-04 (getEnd()) */
1499 MESSAGE("_intersect() : OK");
1500 for ( int i=1; i <= inter.NbPnt(); ++i ){
1501 gp_Pnt tmp = inter.Pnt(i);
1502 MESSAGE("_intersect() : pnt("<<i<<") = ("<<tmp.X()<<","<<tmp.Y()<<","<<tmp.Z()<<" )");
1507 if(MYDEBUG) MESSAGE("_intersect() : KO");
1516 // parameters q : IN, v0: INOUT, v1: INOUT
1517 void SMESH_HexaBlocks::_searchInitialQuadWay( HEXA_NS::Quad* q, HEXA_NS::Vertex*& v0, HEXA_NS::Vertex*& v1 )
1519 if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : begin");
1520 v0 = NULL; v1 = NULL;
1521 if ( q->getNbrParents() != 1 ) return; // q must be a skin quad
1523 HEXA_NS::Vertex* qA = q->getVertex(0);
1524 HEXA_NS::Vertex* qB = q->getVertex(1);
1525 HEXA_NS::Vertex* qC = q->getVertex(2);
1526 HEXA_NS::Vertex* qD = q->getVertex(3);
1528 // searching for vertex on opposed quad
1529 HEXA_NS::Vertex *qAA = NULL, *qBB = NULL, *qCC = NULL, *qDD = NULL;
1530 HEXA_NS::Hexa* h = q->getParent(0);
1531 for( int i=0; i < h->countEdge(); ++i ){
1532 HEXA_NS::Edge* e = h->getEdge(i);
1533 HEXA_NS::Vertex* e0 = e->getVertex(0);
1534 HEXA_NS::Vertex* e1 = e->getVertex(1);
1536 if ( e0 == qA and e1 != qB and e1 != qC and e1 != qD ){
1538 } else if ( e1 == qA and e0 != qB and e0 != qC and e0 != qD ){
1540 } else if ( e0 == qB and e1 != qA and e1 != qC and e1 != qD ){
1542 } else if ( e1 == qB and e0 != qA and e0 != qC and e0 != qD ){
1544 } else if ( e0 == qC and e1 != qA and e1 != qB and e1 != qD ){
1546 } else if ( e1 == qC and e0 != qA and e0 != qB and e0 != qD ){
1548 } else if ( e0 == qD and e1 != qA and e1 != qB and e1 != qC ){
1550 } else if ( e1 == qD and e0 != qA and e0 != qB and e0 != qC ){
1555 // working on final value ( point on CAO ), not on model
1556 SMDS_MeshNode *nA = _node[qA], *nAA = _node[qAA];
1557 SMDS_MeshNode *nB = _node[qB], *nBB = _node[qBB];
1558 SMDS_MeshNode *nC = _node[qC], *nCC = _node[qCC];
1559 SMDS_MeshNode *nD = _node[qD], *nDD = _node[qDD];
1561 gp_Pnt pA( nA->X(), nA->Y(), nA->Z() );
1562 gp_Pnt pB( nB->X(), nB->Y(), nB->Z() );
1563 gp_Pnt pC( nC->X(), nC->Y(), nC->Z() );
1564 gp_Pnt pD( nD->X(), nD->Y(), nD->Z() );
1566 gp_Pnt pAA( nAA->X(), nAA->Y(), nAA->Z() );
1567 gp_Pnt pBB( nBB->X(), nBB->Y(), nBB->Z() );
1568 gp_Pnt pCC( nCC->X(), nCC->Y(), nCC->Z() );
1569 gp_Pnt pDD( nDD->X(), nDD->Y(), nDD->Z() );
1571 gp_Vec AB( pA, pB );
1572 gp_Vec AC( pA, pC );
1573 gp_Vec normP = AB^AC;
1574 gp_Dir dirP( normP );
1576 // building plane for point projection
1577 gp_Pln plnP( gp_Pnt(nA->X(), nA->Y(), nA->Z()), dirP);
1578 TopoDS_Shape sPlnP = BRepBuilderAPI_MakeFace(plnP).Face();
1580 // PAAA is the result of PAA projection
1581 gp_Pnt pAAA = _intersect( pAA, AB, AC, sPlnP );
1582 gp_Pnt pBBB = _intersect( pBB, AB, AC, sPlnP );
1583 gp_Pnt pCCC = _intersect( pCC, AB, AC, sPlnP );
1584 gp_Pnt pDDD = _intersect( pDD, AB, AC, sPlnP );
1586 gp_Dir AA( gp_Vec(pAA, pAAA) );
1587 gp_Dir BB( gp_Vec(pBB, pBBB) );
1588 gp_Dir CC( gp_Vec(pCC, pCCC) );
1589 gp_Dir DD( gp_Vec(pDD, pDDD) );
1591 // eventually, we are able to know if the input quad is a good client!
1592 // exit the fonction otherwise
1593 if ( AA.IsOpposite(BB, HEXA_QUAD_WAY) ) return;
1594 if ( BB.IsOpposite(CC, HEXA_QUAD_WAY) ) return;
1595 if ( CC.IsOpposite(DD, HEXA_QUAD_WAY) ) return;
1597 // ok, give the input quad the good orientation by
1599 if ( !dirP.IsOpposite(AA, HEXA_QUAD_WAY) ) { //OK
1605 if(MYDEBUG) MESSAGE("_searchInitialQuadWay() : end");
1608 SMESH_Group* SMESH_HexaBlocks::_createGroup(HEXA_NS::Group& grHex)
1610 if(MYDEBUG) MESSAGE("_createGroup() : : begin <<<<<<");
1612 std::string aGrName = grHex.getName();
1613 HEXA_NS::EnumGroup grHexKind = grHex.getKind();
1615 if(MYDEBUG) MESSAGE("aGrName"<<aGrName);
1617 SMDSAbs_ElementType aGrType;
1618 switch ( grHexKind ){
1619 case HEXA_NS::HexaCell : aGrType = SMDSAbs_Volume; break;
1620 case HEXA_NS::QuadCell : aGrType = SMDSAbs_Face ; break;
1621 case HEXA_NS::EdgeCell : aGrType = SMDSAbs_Edge ; break;
1622 case HEXA_NS::HexaNode : aGrType = SMDSAbs_Node ; break;
1623 case HEXA_NS::QuadNode : aGrType = SMDSAbs_Node ; break;
1624 case HEXA_NS::EdgeNode : aGrType = SMDSAbs_Node ; break;
1625 case HEXA_NS::VertexNode : aGrType = SMDSAbs_Node ; break;
1626 default : ASSERT(false);
1630 SMESH_Group* aGr = _theMesh->AddGroup(aGrType, aGrName.c_str(), aId);
1632 if(MYDEBUG) MESSAGE("_createGroup() : end >>>>>>>>");
1636 void SMESH_HexaBlocks::_fillGroup(HEXA_NS::Group* grHex )
1638 if(MYDEBUG) MESSAGE("_fillGroup() : : begin <<<<<<");
1640 SMESH_Group* aGr = _createGroup( *grHex );
1641 HEXA_NS::EltBase* grHexElt = NULL;
1642 HEXA_NS::EnumGroup grHexKind = grHex->getKind();
1643 int grHexNbElt = grHex->countElement();
1645 if(MYDEBUG) MESSAGE("_fillGroup() : kind = " << grHexKind);
1646 if(MYDEBUG) MESSAGE("_fillGroup() : count= " << grHexNbElt);
1648 // A)Looking for elements ID
1649 std::vector<const SMDS_MeshElement*> aGrEltIDs;
1651 for ( int n=0; n<grHexNbElt; ++n ){
1652 grHexElt = grHex->getElement(n);
1654 switch ( grHexKind ){
1655 case HEXA_NS::HexaCell:
1657 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1658 if ( _volumesOnHexa.count(h)>0 ){
1659 SMESHVolumes volumes = _volumesOnHexa[h];
1660 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1661 aGrEltIDs.push_back(*aVolume);
1664 if(MYDEBUG) MESSAGE("GROUP OF VOLUME: volume for hexa (id = "<<h->getId()<<") not found");
1668 case HEXA_NS::QuadCell:
1670 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1671 if ( _facesOnQuad.count(q)>0 ){
1672 SMESHFaces faces = _facesOnQuad[q];
1673 for ( SMESHFaces::iterator aFace = faces.begin(); aFace != faces.end(); ++aFace ){
1674 aGrEltIDs.push_back(*aFace);
1677 if(MYDEBUG) MESSAGE("GROUP OF FACE: face for quad (id = "<<q->getId()<<") not found");
1681 case HEXA_NS::EdgeCell:
1683 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1684 if ( _edgesOnEdge.count(e)>0 ){
1685 SMESHEdges edges = _edgesOnEdge[e];
1686 for ( SMESHEdges::iterator anEdge = edges.begin(); anEdge != edges.end(); ++anEdge ){
1687 aGrEltIDs.push_back(*anEdge);
1690 if(MYDEBUG) MESSAGE("GROUP OF Edge: edge for edge (id = "<<e->getId()<<") not found");
1694 case HEXA_NS::HexaNode:
1696 HEXA_NS::Hexa* h = reinterpret_cast<HEXA_NS::Hexa*>(grHexElt);
1697 if ( _volumesOnHexa.count(h)>0 ){
1698 SMESHVolumes volumes = _volumesOnHexa[h];
1699 for ( SMESHVolumes::iterator aVolume = volumes.begin(); aVolume != volumes.end(); ++aVolume ){
1700 SMDS_ElemIteratorPtr aNodeIter = (*aVolume)->nodesIterator();
1701 while( aNodeIter->more() ){
1702 const SMDS_MeshNode* aNode =
1703 dynamic_cast<const SMDS_MeshNode*>( aNodeIter->next() );
1705 aGrEltIDs.push_back(aNode);
1710 if(MYDEBUG) MESSAGE("GROUP OF HEXA NODES: nodes on hexa (id = "<<h->getId()<<") not found");
1714 case HEXA_NS::QuadNode:
1716 HEXA_NS::Quad* q = reinterpret_cast<HEXA_NS::Quad*>(grHexElt);
1717 if ( _quadNodes.count(q)>0 ){
1718 ArrayOfSMESHNodes nodesOnQuad = _quadNodes[q];
1719 for ( ArrayOfSMESHNodes::iterator nodes = nodesOnQuad.begin(); nodes != nodesOnQuad.end(); ++nodes){
1720 for ( SMESHNodes::iterator aNode = nodes->begin(); aNode != nodes->end(); ++aNode){
1721 aGrEltIDs.push_back(*aNode);
1725 if(MYDEBUG) MESSAGE("GROUP OF QUAD NODES: nodes on quad (id = "<<q->getId()<<") not found");
1729 case HEXA_NS::EdgeNode:
1731 HEXA_NS::Edge* e = reinterpret_cast<HEXA_NS::Edge*>(grHexElt);
1732 if ( _nodesOnEdge.count(e)>0 ){
1733 SMESHNodes nodes = _nodesOnEdge[e];
1734 for ( SMESHNodes::iterator aNode = nodes.begin(); aNode != nodes.end(); ++aNode){
1735 aGrEltIDs.push_back(*aNode);
1738 if(MYDEBUG) MESSAGE("GROUP OF EDGE NODES: nodes on edge (id = "<<e->getId()<<") not found");
1742 case HEXA_NS::VertexNode:
1744 HEXA_NS::Vertex* v = reinterpret_cast<HEXA_NS::Vertex*>(grHexElt);
1745 if ( _node.count(v)>0 ){
1746 aGrEltIDs.push_back(_node[v]);
1748 if(MYDEBUG) MESSAGE("GROUP OF VERTEX NODES: nodes for vertex (id = "<<v->getId()<<") not found");
1752 default : ASSERT(false);
1756 // B)Filling the group on SMESH
1757 SMESHDS_Group* aGroupDS = dynamic_cast<SMESHDS_Group*>( aGr->GetGroupDS() );
1759 for ( int i=0; i < aGrEltIDs.size(); i++ ) {
1760 aGroupDS->SMDSGroup().Add( aGrEltIDs[i] );
1763 if(MYDEBUG) MESSAGE("_fillGroup() : end >>>>>>>>");