1 // Copyright (C) 2006-2011 EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "MeshCut_Fonctions.hxx"
21 #include "MeshCut_Globals.hxx"
26 using namespace MESHCUT;
29 //=================================================================================================
30 // intersectionSegmentPlan
31 //=================================================================================================
33 float MESHCUT::longueurSegment(int ngA, int ngB)
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);
48 float MESHCUT::distanceNoeudPlan(float point[3])
50 return normale[0] * point[0] + normale[1] * point[1] + normale[2] * point[2] + d;
53 float MESHCUT::distanceNoeudPlan(int ng)
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);
62 int MESHCUT::positionNoeudPlan(int indiceNoeud)
64 if (distanceNoeudPlan(indiceNoeud + 1) > epsilon)
66 else if (distanceNoeudPlan(indiceNoeud + 1) < -epsilon)
73 * Equation paramétrique de la droite AB: OP = OA + lambda AB
75 * Fonction caractéristique du plan : PI(X,Y,Z) = nx X + ny Y + nz Z + d
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.
80 * Il lui correspond la coordonnée lambda = - PI(OA) / n.AB.
82 * Cette intersection est dans le segment si lambda est dans [0,1]
84 int MESHCUT::intersectionSegmentPlan(int it4, int na)
87 int ngA, ngB; // Numéros des noeuds extrémités AB
88 float lambda, ps; //, ab; // ab = longueur AB
91 // Détermination des ng des extrémités de l'arête passée en argument na
92 int * offset = MAILLAGE1->CNX[TETRA4] + 4 * it4;
124 ERREUR("Edge number superior to 6");
126 string cle1 = int2string(ngA) + (string) "_" + int2string(ngB);
127 string cle2 = int2string(ngB) + (string) "_" + int2string(ngA);
129 if (intersections[cle1])
130 return intersections[cle1];
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];
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");
149 // Calcul du produit scalaire AB.n
151 for (int k = 0; k < 3; k++)
152 ps += (B[k] - A[k]) * normale[k];
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;
164 ERREUR("Error on null scalar product");
166 // PS non nul: l'intersection AB/plan existe
168 lambda = -distanceNoeudPlan(A) / ps;
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;
180 // cout << "création noeud " << indexNouveauxNoeuds << " : " << inter[0] << " " << inter[1] << " " << inter[2]
183 cout << " sortie nouveau noeud, lambda = " << lambda << " , noeud = " << indexNouveauxNoeuds << endl;
184 return indexNouveauxNoeuds;