src/Tools/MeshCut/MeshCut_Fonctions.cxx

   1 // Copyright (C) 2006-2019  EDF R&D
   2 //
   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, or (at your option) any later version.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "MeshCut_Fonctions.hxx"
  21 #include "MeshCut_Globals.hxx"
  22
  23 #include <iostream>
  24 #include <cmath>
  25
  26 using namespace MESHCUT;
  27 using namespace std;
  28
  29 //=================================================================================================
  30 //                                    intersectionSegmentPlan
  31 //=================================================================================================
  32
  33 float MESHCUT::longueurSegment(int ngA, int ngB)
  34 {
  35   float A[3], B[3];
  36   A[0] = MAILLAGE1->XX[ngA - 1];
  37   A[1] = MAILLAGE1->YY[ngA - 1];
  38   A[2] = MAILLAGE1->ZZ[ngA - 1];
  39   B[0] = MAILLAGE1->XX[ngB - 1];
  40   B[1] = MAILLAGE1->YY[ngB - 1];
  41   B[2] = MAILLAGE1->ZZ[ngB - 1];
  42   float dx = B[0] - A[0];
  43   float dy = B[1] - A[1];
  44   float dz = B[2] - A[2];
  45   return sqrt(dx * dx + dy * dy + dz * dz);
  46 }
  47
  48 float MESHCUT::distanceNoeudPlan(float point[3])
  49 {
  50   return normale[0] * point[0] + normale[1] * point[1] + normale[2] * point[2] + d;
  51 }
  52
  53 float MESHCUT::distanceNoeudPlan(int ng)
  54 {
  55   float A[3];
  56   A[0] = MAILLAGE1->XX[ng - 1];
  57   A[1] = MAILLAGE1->YY[ng - 1];
  58   A[2] = MAILLAGE1->ZZ[ng - 1];
  59   return distanceNoeudPlan(A);
  60 }
  61
  62 int MESHCUT::positionNoeudPlan(int indiceNoeud)
  63 {
  64   if (distanceNoeudPlan(indiceNoeud + 1) > epsilon)
  65     return 1;
  66   else if (distanceNoeudPlan(indiceNoeud + 1) < -epsilon)
  67     return -1;
  68   else
  69     return 0;
  70 }
  71
  72 /*!
  73  * Equation paramétrique de la droite AB:    OP = OA + lambda AB
  74  *
  75  * Fonction caractéristique du plan : PI(X,Y,Z) = nx X + ny Y + nz Z + d
  76  *
  77  * Pour un point P de la droite: PI(OP) = PI(OA) + lambda n.AB     avec n=(nx,ny,nz),
  78  * L'intersection AB/plan est donnée par le point P tel que PI(OP)=0.
  79  *
  80  * Il lui correspond la coordonnée   lambda = - PI(OA) / n.AB.
  81  *
  82  * Cette intersection est dans le segment si lambda est dans [0,1]
  83  */
  84 int MESHCUT::intersectionSegmentPlan(int it4, int na)
  85 {
  86
  87   int ngA = -1, ngB = -1; // Numéros des noeuds extrémités AB
  88   float lambda, ps; //, ab; // ab = longueur AB
  89   float A[3], B[3];
  90
  91   // Détermination des ng des extrémités de l'arête passée en argument na
  92   med_int * offset = MAILLAGE1->CNX[TETRA4] + 4 * it4;
  93   if (na == 0)
  94     {
  95       ngA = *(offset + 0);
  96       ngB = *(offset + 1);
  97     }
  98   else if (na == 1)
  99     {
 100       ngA = *(offset + 0);
 101       ngB = *(offset + 2);
 102     }
 103   else if (na == 2)
 104     {
 105       ngA = *(offset + 0);
 106       ngB = *(offset + 3);
 107     }
 108   else if (na == 3)
 109     {
 110       ngA = *(offset + 1);
 111       ngB = *(offset + 2);
 112     }
 113   else if (na == 4)
 114     {
 115       ngA = *(offset + 1);
 116       ngB = *(offset + 3);
 117     }
 118   else if (na == 5)
 119     {
 120       ngA = *(offset + 2);
 121       ngB = *(offset + 3);
 122     }
 123   else
 124     ERREUR("Edge number superior to 6");
 125
 126   string cle1 = int2string(ngA) + (string) "_" + int2string(ngB);
 127   string cle2 = int2string(ngB) + (string) "_" + int2string(ngA);
 128
 129   if (intersections[cle1])
 130     return intersections[cle1];
 131
 132   else
 133     {
 134
 135       A[0] = MAILLAGE1->XX[ngA - 1];
 136       A[1] = MAILLAGE1->YY[ngA - 1];
 137       A[2] = MAILLAGE1->ZZ[ngA - 1];
 138       B[0] = MAILLAGE1->XX[ngB - 1];
 139       B[1] = MAILLAGE1->YY[ngB - 1];
 140       B[2] = MAILLAGE1->ZZ[ngB - 1];
 141
 142       //      // Longueur AB
 143       //      float lx = B[0] - A[0], ly = B[1] - A[1], lz = B[2] - A[2];
 144       //      ab = sqrt(lx * lx + ly * ly + lz * lz);
 145       //      // La longueur maximale théorique est 2 epsilon
 146       //      if (ab < 2 * epsilon * 0.9)
 147       //        ERREUR("Arête de longueur inférieure au minimum théorique 2 epsilon");
 148
 149       // Calcul du produit scalaire AB.n
 150       ps = 0.0;
 151       for (int k = 0; k < 3; k++)
 152         ps += (B[k] - A[k]) * normale[k];
 153       // ps = ps / ab ;
 154
 155       if (debug)
 156         {
 157           cout << "Routine ISP : arête " << na << " -  ngA=" << ngA << " ngB=" << ngB << endl;
 158           cout << "A : " << A[0] << ' ' << A[1] << ' ' << A[2] << endl;
 159           cout << "B : " << B[0] << ' ' << B[1] << ' ' << B[2] << endl;
 160           cout << "N : " << normale[0] << ' ' << normale[1] << ' ' << normale[2] << endl;
 161         }
 162
 163       if (fabs(ps) == 0.0)
 164         ERREUR("Error on null scalar product");
 165
 166       // PS non nul: l'intersection AB/plan existe
 167
 168       lambda = -distanceNoeudPlan(A) / ps;
 169
 170       float inter[3];
 171       for (int k = 0; k < 3; k++)
 172         inter[k] = A[k] + lambda * (B[k] - A[k]);
 173       newXX.push_back(inter[0]);
 174       newYY.push_back(inter[1]);
 175       newZZ.push_back(inter[2]);
 176       indexNouveauxNoeuds++;
 177       intersections[cle1] = indexNouveauxNoeuds;
 178       intersections[cle2] = indexNouveauxNoeuds;
 179
 180       //      cout << "création noeud " << indexNouveauxNoeuds << " : " << inter[0] << " " << inter[1] << " " << inter[2]
 181       //          << endl;
 182       if (debug)
 183         cout << " sortie nouveau noeud, lambda = " << lambda << " , noeud = " << indexNouveauxNoeuds << endl;
 184       return indexNouveauxNoeuds;
 185
 186     }
 187 }
 188