Back to index

salome-smesh  6.5.0
MeshCut_Fonctions.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2006-2012  EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MeshCut_Fonctions.hxx"
00021 #include "MeshCut_Globals.hxx"
00022 
00023 #include <iostream>
00024 #include <cmath>
00025 
00026 using namespace MESHCUT;
00027 using namespace std;
00028 
00029 //=================================================================================================
00030 //                                    intersectionSegmentPlan
00031 //=================================================================================================
00032 
00033 float MESHCUT::longueurSegment(int ngA, int ngB)
00034 {
00035   float A[3], B[3];
00036   A[0] = MAILLAGE1->XX[ngA - 1];
00037   A[1] = MAILLAGE1->YY[ngA - 1];
00038   A[2] = MAILLAGE1->ZZ[ngA - 1];
00039   B[0] = MAILLAGE1->XX[ngB - 1];
00040   B[1] = MAILLAGE1->YY[ngB - 1];
00041   B[2] = MAILLAGE1->ZZ[ngB - 1];
00042   float dx = B[0] - A[0];
00043   float dy = B[1] - A[1];
00044   float dz = B[2] - A[2];
00045   return sqrt(dx * dx + dy * dy + dz * dz);
00046 }
00047 
00048 float MESHCUT::distanceNoeudPlan(float point[3])
00049 {
00050   return normale[0] * point[0] + normale[1] * point[1] + normale[2] * point[2] + d;
00051 }
00052 
00053 float MESHCUT::distanceNoeudPlan(int ng)
00054 {
00055   float A[3];
00056   A[0] = MAILLAGE1->XX[ng - 1];
00057   A[1] = MAILLAGE1->YY[ng - 1];
00058   A[2] = MAILLAGE1->ZZ[ng - 1];
00059   return distanceNoeudPlan(A);
00060 }
00061 
00062 int MESHCUT::positionNoeudPlan(int indiceNoeud)
00063 {
00064   if (distanceNoeudPlan(indiceNoeud + 1) > epsilon)
00065     return 1;
00066   else if (distanceNoeudPlan(indiceNoeud + 1) < -epsilon)
00067     return -1;
00068   else
00069     return 0;
00070 }
00071 
00084 int MESHCUT::intersectionSegmentPlan(int it4, int na)
00085 {
00086 
00087   int ngA, ngB; // Numéros des noeuds extrémités AB
00088   float lambda, ps; //, ab; // ab = longueur AB
00089   float A[3], B[3];
00090 
00091   // Détermination des ng des extrémités de l'arête passée en argument na
00092   int * offset = MAILLAGE1->CNX[TETRA4] + 4 * it4;
00093   if (na == 0)
00094     {
00095       ngA = *(offset + 0);
00096       ngB = *(offset + 1);
00097     }
00098   else if (na == 1)
00099     {
00100       ngA = *(offset + 0);
00101       ngB = *(offset + 2);
00102     }
00103   else if (na == 2)
00104     {
00105       ngA = *(offset + 0);
00106       ngB = *(offset + 3);
00107     }
00108   else if (na == 3)
00109     {
00110       ngA = *(offset + 1);
00111       ngB = *(offset + 2);
00112     }
00113   else if (na == 4)
00114     {
00115       ngA = *(offset + 1);
00116       ngB = *(offset + 3);
00117     }
00118   else if (na == 5)
00119     {
00120       ngA = *(offset + 2);
00121       ngB = *(offset + 3);
00122     }
00123   else
00124     ERREUR("Edge number superior to 6");
00125 
00126   string cle1 = int2string(ngA) + (string) "_" + int2string(ngB);
00127   string cle2 = int2string(ngB) + (string) "_" + int2string(ngA);
00128 
00129   if (intersections[cle1])
00130     return intersections[cle1];
00131 
00132   else
00133     {
00134 
00135       A[0] = MAILLAGE1->XX[ngA - 1];
00136       A[1] = MAILLAGE1->YY[ngA - 1];
00137       A[2] = MAILLAGE1->ZZ[ngA - 1];
00138       B[0] = MAILLAGE1->XX[ngB - 1];
00139       B[1] = MAILLAGE1->YY[ngB - 1];
00140       B[2] = MAILLAGE1->ZZ[ngB - 1];
00141 
00142       //      // Longueur AB
00143       //      float lx = B[0] - A[0], ly = B[1] - A[1], lz = B[2] - A[2];
00144       //      ab = sqrt(lx * lx + ly * ly + lz * lz);
00145       //      // La longueur maximale théorique est 2 epsilon
00146       //      if (ab < 2 * epsilon * 0.9)
00147       //        ERREUR("Arête de longueur inférieure au minimum théorique 2 epsilon");
00148 
00149       // Calcul du produit scalaire AB.n
00150       ps = 0.0;
00151       for (int k = 0; k < 3; k++)
00152         ps += (B[k] - A[k]) * normale[k];
00153       // ps = ps / ab ;
00154 
00155       if (debug)
00156         {
00157           cout << "Routine ISP : arête " << na << " -  ngA=" << ngA << " ngB=" << ngB << endl;
00158           cout << "A : " << A[0] << ' ' << A[1] << ' ' << A[2] << endl;
00159           cout << "B : " << B[0] << ' ' << B[1] << ' ' << B[2] << endl;
00160           cout << "N : " << normale[0] << ' ' << normale[1] << ' ' << normale[2] << endl;
00161         }
00162 
00163       if (fabs(ps) == 0.0)
00164         ERREUR("Error on null scalar product");
00165 
00166       // PS non nul: l'intersection AB/plan existe
00167 
00168       lambda = -distanceNoeudPlan(A) / ps;
00169 
00170       float inter[3];
00171       for (int k = 0; k < 3; k++)
00172         inter[k] = A[k] + lambda * (B[k] - A[k]);
00173       newXX.push_back(inter[0]);
00174       newYY.push_back(inter[1]);
00175       newZZ.push_back(inter[2]);
00176       indexNouveauxNoeuds++;
00177       intersections[cle1] = indexNouveauxNoeuds;
00178       intersections[cle2] = indexNouveauxNoeuds;
00179 
00180       //      cout << "création noeud " << indexNouveauxNoeuds << " : " << inter[0] << " " << inter[1] << " " << inter[2]
00181       //          << endl;
00182       if (debug)
00183         cout << " sortie nouveau noeud, lambda = " << lambda << " , noeud = " << indexNouveauxNoeuds << endl;
00184       return indexNouveauxNoeuds;
00185 
00186     }
00187 }
00188