src/HEXABLOCK/HexElements_check.cxx

   1
   2 // C++ : Controle arguments
   3
   4 // Copyright (C) 2009-2023  CEA, EDF
   5 //
   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.
  10 //
  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.
  15 //
  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
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HexElements.hxx"
  24
  25 #include "HexDocument.hxx"
  26 #include "HexVector.hxx"
  27 #include "HexHexa.hxx"
  28 #include "HexQuad.hxx"
  29 #include "HexEdge.hxx"
  30 #include "HexVertex.hxx"
  31
  32 #include "HexGlobale.hxx"
  33
  34 BEGIN_NAMESPACE_HEXA
  35
  36 double calcul_phimax (double radhole, double radext, bool sphere);
  37
  38 // ======================================================== checkSystem
  39 void Elements::checkSystem (int arg, Vector* vx, Vector* vy, Vector* vz,
  40                             double* ux, double* uy, double* uz)
  41 {
  42    int jer [DIM3];
  43    jer[dir_x] = vx->getUnitVector (ux);
  44    jer[dir_y] = vy->getUnitVector (uy);
  45    jer[dir_z] = vz->getUnitVector (uz);
  46
  47    bool sortie = false;
  48    for (int dd=dir_x ; dd<=dir_z ; dd++)
  49        if (jer[dd]!=HOK)
  50           {
  51           setError (HERR);
  52           sortie = true;
  53           Mess << "Vector number " << dd+1 << " is null" ;
  54           }
  55    if (sortie)
  56        return;
  57
  58    double prod = prod_mixte (ux, uy, uz);
  59    if (prod > -Epsil && prod < Epsil)
  60       {
  61       setError (HERR);
  62       Mess << "These vectors are coplanar" ;
  63       }
  64 }
  65 // ======================================================== checkAxis
  66 void Elements::checkAxis (Vector* vz, double* uz)
  67 {
  68    int ier = vz->getUnitVector (uz);
  69    if (ier != HOK)
  70       {
  71       setError (HERR);
  72       Mess << " Second argument vector is null" ;
  73       }
  74 }
  75 // ======================================================== checkDiams
  76 void Elements::checkDiams (int arg, double rint, double rext)
  77 {
  78    if (rint < 0.0)
  79       {
  80       setError (HERR);
  81       Mess << "Argument number " << arg << " :";
  82       Mess << "The internal radius must be positive";
  83       }
  84    if (rext < 0.0)
  85       {
  86       setError (HERR);
  87       Mess << "Argument number " << arg+1 << " :";
  88       Mess << "The external radius must be positive";
  89       }
  90    if (rext <= rint )
  91       {
  92       setError (HERR);
  93       Mess << "The internal radius (arg" << arg << " = " << rint
  94            << ") must be smaller than the external (arg" << arg+1
  95            << " = " << rext << ")";
  96       }
  97 }
  98 // ======================================================== checkSystem
  99 void Elements::checkSystem (int arg, Vector* vx, Vector* vz,
 100                                      double* ux, double* uz)
 101 {
 102    int ier1 = vx->getUnitVector (ux);
 103    int ier2 = vz->getUnitVector (uz);
 104
 105    if (ier1!=HOK)
 106       {
 107       setError (HERR);
 108       Mess << "First vector is null" ;
 109       return;
 110       }
 111    else if (ier2!=HOK)
 112       {
 113       setError (HERR);
 114       Mess << "Second vector is null" ;
 115       return;
 116       }
 117
 118    Real3 pv;
 119    prod_vectoriel (ux, uz, pv);
 120    double prod = calc_norme (pv);
 121    if (prod > -Epsil && prod < Epsil)
 122       {
 123       setError (HERR);
 124       Mess << "These two vectors are colinear" ;
 125       }
 126 }
 127 // ======================================================== checkSize
 128 void Elements::checkSize (int arg, int nx, int ny, int nz, bool rad)
 129 {
 130    if (nx <= 0)
 131       {
 132       Mess << "Argument number " << arg << " :";
 133       Mess << "This dimension must be positive";
 134       setError (HERR);
 135       }
 136    if (nz <= 0)
 137       {
 138       Mess << "Argument number " << arg+2 << " :";
 139       Mess << "This dimension must be positive";
 140       setError (HERR);
 141       }
 142
 143    if (rad)
 144       {
 145       if (ny < 3)
 146          {
 147          Mess << "Argument number " << arg+1 << " :";
 148          Mess << "The number of sectors must be greather than 3";
 149          setError (HERR);
 150          }
 151       }
 152    else if (nx < 1)
 153       {
 154       Mess << "Argument number " << arg+1 << " :";
 155       Mess << "This dimension must be positive";
 156       setError (HERR);
 157       }
 158 }
 159 // ======================================================== checkVector
 160 void Elements::checkVector (int arg, RealVector& table, int lgmin, bool relative)
 161 {
 162    int nbre = table.size();
 163    if (nbre<lgmin)
 164       {
 165       setError (HERR);
 166       Mess << "Argument number " << arg+1 << " :";
 167       Mess << "Size of real vector must be greather or equal than " << lgmin;
 168       return;
 169       }
 170
 171    for (int nv=0 ; nv<nbre ; nv++)
 172        {
 173        double val = table [nv];
 174
 175        if (relative && (val<=0.0 || val >= 1.0))
 176           {
 177           setError (HERR);
 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;
 182           return;
 183           }
 184
 185        if (nv>0 && val <= table [nv-1])
 186           {
 187           setError (HERR);
 188           Mess << "Argument number " << arg+1 << " :";
 189           Mess << "Vector is not growing" ;
 190           Mess <<    "(value[" << nv   << "] = " << val
 191                << " <= value[" << nv-1 << "] = " << table [nv-1]  ;
 192           return;
 193           }
 194        }
 195 }
 196 // ======================================================== checkPipe
 197 void Elements::checkPipe (int arg, double rint, double rext, double angle,
 198                                                              double hauteur)
 199 {
 200    if (rint <= Epsil)
 201       {
 202       setError (HERR);
 203       Mess << "Argument number " << arg << " :";
 204       Mess << "Internal radius (=" << rint << ") must be greather than 0";
 205       }
 206    if (rext <= Epsil)
 207       {
 208       setError (HERR);
 209       Mess << "Argument number " << arg+1 << " :";
 210       Mess << "External radius (=" << rext << ") must be greather than 0";
 211       }
 212    if (rint >= rext)
 213       {
 214       setError (HERR);
 215       Mess << "Arguments number " << arg << " and " << arg+1 << " :";
 216       Mess << "Internal radius (" << rint
 217            << ") must be smaller than external (=" << rext << ")";
 218       }
 219    if (angle <= Epsil)
 220       {
 221       setError (HERR);
 222       Mess << "Argument number " << arg+2 << " :";
 223       Mess << "Angle (=" << angle << ") must be greather than";
 224       }
 225    if (hauteur <= Epsil)
 226       {
 227       setError (HERR);
 228       Mess << "Argument number " << arg+3 << " :";
 229       Mess << "Height (=" << hauteur << ") must be greather than";
 230       }
 231 }
 232 // ======================================================== checkOrig
 233 void Elements::checkOrig (int arg, Vertex* orig)
 234 {
 235    if (orig==NULL || orig->isBad())
 236       {
 237       setError (HERR);
 238       Mess << "Argument number " << arg << " :";
 239       Mess << "Origin vertex is not valid";
 240       }
 241 }
 242 // ======================================================== checkQuads
 243 void Elements::checkQuads (Quads& tquad)
 244 {
 245    int nbre = tquad.size();
 246    if (nbre==0)
 247       {
 248       setError (HERR);
 249       Mess << "Start quads list is empty";
 250       }
 251
 252    for (int nro=0 ; nro<nbre ; nro++)
 253        {
 254        checkQuad (tquad [nro], nro+1);
 255        }
 256 }
 257 // ======================================================== checkQuad
 258 void Elements::checkQuad (Quad* quad, int nro)
 259 {
 260    if (quad==NULL || quad->isBad())
 261       {
 262       setError (HERR);
 263       if (nro>0)
 264          Mess << "Start quad nr " << nro << " is not valid";
 265       else if (nro<0)
 266          Mess << "Target quad must is not valid";
 267       else
 268          Mess << "Start quad is not valid";
 269       return;
 270       }
 271
 272    if (quad->getNbrParents () == 2)
 273       {
 274       setError (HERR);
 275       if (nro>0)
 276          Mess << "Start quad nr " << nro << " must be an external face";
 277       else if (nro<0)
 278          Mess << "Target quad must be an external face";
 279       else
 280          Mess << "Start quad must be an external face";
 281       return;
 282       }
 283 }
 284 // ======================================================== checkSense
 285 void Elements::checkSense (int nro, Vertex* v1, Vertex* v2, Quad* quad)
 286 {
 287    if (v1==NULL || v1->isBad())
 288       {
 289       setError (HERR);
 290       Mess << "Argument nr " << nro << " : vertex is not valid";
 291       return;
 292       }
 293
 294    if (v2==NULL || v2->isBad())
 295       {
 296       setError (HERR);
 297       Mess << "Argument nr " << nro+1 << " : vertex is not valid";
 298       return;
 299       }
 300    Edge* edge = quad -> findEdge (v1, v2);
 301    if (edge!=NULL)
 302        return;
 303
 304    cpchar where = nro < 4 ? "start" : "target";
 305    setError (HERR);
 306    Mess << "Arguments nr " << nro << " and " << nro+2
 307         << " : these vertices do not define an edge of the "
 308         << where << " quad";
 309 }
 310 // ======================================================== checkPhi
 311 int Elements::checkPhi (bool sphere, double* orig, double* norm, double rmax,
 312                          double rhole, Vertex* plan, double& phi0, double& phi1)
 313 {
 314    int    ier  = HOK;
 315    double beta = -M_PI/2;
 316
 317    if (plan!=NULL)
 318       {
 319       Real3 oplan;
 320       plan->getPoint (oplan);
 321
 322       Real3  om   = { oplan[0]-orig[0], oplan[1]-orig[1], oplan[2]-orig[2] };
 323       double hmin = prod_scalaire (norm, om);
 324       if (hmin >= rmax)
 325          {
 326          ier = HERR;
 327          setError (ier);
 328          Mess << "Sphere is under the horizontal plan";
 329          }
 330       else if (hmin > -rmax)
 331          beta = asin (hmin/rmax);
 332       }
 333
 334    phi1 = calcul_phimax (rhole, rmax, sphere);
 335    phi0 = std::max (-phi1, beta);
 336
 337    phi1 = rad2degres (phi1);
 338    phi0 = rad2degres (phi0);
 339    return ier;
 340 }
 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)
 346       AH = AB . ua
 347       (mesure algebrique de AH = produit scalaire de AB par ua
 348
 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)
 353 {
 354    int ier = HOK;
 355    left = right = true;
 356
 357    Real3 ab, cprim;
 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);
 361
 362    for (int dd=dir_x ; dd<=dir_z ; dd++)
 363        {
 364        center [dd] = pa [dd] + ah * ua [dd];
 365        cprim  [dd] = pb [dd] + bh * ub [dd];
 366        }
 367                               // Controles :
 368    ah = fabs (ah);
 369    bh = fabs (bh);
 370
 371    double dist  = calc_distance (center, cprim);
 372    double dcmin = lga*0.05;
 373    double lgmin = (ah - rb)*1.1;
 374
 375    if (ah < ra )
 376       {
 377       ier = HERR;
 378       setError (ier);
 379       Mess << "The base of the big cylinder is included in the smaller" ;
 380       }
 381    else if (lga < lgmin)
 382       {
 383       ier = HERR;
 384       setError (ier);
 385       Mess << "The big cylinder is not long enough too reach the smaller" ;
 386       Mess << "Actual lengh   = " << lga;
 387       Mess << "Minimal length = " << lgmin;
 388       }
 389    else if (lgb < bh - ra)
 390       {
 391       ier = HERR;
 392       setError (ier);
 393       Mess << "The small cylinder is not long enough too reach the large" ;
 394       Mess << "Actual lengh   = " << lgb;
 395       Mess << "Minimal length = " << bh-rb;
 396       }
 397    else if (dist > dcmin)
 398       {
 399       ier = HERR;
 400       setError (ier);
 401       Mess << "Cylinders axes are not secant";
 402       Mess << "Distance between axes = " << dist;
 403       }
 404    else
 405       {
 406       if (bh < rb)
 407          left = false;
 408       if (bh + rb > lgb )
 409          right = false;
 410       }
 411
 412    if (NOT left && NOT right)
 413       {
 414       ier = HERR;
 415       setError (ier);
 416       Mess << "The small cylinder is included in the large" ;
 417       }
 418
 419    return ier;
 420 }
 421 // ======================================================== checkDisco
 422 void Elements::checkDisco (Hexa* cell, Vertex* element)
 423 {
 424    if (BadElement (cell))
 425       {
 426       Mess << "Hexaedron is not valid";
 427       setError (HERR);
 428       return;
 429       }
 430
 431    if (BadElement (element))
 432       {
 433       Mess << "Vertex is not valid";
 434       setError (HERR);
 435       return;
 436       }
 437
 438    int node = cell->findVertex (element);
 439    if (node==NOTHING)
 440       {
 441       setError (HERR);
 442       Mess << "Vertex doesn't belong to Hexaedron";
 443       return;
 444       }
 445 }
 446 // ======================================================== checkDisco
 447 void Elements::checkDisco (Hexa* cell, Edge* element)
 448 {
 449    if (BadElement (cell))
 450       {
 451       Mess << "Hexaedron is not valid";
 452       setError (HERR);
 453       return;
 454       }
 455
 456    if (BadElement (element))
 457       {
 458       Mess << "Edge is not valid";
 459       setError (HERR);
 460       return;
 461       }
 462
 463    int node = cell->findEdge (element);
 464    if (node==NOTHING)
 465       {
 466       setError (HERR);
 467       Mess << "Edge doesn't belong to Hexaedron";
 468       return;
 469       }
 470 }
 471 // ======================================================== checkDisco
 472 void Elements::checkDisco (Hexa* cell, Quad* element)
 473 {
 474    if (BadElement (cell))
 475       {
 476       Mess << "Hexaedron is not valid";
 477       setError (HERR);
 478       return;
 479       }
 480
 481    if (BadElement (element))
 482       {
 483       Mess << "Quad is not valid";
 484       setError (HERR);
 485       return;
 486       }
 487
 488    int node = cell->findQuad (element);
 489    if (node==NOTHING)
 490       {
 491       setError (HERR);
 492       Mess << "Quadrangle doesn't belong to Hexaedron";
 493       return;
 494       }
 495 }
 496 // ======================================================== checkDisco
 497 void Elements::checkDisco (Hexas& thexas, Edges& tedges)
 498 {
 499    int nbedges = tedges.size();
 500    int nbhexas = thexas.size();
 501
 502    if (nbhexas != nbedges)
 503       {
 504       Mess << " **** Error in Document::disconnectEdges";
 505       Mess << " **** Number of Edges and number of Hexas are different";
 506       setError (HERR);
 507       return;
 508       }
 509
 510    for (int nro=0 ; nro<nbedges ; nro++)
 511        {
 512        Edge*  edge = tedges [nro];
 513        Hexa*  hexa = thexas [nro];
 514        if (BadElement (edge))
 515           {
 516           Mess << " **** Eddge number " << nro+1 << " is incorrect";
 517           setError (HERR);
 518           return;
 519           }
 520
 521        if (BadElement (hexa))
 522           {
 523           Mess << " **** Hexa number " << nro+1 << " is incorrect";
 524           setError (HERR);
 525           return;
 526           }
 527
 528        int ned = hexa->findEdge (edge);
 529        if (ned==NOTHING)
 530           {
 531           Mess << " **** Edge number " << nro+1
 532                << " doesnt belong to corresponding hexa";
 533           setError (HERR);
 534           return;
 535           }
 536        if (nro>0)
 537           {
 538           Vertex* vertex = edge->commonVertex (tedges[nro-1]);
 539           if (vertex==NULL)
 540              {
 541              setError (HERR);
 542              Mess << " **** Edge number " << nro
 543                   << " doesnt intesect next edge";
 544              }
 545           }
 546        }
 547 }
 548 // ====================================================== checkContour
 549 void Elements::checkContour (Quads& tquad, Vertex* v1, Vertex* v2, bool target,
 550                              Edges& tedge)
 551 {
 552    tedge.clear ();
 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" ;
 557
 558    Edge*  edge1 = el_root->findEdge (v1, v2);
 559    if (BadElement (edge1))
 560       {
 561       setError (HERR);
 562       Mess << nmedge;
 563       return;
 564       }
 565
 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)
 570         {
 571         Quad* quad = tquad [nq];
 572         if (BadElement (quad))
 573            {
 574            setError (HERR);
 575            Mess << who << " quad nr " << nq+1 << " is not valid";
 576            return;
 577            }
 578          else if (target && quad->getNbrParents() != 1)
 579            {
 580            setError (HERR);
 581            Mess << " Target quad nr " << nq+1
 582                 << " must be an external face";
 583            return;
 584            }
 585
 586         for (int ned=0 ; ned<QUAD4 ; ++ned)
 587             {
 588             Edge* edge = quad->getEdge (ned);
 589             edge_count [edge] += 1;
 590             }
 591         }
 592
 593    int pos1 = edge_count [edge1];
 594    if (pos1==0)
 595        {
 596        setError (HERR);
 597        Mess << nmedge << " of the " << who << " quads";
 598        return;
 599        }
 600    else if (pos1==2)
 601        {
 602        setError (HERR);
 603        Mess << nmedge << " of the " << who << " contour";
 604        return;
 605        }
 606
 607    tedge.push_back (edge1);
 608    Vertex* vlast = v2;
 609    Edge*   elast = edge1;
 610    while (vlast != v1)
 611          {
 612          int   nbre  = vlast->getNbrParents();
 613          Edge* enext = NULL;
 614          for (int ned=0 ; ned<nbre && enext == NULL; ++ned)
 615              {
 616              Edge* edge = vlast->getParent(ned);
 617              if (edge != elast && edge_count [edge]==1)
 618                  enext = edge;
 619              }
 620          if (enext==NULL)
 621             {
 622             setError (HERR);
 623             Mess << who << " as an unclosed contour";
 624             return;
 625             }
 626          tedge.push_back (enext);
 627          vlast = enext->opposedVertex (vlast);
 628          elast = enext;
 629          }
 630 }
 631 // ====================================================== checkContour
 632 void Elements::checkContour (Quads& tquad, Vertex* v1, Vertex* v2, bool target,
 633                              Vertices& tvertex)
 634 {
 635    tvertex.clear ();
 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" ;
 640
 641    Edge*  edge1 = el_root->findEdge (v1, v2);
 642    if (BadElement (edge1))
 643       {
 644       setError (HERR);
 645       Mess << nmedge;
 646       return;
 647       }
 648
 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)
 653         {
 654         Quad* quad = tquad [nq];
 655         if (BadElement (quad))
 656            {
 657            setError (HERR);
 658            Mess << who << " quad nr " << nq+1 << " is not valid";
 659            return;
 660            }
 661          else if (target && quad->getNbrParents() != 1)
 662            {
 663            setError (HERR);
 664            Mess << " Target quad nr " << nq+1
 665                 << " must be an external face";
 666            return;
 667            }
 668
 669         for (int ned=0 ; ned<QUAD4 ; ++ned)
 670             {
 671             Edge* edge = quad->getEdge (ned);
 672             edge_count [edge] += 1;
 673             }
 674         }
 675
 676    int pos1 = edge_count [edge1];
 677    if (pos1==0)
 678        {
 679        setError (HERR);
 680        Mess << nmedge << " of the " << who << " quads";
 681        return;
 682        }
 683    else if (pos1==2)
 684        {
 685        setError (HERR);
 686        Mess << nmedge << " of the " << who << " contour";
 687        return;
 688        }
 689
 690    tvertex.push_back (v1);
 691    Vertex* vlast = v2;
 692    Edge*   elast = edge1;
 693    while (vlast != v1)
 694          {
 695          tvertex.push_back (vlast);
 696          int   nbre  = vlast->getNbrParents();
 697          Edge* enext = NULL;
 698          for (int ned=0 ; ned<nbre && enext == NULL; ++ned)
 699              {
 700              Edge* edge = vlast->getParent(ned);
 701              if (edge != elast && edge_count [edge]==1)
 702                  enext = edge;
 703              }
 704          if (enext==NULL)
 705             {
 706             setError (HERR);
 707             Mess << who << " as an unclosed contour";
 708             return;
 709             }
 710          vlast = enext->opposedVertex (vlast);
 711          elast = enext;
 712          }
 713 }
 714 // ======================================================== calcul_phimax
 715 double calcul_phimax (double radhole, double radext, bool sphere)
 716 {
 717    double phi = asin (radhole/radext);
 718    if (sphere)
 719        phi = 2*phi;
 720    phi = M_PI*DEMI - phi;
 721    return phi;
 722 }
 723
 724 END_NAMESPACE_HEXA