Back to index

salome-med  6.5.0
InterpKernelMeshQuality.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, 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 "InterpKernelMeshQuality.hxx"
00021 
00022 #include <cmath>
00023 #include <limits>
00024 #include <numeric>
00025 #include <algorithm>
00026 
00027 double INTERP_KERNEL::quadSkew(const double *coo)
00028 {
00029   double pa0[3]={
00030     coo[3]+coo[6]-coo[0]-coo[9],
00031     coo[4]+coo[7]-coo[1]-coo[10],
00032     coo[5]+coo[8]-coo[2]-coo[11]
00033   };
00034   double pa1[3]={
00035     coo[6]+coo[9]-coo[0]-coo[3],
00036     coo[7]+coo[10]-coo[1]-coo[4],
00037     coo[8]+coo[11]-coo[2]-coo[5],
00038   };
00039   double l0=sqrt(pa0[0]*pa0[0]+pa0[1]*pa0[1]+pa0[2]*pa0[2]);
00040   double l1=sqrt(pa1[0]*pa1[0]+pa1[1]*pa1[1]+pa1[2]*pa1[2]);
00041   if(l0<1.e-15)
00042     return 0.;
00043   if(l1<1.e-15)
00044     return 0.;
00045   pa0[0]/=l0; pa0[1]/=l0; pa0[2]/=l0;
00046   pa1[0]/=l1; pa1[1]/=l1; pa1[2]/=l1;
00047   return pa0[0]*pa1[0]+pa0[1]*pa1[1]+pa0[2]*pa1[2];
00048 }
00049 
00050 double INTERP_KERNEL::quadEdgeRatio(const double *coo)
00051 {
00052   double a2=(coo[3]-coo[0])*(coo[3]-coo[0])+(coo[4]-coo[1])*(coo[4]-coo[1])+(coo[5]-coo[2])*(coo[5]-coo[2]);
00053   double b2=(coo[6]-coo[3])*(coo[6]-coo[3])+(coo[7]-coo[4])*(coo[7]-coo[4])+(coo[8]-coo[5])*(coo[8]-coo[5]);
00054   double c2=(coo[9]-coo[6])*(coo[9]-coo[6])+(coo[10]-coo[7])*(coo[10]-coo[7])+(coo[11]-coo[8])*(coo[11]-coo[8]);
00055   double d2=(coo[0]-coo[9])*(coo[0]-coo[9])+(coo[1]-coo[10])*(coo[1]-coo[10])+(coo[2]-coo[11])*(coo[2]-coo[11]);
00056   double mab=a2<b2?a2:b2;
00057   double Mab=a2<b2?b2:a2;
00058   double mcd=c2<d2?c2:d2;
00059   double Mcd=c2<d2?d2:c2;
00060   double m2=mab<mcd?mab:mcd;
00061   double M2=Mab>Mcd?Mab:Mcd;
00062   if(m2>1.e-15)
00063     return sqrt(M2/m2);
00064   else
00065     return std::numeric_limits<double>::max();
00066 }
00067 
00068 double INTERP_KERNEL::quadAspectRatio(const double *coo)
00069 {
00070   double a=sqrt((coo[3]-coo[0])*(coo[3]-coo[0])+(coo[4]-coo[1])*(coo[4]-coo[1])+(coo[5]-coo[2])*(coo[5]-coo[2]));
00071   double b=sqrt((coo[6]-coo[3])*(coo[6]-coo[3])+(coo[7]-coo[4])*(coo[7]-coo[4])+(coo[8]-coo[5])*(coo[8]-coo[5]));
00072   double c=sqrt((coo[9]-coo[6])*(coo[9]-coo[6])+(coo[10]-coo[7])*(coo[10]-coo[7])+(coo[11]-coo[8])*(coo[11]-coo[8]));
00073   double d=sqrt((coo[0]-coo[9])*(coo[0]-coo[9])+(coo[1]-coo[10])*(coo[1]-coo[10])+(coo[2]-coo[11])*(coo[2]-coo[11]));
00074   double ma=a>b?a:b;
00075   double mb=c>d?c:d;
00076   double hm=ma>mb?ma:mb;
00077   double ab[3]={(coo[4]-coo[1])*(coo[8]-coo[5])-(coo[7]-coo[4])*(coo[5]-coo[2]),
00078                 (coo[5]-coo[2])*(coo[6]-coo[3])-(coo[3]-coo[0])*(coo[8]-coo[5]),
00079                 (coo[3]-coo[0])*(coo[7]-coo[4])-(coo[4]-coo[1])*(coo[6]-coo[3])};
00080   double cd[3]={(coo[10]-coo[7])*(coo[2]-coo[11])-(coo[1]-coo[10])*(coo[11]-coo[8]),
00081                 (coo[11]-coo[8])*(coo[0]-coo[9])-(coo[9]-coo[6])*(coo[2]-coo[11]),
00082                 (coo[9]-coo[6])*(coo[1]-coo[10])-(coo[10]-coo[7])*(coo[0]-coo[9])};
00083   double e=sqrt(ab[0]*ab[0]+ab[1]*ab[1]+ab[2]*ab[2])+sqrt(cd[0]*cd[0]+cd[1]*cd[1]+cd[2]*cd[2]);
00084   if(d>1e-15)
00085     return 0.5*(a+b+c+d)*hm/e;
00086   else
00087     return std::numeric_limits<double>::max();
00088 }
00089 
00090 double INTERP_KERNEL::quadWarp(const double *coo)
00091 {
00092   double e0[3]={coo[3]-coo[0],coo[4]-coo[1],coo[5]-coo[2]};
00093   double e1[3]={coo[6]-coo[3],coo[7]-coo[4],coo[8]-coo[5]};
00094   double e2[3]={coo[9]-coo[6],coo[10]-coo[7],coo[11]-coo[8]};
00095   double e3[3]={coo[0]-coo[9],coo[1]-coo[10],coo[2]-coo[11]};
00096   
00097   double n0[3]={e3[1]*e0[2]-e3[2]*e0[1],e3[2]*e0[0]-e3[0]*e0[2],e3[0]*e0[1]-e3[1]*e0[0]};
00098   double n1[3]={e0[1]*e1[2]-e0[2]*e1[1],e0[2]*e1[0]-e0[0]*e1[2],e0[0]*e1[1]-e0[1]*e1[0]};
00099   double n2[3]={e1[1]*e2[2]-e1[2]*e2[1],e1[2]*e2[0]-e1[0]*e2[2],e1[0]*e2[1]-e1[1]*e2[0]};
00100   double n3[3]={e2[1]*e3[2]-e2[2]*e3[1],e2[2]*e3[0]-e2[0]*e3[2],e2[0]*e3[1]-e2[1]*e3[0]};
00101 
00102   double l0=sqrt(n0[0]*n0[0]+n0[1]*n0[1]+n0[2]*n0[2]);
00103   double l1=sqrt(n1[0]*n1[0]+n1[1]*n1[1]+n1[2]*n1[2]);
00104   double l2=sqrt(n2[0]*n2[0]+n2[1]*n2[1]+n2[2]*n2[2]);
00105   double l3=sqrt(n3[0]*n3[0]+n3[1]*n3[1]+n3[2]*n3[2]);
00106 
00107   if(l0<1.e-15 || l1<1.e-15 || l2<1.e-15 || l3<1e-15)
00108     return std::numeric_limits<double>::min();
00109 
00110   double warp=std::min(n0[0]/l0*n2[0]/l2+n0[1]/l0*n2[1]/l2+n0[2]/l0*n2[2]/l2,n1[0]/l1*n3[0]/l3+n1[1]/l1*n3[1]/l3+n1[2]/l1*n3[2]/l3);
00111   return warp*warp*warp;
00112 }
00113 
00114 double INTERP_KERNEL::triEdgeRatio(const double *coo)
00115 {
00116   double a2=(coo[3]-coo[0])*(coo[3]-coo[0])+(coo[4]-coo[1])*(coo[4]-coo[1])+(coo[5]-coo[2])*(coo[5]-coo[2]);
00117   double b2=(coo[6]-coo[3])*(coo[6]-coo[3])+(coo[7]-coo[4])*(coo[7]-coo[4])+(coo[8]-coo[5])*(coo[8]-coo[5]);
00118   double c2=(coo[0]-coo[6])*(coo[0]-coo[6])+(coo[1]-coo[7])*(coo[1]-coo[7])+(coo[2]-coo[8])*(coo[2]-coo[8]);
00119   double mab=a2<b2?a2:b2;
00120   double Mab=a2<b2?b2:a2;
00121   double m2=c2>mab?mab:c2;
00122   double M2=c2>Mab?c2:Mab;
00123   if(m2>1.e-15)
00124     return sqrt(M2/m2);
00125   else
00126     return std::numeric_limits<double>::max();
00127 }
00128 
00129 double INTERP_KERNEL::triAspectRatio(const double *coo)
00130 {
00131   double a=sqrt((coo[3]-coo[0])*(coo[3]-coo[0])+(coo[4]-coo[1])*(coo[4]-coo[1])+(coo[5]-coo[2])*(coo[5]-coo[2]));
00132   double b=sqrt((coo[6]-coo[3])*(coo[6]-coo[3])+(coo[7]-coo[4])*(coo[7]-coo[4])+(coo[8]-coo[5])*(coo[8]-coo[5]));
00133   double c=sqrt((coo[0]-coo[6])*(coo[0]-coo[6])+(coo[1]-coo[7])*(coo[1]-coo[7])+(coo[2]-coo[8])*(coo[2]-coo[8]));
00134  
00135   double hm=a>b?a:b;
00136   hm=hm>c?hm:c;
00137 
00138   double ab[3]={(coo[4]-coo[1])*(coo[8]-coo[5])-(coo[7]-coo[4])*(coo[5]-coo[2]),
00139                 (coo[5]-coo[2])*(coo[6]-coo[3])-(coo[3]-coo[0])*(coo[8]-coo[5]),
00140                 (coo[3]-coo[0])*(coo[7]-coo[4])-(coo[4]-coo[1])*(coo[6]-coo[3])};
00141   double d=sqrt(ab[0]*ab[0]+ab[1]*ab[1]+ab[2]*ab[2]);
00142   static const double normalizeCoeff=sqrt(3.)/6.;
00143   if(d>1.e-15) 
00144     return normalizeCoeff*hm*(a+b+c)/d;
00145   else
00146     return std::numeric_limits<double>::max();
00147 }
00148 
00149 double INTERP_KERNEL::tetraEdgeRatio(const double *coo)
00150 {
00151   double a[3]={coo[3]-coo[0],coo[4]-coo[1],coo[5]-coo[2]};
00152   double b[3]={coo[6]-coo[3],coo[7]-coo[4],coo[8]-coo[5]};
00153   double c[3]={coo[0]-coo[6],coo[1]-coo[7],coo[2]-coo[8]};
00154   double d[3]={coo[9]-coo[0],coo[10]-coo[1],coo[11]-coo[2]};
00155   double e[3]={coo[9]-coo[3],coo[10]-coo[4],coo[11]-coo[5]};
00156   double f[3]={coo[9]-coo[6],coo[10]-coo[7],coo[11]-coo[8]};
00157   
00158   double l2[6]=
00159     {a[0]*a[0]+a[1]*a[1]+a[2]*a[2],
00160      b[0]*b[0]+b[1]*b[1]+b[2]*b[2],
00161      c[0]*c[0]+c[1]*c[1]+c[2]*c[2],
00162      d[0]*d[0]+d[1]*d[1]+d[2]*d[2],
00163      e[0]*e[0]+e[1]*e[1]+e[2]*e[2],
00164      f[0]*f[0]+f[1]*f[1]+f[2]*f[2]};
00165 
00166   double M2=*std::max_element(l2,l2+6);
00167   double m2=*std::min_element(l2,l2+6);
00168   if(m2>1e-15)
00169     return sqrt(M2/m2);
00170   else
00171     return std::numeric_limits<double>::max();
00172 }
00173 
00174 double INTERP_KERNEL::tetraAspectRatio(const double *coo)
00175 {
00176   static const double normalizeCoeff=sqrt(6.)/12.;
00177   double ab[3]={coo[3]-coo[0],coo[4]-coo[1],coo[5]-coo[2]};
00178   double ac[3]={coo[6]-coo[0],coo[7]-coo[1],coo[8]-coo[2]};
00179   double ad[3]={coo[9]-coo[0],coo[10]-coo[1],coo[11]-coo[2]};
00180   double detTet=(ab[0]*(ac[1]*ad[2]-ac[2]*ad[1]))+(ab[1]*(ac[2]*ad[0]-ac[0]*ad[2]))+(ab[2]*(ac[0]*ad[1]-ac[1]*ad[0]));
00181   //if(detTet<1.e-15)
00182   //  return std::numeric_limits<double>::max();
00183   double bc[3]={coo[6]-coo[3],coo[7]-coo[4],coo[8]-coo[5]};
00184   double bd[3]={coo[9]-coo[3],coo[10]-coo[4],coo[11]-coo[5]};
00185   double cd[3]={coo[9]-coo[6],coo[10]-coo[7],coo[11]-coo[8]};
00186 
00187   double ab2=ab[0]*ab[0]+ab[1]*ab[1]+ab[2]*ab[2];
00188   double bc2=bc[0]*bc[0]+bc[1]*bc[1]+bc[2]*bc[2];
00189   double ac2=ac[0]*ac[0]+ac[1]*ac[1]+ac[2]*ac[2];
00190   double ad2=ad[0]*ad[0]+ad[1]*ad[1]+ad[2]*ad[2];
00191   double bd2=bd[0]*bd[0]+bd[1]*bd[1]+bd[2]*bd[2];
00192   double cd2=cd[0]*cd[0]+cd[1]*cd[1]+cd[2]*cd[2];
00193 
00194   double A=ab2>bc2?ab2:bc2;
00195   double B=ac2>ad2?ac2:ad2;
00196   double C=bd2>cd2?bd2:cd2;
00197   double D=A>B?A:B;
00198   double hm=D>C?sqrt(D):sqrt(C);
00199 
00200   bd[0]=ab[1]*bc[2]-ab[2]*bc[1]; bd[1]=ab[2]*bc[0]-ab[0]*bc[2]; bd[2]=ab[0]*bc[1]-ab[1]*bc[0];
00201   A=sqrt(bd[0]*bd[0]+bd[1]*bd[1]+bd[2]*bd[2]);
00202   bd[0]=ab[1]*ad[2]-ab[2]*ad[1]; bd[1]=ab[2]*ad[0]-ab[0]*ad[2]; bd[2]=ab[0]*ad[1]-ab[1]*ad[0];
00203   B=sqrt(bd[0]*bd[0]+bd[1]*bd[1]+bd[2]*bd[2]);
00204   bd[0]=ac[1]*ad[2]-ac[2]*ad[1]; bd[1]=ac[2]*ad[0]-ac[0]*ad[2]; bd[2]=ac[0]*ad[1]-ac[1]*ad[0];
00205   C=sqrt(bd[0]*bd[0]+bd[1]*bd[1]+bd[2]*bd[2]);
00206   bd[0]=bc[1]*cd[2]-bc[2]*cd[1]; bd[1]=bc[2]*cd[0]-bc[0]*cd[2]; bd[2]=bc[0]*cd[1]-bc[1]*cd[0];
00207   D=sqrt(bd[0]*bd[0]+bd[1]*bd[1]+bd[2]*bd[2]);
00208   return normalizeCoeff*hm*(A+B+C+D)/fabs(detTet);
00209 }