src/Tools/MeshCut/MeshCut_Fonctions.cxx

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