1 #include "MeshCut_Fonctions.hxx"
2 #include "MeshCut_Globals.hxx"
7 using namespace MESHCUT;
10 //=================================================================================================
11 // intersectionSegmentPlan
12 //=================================================================================================
14 float MESHCUT::longueurSegment(int ngA, int ngB)
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);
29 float MESHCUT::distanceNoeudPlan(float point[3])
31 return normale[0] * point[0] + normale[1] * point[1] + normale[2] * point[2] + d;
34 float MESHCUT::distanceNoeudPlan(int ng)
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);
43 int MESHCUT::positionNoeudPlan(int indiceNoeud)
45 if (distanceNoeudPlan(indiceNoeud + 1) > epsilon)
47 else if (distanceNoeudPlan(indiceNoeud + 1) < -epsilon)
54 * Equation paramétrique de la droite AB: OP = OA + lambda AB
56 * Fonction caractéristique du plan : PI(X,Y,Z) = nx X + ny Y + nz Z + d
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.
61 * Il lui correspond la coordonnée lambda = - PI(OA) / n.AB.
63 * Cette intersection est dans le segment si lambda est dans [0,1]
65 int MESHCUT::intersectionSegmentPlan(int it4, int na)
68 int ngA, ngB; // Numéros des noeuds extrémités AB
69 float lambda, ps; //, ab; // ab = longueur AB
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;
105 ERREUR("Edge number superior to 6");
107 string cle1 = int2string(ngA) + (string) "_" + int2string(ngB);
108 string cle2 = int2string(ngB) + (string) "_" + int2string(ngA);
110 if (intersections[cle1])
111 return intersections[cle1];
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];
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");
130 // Calcul du produit scalaire AB.n
132 for (int k = 0; k < 3; k++)
133 ps += (B[k] - A[k]) * normale[k];
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;
145 ERREUR("Error on null scalar product");
147 // PS non nul: l'intersection AB/plan existe
149 lambda = -distanceNoeudPlan(A) / ps;
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;
161 // cout << "création noeud " << indexNouveauxNoeuds << " : " << inter[0] << " " << inter[1] << " " << inter[2]
164 cout << " sortie nouveau noeud, lambda = " << lambda << " , noeud = " << indexNouveauxNoeuds << endl;
165 return indexNouveauxNoeuds;