2 // C++ : Controle arguments
4 // Copyright (C) 2009-2023 CEA, EDF
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, or (at your option) any later version.
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
23 #include "HexElements.hxx"
25 #include "HexDocument.hxx"
26 #include "HexVector.hxx"
27 #include "HexHexa.hxx"
28 #include "HexQuad.hxx"
29 #include "HexEdge.hxx"
30 #include "HexVertex.hxx"
32 #include "HexGlobale.hxx"
36 double calcul_phimax (double radhole, double radext, bool sphere);
38 // ======================================================== checkSystem
39 void Elements::checkSystem (int arg, Vector* vx, Vector* vy, Vector* vz,
40 double* ux, double* uy, double* uz)
43 jer[dir_x] = vx->getUnitVector (ux);
44 jer[dir_y] = vy->getUnitVector (uy);
45 jer[dir_z] = vz->getUnitVector (uz);
48 for (int dd=dir_x ; dd<=dir_z ; dd++)
53 Mess << "Vector number " << dd+1 << " is null" ;
58 double prod = prod_mixte (ux, uy, uz);
59 if (prod > -Epsil && prod < Epsil)
62 Mess << "These vectors are coplanar" ;
65 // ======================================================== checkAxis
66 void Elements::checkAxis (Vector* vz, double* uz)
68 int ier = vz->getUnitVector (uz);
72 Mess << " Second argument vector is null" ;
75 // ======================================================== checkDiams
76 void Elements::checkDiams (int arg, double rint, double rext)
81 Mess << "Argument number " << arg << " :";
82 Mess << "The internal radius must be positive";
87 Mess << "Argument number " << arg+1 << " :";
88 Mess << "The external radius must be positive";
93 Mess << "The internal radius (arg" << arg << " = " << rint
94 << ") must be smaller than the external (arg" << arg+1
95 << " = " << rext << ")";
98 // ======================================================== checkSystem
99 void Elements::checkSystem (int arg, Vector* vx, Vector* vz,
100 double* ux, double* uz)
102 int ier1 = vx->getUnitVector (ux);
103 int ier2 = vz->getUnitVector (uz);
108 Mess << "First vector is null" ;
114 Mess << "Second vector is null" ;
119 prod_vectoriel (ux, uz, pv);
120 double prod = calc_norme (pv);
121 if (prod > -Epsil && prod < Epsil)
124 Mess << "These two vectors are colinear" ;
127 // ======================================================== checkSize
128 void Elements::checkSize (int arg, int nx, int ny, int nz, bool rad)
132 Mess << "Argument number " << arg << " :";
133 Mess << "This dimension must be positive";
138 Mess << "Argument number " << arg+2 << " :";
139 Mess << "This dimension must be positive";
147 Mess << "Argument number " << arg+1 << " :";
148 Mess << "The number of sectors must be greather than 3";
154 Mess << "Argument number " << arg+1 << " :";
155 Mess << "This dimension must be positive";
159 // ======================================================== checkVector
160 void Elements::checkVector (int arg, RealVector& table, int lgmin, bool relative)
162 int nbre = table.size();
166 Mess << "Argument number " << arg+1 << " :";
167 Mess << "Size of real vector must be greather or equal than " << lgmin;
171 for (int nv=0 ; nv<nbre ; nv++)
173 double val = table [nv];
175 if (relative && (val<=0.0 || val >= 1.0))
178 Mess << "Argument number " << arg+1 << " :";
179 Mess << "This vector contains relative lengths" ;
180 Mess << "Values must be between 0 and 1" ;
181 Mess << "(value[" << nv << "] = " << val;
185 if (nv>0 && val <= table [nv-1])
188 Mess << "Argument number " << arg+1 << " :";
189 Mess << "Vector is not growing" ;
190 Mess << "(value[" << nv << "] = " << val
191 << " <= value[" << nv-1 << "] = " << table [nv-1] ;
196 // ======================================================== checkPipe
197 void Elements::checkPipe (int arg, double rint, double rext, double angle,
203 Mess << "Argument number " << arg << " :";
204 Mess << "Internal radius (=" << rint << ") must be greather than 0";
209 Mess << "Argument number " << arg+1 << " :";
210 Mess << "External radius (=" << rext << ") must be greather than 0";
215 Mess << "Arguments number " << arg << " and " << arg+1 << " :";
216 Mess << "Internal radius (" << rint
217 << ") must be smaller than external (=" << rext << ")";
222 Mess << "Argument number " << arg+2 << " :";
223 Mess << "Angle (=" << angle << ") must be greather than";
225 if (hauteur <= Epsil)
228 Mess << "Argument number " << arg+3 << " :";
229 Mess << "Height (=" << hauteur << ") must be greather than";
232 // ======================================================== checkOrig
233 void Elements::checkOrig (int arg, Vertex* orig)
235 if (orig==NULL || orig->isBad())
238 Mess << "Argument number " << arg << " :";
239 Mess << "Origin vertex is not valid";
242 // ======================================================== checkQuads
243 void Elements::checkQuads (Quads& tquad)
245 int nbre = tquad.size();
249 Mess << "Start quads list is empty";
252 for (int nro=0 ; nro<nbre ; nro++)
254 checkQuad (tquad [nro], nro+1);
257 // ======================================================== checkQuad
258 void Elements::checkQuad (Quad* quad, int nro)
260 if (quad==NULL || quad->isBad())
264 Mess << "Start quad nr " << nro << " is not valid";
266 Mess << "Target quad must is not valid";
268 Mess << "Start quad is not valid";
272 if (quad->getNbrParents () == 2)
276 Mess << "Start quad nr " << nro << " must be an external face";
278 Mess << "Target quad must be an external face";
280 Mess << "Start quad must be an external face";
284 // ======================================================== checkSense
285 void Elements::checkSense (int nro, Vertex* v1, Vertex* v2, Quad* quad)
287 if (v1==NULL || v1->isBad())
290 Mess << "Argument nr " << nro << " : vertex is not valid";
294 if (v2==NULL || v2->isBad())
297 Mess << "Argument nr " << nro+1 << " : vertex is not valid";
300 Edge* edge = quad -> findEdge (v1, v2);
304 cpchar where = nro < 4 ? "start" : "target";
306 Mess << "Arguments nr " << nro << " and " << nro+2
307 << " : these vertices do not define an edge of the "
310 // ======================================================== checkPhi
311 int Elements::checkPhi (bool sphere, double* orig, double* norm, double rmax,
312 double rhole, Vertex* plan, double& phi0, double& phi1)
315 double beta = -M_PI/2;
320 plan->getPoint (oplan);
322 Real3 om = { oplan[0]-orig[0], oplan[1]-orig[1], oplan[2]-orig[2] };
323 double hmin = prod_scalaire (norm, om);
328 Mess << "Sphere is under the horizontal plan";
330 else if (hmin > -rmax)
331 beta = asin (hmin/rmax);
334 phi1 = calcul_phimax (rhole, rmax, sphere);
335 phi0 = std::max (-phi1, beta);
337 phi1 = rad2degres (phi1);
338 phi0 = rad2degres (phi0);
341 // ======================================================== checkInter
342 /* === Intersection de 2 cylindres (A,ua) et (B,ub)
343 ua et ub : vecteurs directeurs unitaires.
344 A a le gros diametre, B le petit
345 H = proj ortho de B sur (A, ua)
347 (mesure algebrique de AH = produit scalaire de AB par ua
349 ---------------------------------------------------------------- */
350 int Elements::checkInter (double* pa, double* ua, double ra, double lga,
351 double* pb, double* ub, double rb, double lgb,
352 double* center, bool& left, bool& right)
358 calc_vecteur (pa, pb, ab);
359 double ah = prod_scalaire (ab, ua); // ah > 0 ???? : A voir avec FK
360 double bh = -prod_scalaire (ab, ub);
362 for (int dd=dir_x ; dd<=dir_z ; dd++)
364 center [dd] = pa [dd] + ah * ua [dd];
365 cprim [dd] = pb [dd] + bh * ub [dd];
371 double dist = calc_distance (center, cprim);
372 double dcmin = lga*0.05;
373 double lgmin = (ah - rb)*1.1;
379 Mess << "The base of the big cylinder is included in the smaller" ;
381 else if (lga < lgmin)
385 Mess << "The big cylinder is not long enough too reach the smaller" ;
386 Mess << "Actual lengh = " << lga;
387 Mess << "Minimal length = " << lgmin;
389 else if (lgb < bh - ra)
393 Mess << "The small cylinder is not long enough too reach the large" ;
394 Mess << "Actual lengh = " << lgb;
395 Mess << "Minimal length = " << bh-rb;
397 else if (dist > dcmin)
401 Mess << "Cylinders axes are not secant";
402 Mess << "Distance between axes = " << dist;
412 if (NOT left && NOT right)
416 Mess << "The small cylinder is included in the large" ;
421 // ======================================================== checkDisco
422 void Elements::checkDisco (Hexa* cell, Vertex* element)
424 if (BadElement (cell))
426 Mess << "Hexaedron is not valid";
431 if (BadElement (element))
433 Mess << "Vertex is not valid";
438 int node = cell->findVertex (element);
442 Mess << "Vertex doesn't belong to Hexaedron";
446 // ======================================================== checkDisco
447 void Elements::checkDisco (Hexa* cell, Edge* element)
449 if (BadElement (cell))
451 Mess << "Hexaedron is not valid";
456 if (BadElement (element))
458 Mess << "Edge is not valid";
463 int node = cell->findEdge (element);
467 Mess << "Edge doesn't belong to Hexaedron";
471 // ======================================================== checkDisco
472 void Elements::checkDisco (Hexa* cell, Quad* element)
474 if (BadElement (cell))
476 Mess << "Hexaedron is not valid";
481 if (BadElement (element))
483 Mess << "Quad is not valid";
488 int node = cell->findQuad (element);
492 Mess << "Quadrangle doesn't belong to Hexaedron";
496 // ======================================================== checkDisco
497 void Elements::checkDisco (Hexas& thexas, Edges& tedges)
499 int nbedges = tedges.size();
500 int nbhexas = thexas.size();
502 if (nbhexas != nbedges)
504 Mess << " **** Error in Document::disconnectEdges";
505 Mess << " **** Number of Edges and number of Hexas are different";
510 for (int nro=0 ; nro<nbedges ; nro++)
512 Edge* edge = tedges [nro];
513 Hexa* hexa = thexas [nro];
514 if (BadElement (edge))
516 Mess << " **** Eddge number " << nro+1 << " is incorrect";
521 if (BadElement (hexa))
523 Mess << " **** Hexa number " << nro+1 << " is incorrect";
528 int ned = hexa->findEdge (edge);
531 Mess << " **** Edge number " << nro+1
532 << " doesnt belong to corresponding hexa";
538 Vertex* vertex = edge->commonVertex (tedges[nro-1]);
542 Mess << " **** Edge number " << nro
543 << " doesnt intesect next edge";
548 // ====================================================== checkContour
549 void Elements::checkContour (Quads& tquad, Vertex* v1, Vertex* v2, bool target,
553 cpchar who = target ? "Target" : "Pattern";
554 std::string nmedge = target ? "Vertices of target (args 5 and 6)"
555 : "Vertices of pattern (args 3 and 4)" ;
556 nmedge += "don't define an edge" ;
558 Edge* edge1 = el_root->findEdge (v1, v2);
559 if (BadElement (edge1))
566 std::map <Edge*, int> edge_count;
567 std::map <Edge*, int> :: iterator iter;
568 int nbre = tquad.size();
569 for (int nq=0 ; nq<nbre ; ++nq)
571 Quad* quad = tquad [nq];
572 if (BadElement (quad))
575 Mess << who << " quad nr " << nq+1 << " is not valid";
578 else if (target && quad->getNbrParents() != 1)
581 Mess << " Target quad nr " << nq+1
582 << " must be an external face";
586 for (int ned=0 ; ned<QUAD4 ; ++ned)
588 Edge* edge = quad->getEdge (ned);
589 edge_count [edge] += 1;
593 int pos1 = edge_count [edge1];
597 Mess << nmedge << " of the " << who << " quads";
603 Mess << nmedge << " of the " << who << " contour";
607 tedge.push_back (edge1);
612 int nbre = vlast->getNbrParents();
614 for (int ned=0 ; ned<nbre && enext == NULL; ++ned)
616 Edge* edge = vlast->getParent(ned);
617 if (edge != elast && edge_count [edge]==1)
623 Mess << who << " as an unclosed contour";
626 tedge.push_back (enext);
627 vlast = enext->opposedVertex (vlast);
631 // ====================================================== checkContour
632 void Elements::checkContour (Quads& tquad, Vertex* v1, Vertex* v2, bool target,
636 cpchar who = target ? "Target" : "Pattern";
637 std::string nmedge = target ? "Vertices of target (args 4 and 6)"
638 : "Vertices of pattern (args 3 and 5)" ;
639 nmedge += "don't define an edge" ;
641 Edge* edge1 = el_root->findEdge (v1, v2);
642 if (BadElement (edge1))
649 std::map <Edge*, int> edge_count;
650 std::map <Edge*, int> :: iterator iter;
651 int nbre = tquad.size();
652 for (int nq=0 ; nq<nbre ; ++nq)
654 Quad* quad = tquad [nq];
655 if (BadElement (quad))
658 Mess << who << " quad nr " << nq+1 << " is not valid";
661 else if (target && quad->getNbrParents() != 1)
664 Mess << " Target quad nr " << nq+1
665 << " must be an external face";
669 for (int ned=0 ; ned<QUAD4 ; ++ned)
671 Edge* edge = quad->getEdge (ned);
672 edge_count [edge] += 1;
676 int pos1 = edge_count [edge1];
680 Mess << nmedge << " of the " << who << " quads";
686 Mess << nmedge << " of the " << who << " contour";
690 tvertex.push_back (v1);
695 tvertex.push_back (vlast);
696 int nbre = vlast->getNbrParents();
698 for (int ned=0 ; ned<nbre && enext == NULL; ++ned)
700 Edge* edge = vlast->getParent(ned);
701 if (edge != elast && edge_count [edge]==1)
707 Mess << who << " as an unclosed contour";
710 vlast = enext->opposedVertex (vlast);
714 // ======================================================== calcul_phimax
715 double calcul_phimax (double radhole, double radext, bool sphere)
717 double phi = asin (radhole/radext);
720 phi = M_PI*DEMI - phi;