2 // C++ : Gestion des Quadrangles
4 // Copyright (C) 2009-2013 CEA/DEN, EDF R&D
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 #include "HexQuad.hxx"
24 #include "HexDocument.hxx"
25 #include "HexHexa.hxx"
26 #include "HexElements.hxx"
27 #include "HexGlobale.hxx"
29 #include "HexXmlWriter.hxx"
30 #include "HexNewShape.hxx"
31 #include "HexFaceShape.hxx"
37 // ======================================================== Constructeur
38 Quad::Quad (Vertex* va, Vertex* vb, Vertex* vc, Vertex* vd)
39 : EltBase (va->dad(), EL_QUAD)
46 q_orientation = Q_UNDEFINED;
48 for (int nro=0 ; nro<QUAD4 ; nro++)
50 q_edge [nro] = new Edge (q_vertex[nro],
51 q_vertex[(nro+1) MODULO QUAD4]);
53 if (BadElement (q_vertex [nro]) || BadElement (q_edge [nro]))
56 for (int nv=nro+1 ; nv<QUAD4 ; nv++)
57 if (q_vertex[nv] == q_vertex[nro])
61 if (el_root != NULL && el_status==HOK)
62 el_root->addQuad (this);
65 // ======================================================== Constructeur bis
66 Quad::Quad (Edge* ea, Edge* eb, Edge* ec, Edge* ed)
67 : EltBase (ea->dad(), EL_QUAD)
74 q_orientation = Q_UNDEFINED;
76 for (int nro=0 ; nro<QUAD4 ; nro++)
78 if (BadElement (q_edge[nro]))
82 for (int nv=nro+1 ; nv<QUAD4 ; nv++)
83 if (q_edge[nv] == q_edge[nro])
90 // Cond necessaire : ea disjoint de ec (opposes)
91 int nc = ea->inter (ec);
114 for (int nro=0 ; nro<QUAD4 ; nro++)
116 int prec = (nro+1) MODULO QUAD4;
118 int nc = q_edge[nro] -> inter (q_edge[prec]);
120 node = q_edge[nro]->getVertex (nc);
123 q_vertex [prec] = node;
128 printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
129 printf (" +++ Quadrangle impossible \n");
130 printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
132 printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
134 for (int ned=0; ned<QUAD4; ned++)
136 q_edge[ned]->dumpPlus ();
138 HexDump (q_vertex[0]);
139 HexDump (q_vertex[1]);
140 HexDump (q_vertex[2]);
141 HexDump (q_vertex[3]);
143 printf (" +++++++++++++++++++++++++++++++++++++++++++ \n");
144 fatal_error ("Quadrangle impossible");
147 if (el_root != NULL && el_status==HOK)
148 el_root->addQuad (this);
151 // ======================================================== Constructeur ter
152 Quad::Quad (Quad* other)
153 : EltBase (other->dad(), EL_QUAD)
155 for (int nro=0 ; nro<QUAD4 ; nro++)
158 q_vertex [nro] = NULL;
160 q_orientation = Q_UNDEFINED;
163 if (el_root != NULL && el_status==HOK)
164 el_root->addQuad (this);
166 // ============================================================ getEdge
167 Edge* Quad::getEdge (int nro)
170 if (nro >=0 && nro < QUAD4 && el_status == HOK && q_edge [nro]->isValid())
173 DumpStart ("getEdge", nro);
177 // ============================================================ getVertex
178 Vertex* Quad::getVertex (int nro)
181 if (nro >=0 && nro < QUAD4 && el_status == HOK && q_vertex [nro]->isValid())
182 elt = q_vertex [nro];
184 DumpStart ("getVertex", nro);
188 // ========================================================= majReferences
189 void Quad::majReferences ()
191 for (int nro=0 ; nro<QUAD4 ; nro++)
192 q_edge [nro] -> addParent (this);
194 // ========================================================= getParent
195 Hexa* Quad::getParent (int nro)
197 return static_cast <Hexa*> (getFather (nro));
199 // ======================================================== anaMerge
200 int Quad::anaMerge (Vertex* v1, Vertex* v2, Vertex* tv1[], Edge* te1[])
203 for (int nro=0 ; orig == NOTHING && nro < QUAD4 ; nro++)
204 if (q_vertex [nro] == v1)
210 int nsp1 = (orig+1) MODULO QUAD4;
211 int nsm1 = (orig+QUAD4-1) MODULO QUAD4;
213 if (q_vertex [nsp1] == v2)
215 for (int nro=0 ; nro < QUAD4 ; nro++)
217 tv1 [nro] = q_vertex [(orig+nro) MODULO QUAD4];
218 te1 [nro] = q_edge [(orig+nro) MODULO QUAD4];
221 else if (q_vertex [nsm1] == v2)
223 for (int nro=0 ; nro < QUAD4 ; nro++)
225 tv1 [nro] = q_vertex [(orig+QUAD4-nro) MODULO QUAD4];
226 te1 [nro] = q_edge [(orig+QUAD4-nro) MODULO QUAD4];
234 // ======================================================== ordoVertex
235 int Quad::ordoVertex (Vertex* v1, Vertex* v2, Vertex* tver[])
238 for (int nro=0 ; orig == NOTHING && nro < QUAD4 ; nro++)
239 if (q_vertex [nro] == v1)
245 int nsp1 = (orig+1) MODULO QUAD4;
246 int nsm1 = (orig+QUAD4-1) MODULO QUAD4;
248 if (q_vertex [nsp1] == v2)
250 for (int nro=0 ; nro < QUAD4 ; nro++)
251 tver [nro] = q_vertex [(orig+nro) MODULO QUAD4];
253 else if (q_vertex [nsm1] == v2)
255 for (int nro=0 ; nro < QUAD4 ; nro++)
256 tver [nro] = q_vertex [(orig+QUAD4-nro) MODULO QUAD4];
263 // ======================================================== ordonner
264 int Quad::ordonner (Vertex* v1, Vertex* v2, Vertex* tver[], Edge* ted[])
266 tver [0] = tver [1] = tver [2] = tver [3] = NULL;
267 ted [0] = ted [1] = ted [2] = ted [3] = NULL;
269 int ier = ordoVertex (v1, v2, tver);
273 for (int nro=0 ; nro < QUAD4 ; nro++)
274 ted [nro] = findEdge (tver[nro], tver [(nro+1) MODULO QUAD4]);
278 // ======================================================== getBrother
279 Quad* Quad::getBrother (StrOrient* orient)
281 /* *****************************
282 printf (" getBrother ");
284 printf (" .. Base : ");
285 orient->v21->printName();
286 orient->v22->printName();
287 printf ("dir=%d, arete=", orient->dir);
288 ***************************** */
290 int n21 = indexVertex (orient->v21);
291 int n22 = indexVertex (orient->v22);
293 int sens = n22 - n21;
294 if (sens > 1) sens -= QUAD4;
295 if (sens < -1) sens += QUAD4;
296 if (sens*sens !=1) return NULL;
300 case OR_LEFT : n22 = n21 - sens;
302 case OR_RIGHT : n21 = n22 + sens;
304 case OR_FRONT : n21 += 2;
310 n21 = (n21 + QUAD4) MODULO QUAD4;
311 n22 = (n22 + QUAD4) MODULO QUAD4;
313 orient->v21 = q_vertex [n21];
314 orient->v22 = q_vertex [n22];
316 Edge* arete = findEdge (orient->v21, orient->v22);
317 // arete->printName("\n");
319 int nbfreres = arete->getNbrParents ();
321 for (int nq = 0 ; nq < nbfreres ; nq++)
323 Quad* next = arete->getParent (nq);
324 if (next!=NULL && next != this )
326 int nbp = next->getNbrParents();
327 Hexa* dad = next->getParent(0);
328 int mark = next->getMark();
329 int mark2 = dad ? dad->getMark() : IS_NONE;
331 if (nbp == 1 && mark2 != IS_MARRIED && mark == IS_NONE)
333 // if (nbp <= 1 && mark == IS_NONE)
339 // ======================================================== coupler
340 int Quad::coupler (Quad* other, StrOrient* orient, Elements* table)
346 cout << " Quads::coupler " << el_name << " -> " << other->getName ()
350 Hexa* hexa = other->getParent(0);
352 setMark (IS_MARRIED);
353 other->setMark (IS_MARRIED);
355 hexa->setMark (IS_MARRIED);
357 for (int ned = 0 ; ned < QUAD4 ; ned++)
359 Edge* arete = q_edge[ned];
360 int nbfreres = arete ->getNbrParents ();
361 for (int nq = 0 ; nq < nbfreres ; nq++)
363 Quad* next = arete->getParent (nq);
364 if (next!=NULL && next != this && next->getMark() > 0)
366 StrOrient new_ori (orient);
367 new_ori.dir = OR_FRONT;
368 Vertex* va = arete->getVertex (V_AMONT);
369 Vertex* vb = arete->getVertex (V_AVAL);
371 // On voit si un point de repere est conserve
372 if (va == orient->v11)
375 new_ori.dir += OR_LEFT;
377 else if (vb == orient->v11)
380 new_ori.dir += OR_LEFT;
383 if (va == orient->v12)
386 new_ori.dir += OR_RIGHT;
388 else if (vb == orient->v12)
391 new_ori.dir += OR_RIGHT;
394 if (new_ori.dir == OR_FRONT)
396 if (definedBy (va, orient->v11))
408 int nro = next->getMark ();
409 Quad* beauf = other->getBrother (&new_ori);
410 int ier = table->coupler (nro, beauf, &new_ori);
413 ier = next->coupler (beauf, &new_ori, table);
421 // ======================================================== getOpposVertex
422 Vertex* Quad::getOpposVertex (Vertex* start)
424 int na = indexVertex (start);
427 return q_vertex [(na+2) MODULO QUAD4];
429 // ======================================================== getOpposEdge
430 Edge* Quad::getOpposEdge (Edge* start, int& sens)
433 int na = indexVertex (start->getVertex (V_AMONT));
434 int nb = indexVertex (start->getVertex (V_AVAL));
436 Vertex* vaprim = q_vertex [(nb+2) MODULO QUAD4];
437 Vertex* vbprim = q_vertex [(na+2) MODULO QUAD4];
439 for (int ned = 0 ; ned < QUAD4 ; ned++)
441 if ( q_edge[ned]->getVertex(V_AMONT) == vaprim
442 && q_edge[ned]->getVertex(V_AVAL ) == vbprim)
447 else if ( q_edge[ned]->getVertex(V_AMONT) == vbprim
448 && q_edge[ned]->getVertex(V_AVAL ) == vaprim)
454 // TODO : traiter l'erreur
455 cout << " ... Probleme dans Quad::getOpposedEdge :" << endl;
456 HexDisplay (el_name);
464 for (int ned = 0 ; ned < QUAD4 ; ned++)
469 // ========================================================= saveXml
470 void Quad::saveXml (XmlWriter* xml)
475 for (int nro=0 ; nro<QUAD4 ; nro++)
477 if (nro>0) edges += " ";
478 edges += q_edge[nro]->getName(buffer);
481 xml->openMark ("Quad");
482 xml->addAttribute ("id", getName (buffer));
483 xml->addAttribute ("edges", edges);
485 xml->addAttribute ("name", el_name);
488 int nbass = tab_assoc.size();
489 for (int nro=0 ; nro<nbass ; nro++)
490 if (tab_assoc[nro] != NULL)
491 tab_assoc[nro]->saveXml (xml);
493 // ======================================================== replaceEdge
494 void Quad::replaceEdge (Edge* old, Edge* par)
496 for (int nro=0 ; nro<QUAD4 ; nro++)
498 if (q_edge[nro]==old)
505 printf (" [%d], ", nro);
506 old->printName (" est remplace par ");
507 par->printName ("\n");
512 // ======================================================== replaceVertex
513 void Quad::replaceVertex (Vertex* old, Vertex* par)
515 for (int nro=0 ; nro<QUAD4 ; nro++)
517 if (q_vertex [nro]==old)
519 q_vertex [nro] = par;
524 printf (" [%d], ", nro);
525 old->printName (" est remplace par ");
526 par->printName ("\n");
531 // ======================================================== dump
537 printf ("*** deleted ***)\n");
541 for (int nro=0 ; nro<QUAD4 ; nro++)
542 PrintName (q_edge[nro]);
546 for (int nro=0 ; nro<QUAD4 ; nro++)
547 PrintName (q_vertex[nro]);
553 printf (" -> (%g, %g, %g)\n", cg[dir_x], cg[dir_y], cg[dir_z]);
555 // ======================================================== dumpPlus
556 void Quad::dumpPlus ()
562 for (int nro=0 ; nro < QUAD4 ; nro++)
564 Vertex* pv = q_vertex[nro];
569 printf ( " (%g, %g, %g)\n", pv->getX(), pv->getY(), pv->getZ());
577 // ======================================================== getOpposEdge (2)
578 Edge* Quad::getOpposEdge (Edge* start)
580 int na = indexEdge (start);
583 return q_edge [(na+2) MODULO QUAD4];
585 // ======================================================== getPerpendicular
586 Edge* Quad::getPerpendicular (Edge* arete, Vertex* node)
588 int na = indexEdge (arete);
592 int nv = arete->index (node);
595 Edge* perp = q_edge [(na+1) MODULO QUAD4];
597 nv = perp->index (node);
601 perp = q_edge [(na+3) MODULO QUAD4];
602 nv = perp->index (node);
608 static cpchar t_ori[] = {"Q_INSIDE", "Q_DIRECT", "Q_INVERSE", "Q_UNDEF"};
609 // ======================================================== setOrientation
610 void Quad::setOrientation (int ori)
613 if (db && (ori==Q_DIRECT || ori==Q_INVERSE))
614 printf (" %s = %s\n", el_name.c_str(), t_ori [ q_orientation ]);
616 // ======================================================== setOrientation
617 int Quad::setOrientation ()
619 q_orientation = Q_INSIDE;
620 if (getNbrParents() != 1)
621 return q_orientation;
623 Real3 cg, orig, pi, pj, vi, vj, vk;
625 Hexa* hexa = getParent(0);
626 hexa->getCenter (cg);
628 /********************************************************************
631 for (int np=0 ; np < QUAD4 ; np++)
633 q_vertex [np ] -> getPoint (orig);
634 q_vertex [(np+1) % 4] -> getPoint (pi);
635 q_vertex [(np+3) % 4] -> getPoint (pj);
637 calc_vecteur (orig, pi, vi);
638 calc_vecteur (orig, pj, vj);
639 calc_vecteur (orig, cg, vk);
640 double pmixte = prod_mixte (vi, vj, vk);
641 q_orientation = pmixte > ZEROR ? Q_DIRECT : Q_INVERSE;
642 if (pmixte>0) printf (">");
648 ******************************************************************* */
649 q_vertex [0] -> getPoint (orig);
650 q_vertex [1] -> getPoint (pi);
651 q_vertex [3] -> getPoint (pj);
653 calc_vecteur (orig, pi, vi);
654 calc_vecteur (orig, pj, vj);
655 calc_vecteur (cg, orig, vk);
657 double pmixte = prod_mixte (vi, vj, vk);
658 q_orientation = pmixte > ZEROR ? Q_DIRECT : Q_INVERSE;
660 printf (" %s := %s\n", el_name.c_str(), t_ori [ q_orientation ]);
661 return q_orientation;
663 // ========================================================== clearAssociation
664 void Quad::clearAssociation ()
667 is_associated = false;
669 // ========================================================== addAssociation
670 int Quad::addAssociation (NewShape* geom, int subid)
675 FaceShape* face = geom->findFace (subid);
676 int ier = addAssociation (face);
680 // ========================================================== addAssociation
681 int Quad::addAssociation (FaceShape* face)
686 face->addAssociation (this);
687 tab_assoc.push_back (face);
688 is_associated = true;
691 // ========================================================== getAssociation
692 FaceShape* Quad::getAssociation (int nro)
694 if (nro < 0 || nro >= (int)tab_assoc.size())
697 return tab_assoc [nro];
699 // ======================================================== commonEdge
700 Edge* Quad::commonEdge (Quad* other)
702 for (int ne1=0 ; ne1<QUAD4 ; ne1++)
703 for (int ne2=0 ; ne2<QUAD4 ; ne2++)
704 if (q_edge [ne1] == other->q_edge [ne2])
709 // ======================================================== Inter
710 int Quad::inter (Quad* other, int& nother)
712 for (int ne1=0 ; ne1<QUAD4 ; ne1++)
713 for (int ne2=0 ; ne2<QUAD4 ; ne2++)
714 if (q_edge [ne1] == other->q_edge [ne2])
723 // ======================================================== Inter (2)
724 Edge* Quad::inter (Quad* other)
726 for (int ne1=0 ; ne1<QUAD4 ; ne1++)
727 for (int ne2=0 ; ne2<QUAD4 ; ne2++)
728 if (q_edge [ne1] == other->q_edge [ne2])
732 // ============================================================ definedBy (v)
733 bool Quad::definedBy (Vertex* v1, Vertex* v2)
735 for (int n1=0 ; n1< QUAD4 ; n1++)
736 if (v1 == q_vertex[n1] && v2 == q_vertex[(n1+2) MODULO QUAD4])
741 // ============================================================ definedBy (e)
742 bool Quad::definedBy (Edge* e1, Edge* e2)
744 if (e1==e2 || BadElement (e1) || BadElement (e2))
747 bool f1=false, f2=false;
748 for (int ned=0 ; ned< QUAD4 ; ned++)
749 if (e1 == q_edge[ned]) f1 = true;
750 else if (e2 == q_edge[ned]) f2 = true;
751 // if (e1 == q_edge[ned] && e2 == q_edge[(ned+2) MODULO QUAD4]) return true;
755 // =============================================================== findEdge
756 Edge* Quad::findEdge (Vertex* v1, Vertex* v2)
758 for (int nro=0 ; nro< QUAD4 ; nro++)
760 Vertex* va = q_edge[nro]->getVertex(V_AMONT) ;
761 Vertex* vb = q_edge[nro]->getVertex(V_AVAL) ;
762 if ((v1==va && v2==vb) || (v1==vb && v2==va))
768 // =============================================================== indexVertex
769 int Quad::indexVertex (Vertex* elt)
771 for (int n1=0 ; n1< QUAD4 ; n1++)
772 if (elt == q_vertex[n1])
777 // =============================================================== indexEdge
778 int Quad::indexEdge (Edge* elt)
780 for (int n1=0 ; n1< QUAD4 ; n1++)
781 if (elt == q_edge[n1])
786 // =============================================================== setColor
787 void Quad::setColor (double val)
789 for (int n1=0 ; n1< QUAD4 ; n1++)
790 q_vertex[n1] -> setColor (val);
792 // =============================================================== duplicate
793 void Quad::duplicate ()
795 q_orientation = Q_UNDEFINED;
796 q_clone = new Quad (GetClone (q_edge [E_A]),
797 GetClone (q_edge [E_B]),
798 GetClone (q_edge [E_C]),
799 GetClone (q_edge [E_D]));
800 q_clone->tab_assoc = tab_assoc;
802 // ============================================================ nearestVertex
803 Vertex* Quad::nearestVertex (Vertex* other)
805 if (BadElement (other))
810 int nbre = countVertex ();
811 for (int nro=0 ; nro<nbre ; nro++)
813 Vertex* vert = getVertex (nro);
814 double dist = other->dist2 (vert);
823 // =============================================================== reorienter
824 void Quad::reorienter ()
826 if (q_orientation != Q_INVERSE)
829 Edge* edge = q_edge [E_B];
830 q_edge [E_B] = q_edge [E_D];
833 q_vertex [E_A] = q_edge [E_D]->commonVertex(q_edge [E_A]);
834 q_vertex [E_B] = q_edge [E_A]->commonVertex(q_edge [E_B]);
835 q_vertex [E_C] = q_edge [E_B]->commonVertex(q_edge [E_C]);
836 q_vertex [E_D] = q_edge [E_C]->commonVertex(q_edge [E_D]);
839 printf (" %s est reoriente\n", el_name.c_str());
840 q_orientation = Q_DIRECT;
842 // =============================================================== getCenter
843 double* Quad::getCenter (double* center)
845 center[dir_x] = center[dir_y] = center[dir_z] = 0;
846 if (BadElement (this))
849 for (int nv=0 ; nv<QUAD4 ; nv++)
851 if (BadElement (q_vertex [nv]))
853 center [dir_x] += q_vertex[nv]->getX()/4;
854 center [dir_y] += q_vertex[nv]->getY()/4;
855 center [dir_z] += q_vertex[nv]->getZ()/4;
859 // =============================================================== dist2
860 double Quad::dist2 (double* point)
864 double d2 = carre (point[dir_x] - center[dir_x])
865 + carre (point[dir_y] - center[dir_y])
866 + carre (point[dir_z] - center[dir_z]) ;
869 // =============================================================== opposedHexa
870 Hexa* Quad::opposedHexa (Hexa* hexa)
872 int nbre = getNbrParents ();
873 for (int nro=0 ; nro <nbre ; ++nro)
875 Hexa* dad = getParent (nro);
876 if (dad!= NULL && dad->isValid() && dad != hexa)