src/INTERP_KERNEL/DiameterCalculator.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
   2 //
   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, or (at your option) any later version.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "DiameterCalculator.hxx"
  22 #include "InterpKernelException.hxx"
  23 #include "CellModel.hxx"
  24
  25 #include <algorithm>
  26 #include <sstream>
  27 #include <cmath>
  28
  29 using namespace INTERP_KERNEL;
  30
  31 NormalizedCellType DiameterCalulatorTRI3S2::TYPE=NORM_TRI3;
  32
  33 NormalizedCellType DiameterCalulatorTRI3S3::TYPE=NORM_TRI3;
  34
  35 NormalizedCellType DiameterCalulatorQUAD4S2::TYPE=NORM_QUAD4;
  36
  37 NormalizedCellType DiameterCalulatorQUAD4S3::TYPE=NORM_QUAD4;
  38
  39 NormalizedCellType DiameterCalulatorTETRA4::TYPE=NORM_TETRA4;
  40
  41 NormalizedCellType DiameterCalulatorHEXA8::TYPE=NORM_HEXA8;
  42
  43 NormalizedCellType DiameterCalulatorPENTA6::TYPE=NORM_PENTA6;
  44
  45 NormalizedCellType DiameterCalulatorPYRA5::TYPE=NORM_PYRA5;
  46
  47 inline double SqNormV2(const double tab[2])
  48 {
  49   return tab[0]*tab[0]+tab[1]*tab[1];
  50 }
  51
  52 inline void DiffV2(const double a[2], const double b[2], double c[2])
  53 {
  54   c[0]=a[0]-b[0];
  55   c[1]=a[1]-b[1];
  56 }
  57
  58 inline double SqNormV3(const double tab[3])
  59 {
  60   return tab[0]*tab[0]+tab[1]*tab[1]+tab[2]*tab[2];
  61 }
  62
  63 inline void DiffV3(const double a[3], const double b[3], double c[3])
  64 {
  65   c[0]=a[0]-b[0];
  66   c[1]=a[1]-b[1];
  67   c[2]=a[2]-b[2];
  68 }
  69
  70 template<class Evaluator>
  71 void ComputeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr)
  72 {
  73   Evaluator evtor;
  74   NormalizedCellType ct(Evaluator::TYPE);
  75   int cti((int) ct);
  76   for(const int *it=bgIds;it!=endIds;it++)
  77     {
  78       int offset(indPtr[*it]);
  79       if(connPtr[offset]==cti)
  80         resPtr[*it]=evtor.computeForOneCellInternal(connPtr+offset+1,connPtr+indPtr[(*it)+1],coordsPtr);
  81       else
  82         {
  83           std::ostringstream oss; oss << "DiameterCalculator::computeForListOfCellIdsUMeshFrmt : invalid nodal connectivity format at cell # " << *it <<  " !";
  84           throw Exception(oss.str().c_str());
  85         }
  86     }
  87 }
  88
  89 template<class Evaluator>
  90 void ComputeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr)
  91 {
  92   Evaluator evtor;
  93   NormalizedCellType ct(Evaluator::TYPE);
  94   int cti((int) ct);
  95   for(int it=bgId;it<endId;it++)
  96     {
  97       int offset(indPtr[it]);
  98       if(connPtr[offset]==cti)
  99         resPtr[it]=evtor.computeForOneCellInternal(connPtr+offset+1,connPtr+indPtr[it+1],coordsPtr);
 100       else
 101         {
 102           std::ostringstream oss; oss << "DiameterCalculator::computeForListOfCellIdsUMeshFrmt : invalid nodal connectivity format at cell # " << it <<  " !";
 103           throw Exception(oss.str().c_str());
 104         }
 105     }
 106 }
 107
 108 template<class Evaluator>
 109 void ComputeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr)
 110 {
 111   Evaluator evtor;
 112   NormalizedCellType ct(Evaluator::TYPE);
 113   const CellModel& cm(CellModel::GetCellModel(ct));
 114   unsigned nbNodes(cm.getNumberOfNodes());
 115   const int *ptr(connPtr);
 116   for(int i=0;i<nbOfCells;i++,ptr+=nbNodes,resPtr++)
 117     *resPtr=evtor.computeForOneCellInternal(ptr,ptr+nbNodes,coordsPtr);
 118 }
 119
 120 double DiameterCalulatorTRI3S2::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 121 {
 122   if(std::distance(bg,endd)==3)
 123     {
 124       const double *a(coordsPtr+2*bg[0]),*b(coordsPtr+2*bg[1]),*c(coordsPtr+2*bg[2]);
 125       double l0[2],l1[2],l2[2];
 126       DiffV2(a,b,l0);
 127       DiffV2(a,c,l1);
 128       DiffV2(b,c,l2);
 129       double res(std::max(SqNormV2(l0),SqNormV2(l1)));
 130       res=std::max(res,SqNormV2(l2));
 131       return std::sqrt(res);
 132     }
 133   else
 134     throw Exception("DiameterCalulatorTRI3S2::computeForOneCellInternal : input connectivity must be of size 3 !");
 135 }
 136
 137 void DiameterCalulatorTRI3S2::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 138 {
 139   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S2>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 140 }
 141
 142 void DiameterCalulatorTRI3S2::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 143 {
 144   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S2>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 145 }
 146
 147 void DiameterCalulatorTRI3S2::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 148 {
 149   ComputeFor1SGTUMeshFrmt<DiameterCalulatorTRI3S2>(nbOfCells,connPtr,coordsPtr,resPtr);
 150 }
 151
 152 double DiameterCalulatorTRI3S3::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 153 {
 154   if(std::distance(bg,endd)==3)
 155     {
 156       const double *a(coordsPtr+3*bg[0]),*b(coordsPtr+3*bg[1]),*c(coordsPtr+3*bg[2]);
 157       double l0[3],l1[3],l2[3];
 158       DiffV3(a,b,l0);
 159       DiffV3(a,c,l1);
 160       DiffV3(b,c,l2);
 161       double res(std::max(SqNormV3(l0),SqNormV3(l1)));
 162       res=std::max(res,SqNormV3(l2));
 163       return std::sqrt(res);
 164     }
 165   else
 166     throw Exception("DiameterCalulatorTRI3S2::computeForOneCellInternal : input connectivity must be of size 3 !");
 167 }
 168
 169 void DiameterCalulatorTRI3S3::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 170 {
 171   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S3>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 172 }
 173
 174 void DiameterCalulatorTRI3S3::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 175 {
 176   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S3>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 177 }
 178
 179 void DiameterCalulatorTRI3S3::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 180 {
 181   ComputeFor1SGTUMeshFrmt<DiameterCalulatorTRI3S3>(nbOfCells,connPtr,coordsPtr,resPtr);
 182 }
 183
 184 double DiameterCalulatorQUAD4S2::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 185 {
 186   if(std::distance(bg,endd)==4)
 187     {
 188       const double *a(coordsPtr+2*bg[0]),*b(coordsPtr+2*bg[1]),*c(coordsPtr+2*bg[2]),*d(coordsPtr+2*bg[3]);
 189       double l0[2],l1[2];
 190       DiffV2(a,c,l0);
 191       DiffV2(b,d,l1);
 192       return std::sqrt(std::max(SqNormV2(l0),SqNormV2(l1)));
 193     }
 194   else
 195     throw Exception("DiameterCalulatorQUAD4S2::computeForOneCellInternal : input connectivity must be of size 4 !");
 196 }
 197
 198 void DiameterCalulatorQUAD4S2::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 199 {
 200   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S2>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 201 }
 202
 203 void DiameterCalulatorQUAD4S2::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 204 {
 205   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S2>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 206 }
 207
 208 void DiameterCalulatorQUAD4S2::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 209 {
 210   ComputeFor1SGTUMeshFrmt<DiameterCalulatorQUAD4S2>(nbOfCells,connPtr,coordsPtr,resPtr);
 211 }
 212
 213 double DiameterCalulatorQUAD4S3::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 214 {
 215   if(std::distance(bg,endd)==4)
 216     {
 217       const double *a(coordsPtr+3*bg[0]),*b(coordsPtr+3*bg[1]),*c(coordsPtr+3*bg[2]),*d(coordsPtr+3*bg[3]);
 218       double l0[3],l1[3];
 219       DiffV3(a,c,l0);
 220       DiffV3(b,d,l1);
 221       return std::sqrt(std::max(SqNormV3(l0),SqNormV3(l1)));
 222     }
 223   else
 224     throw Exception("DiameterCalulatorQUAD4S3::computeForOneCellInternal : input connectivity must be of size 4 !");
 225 }
 226
 227 void DiameterCalulatorQUAD4S3::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 228 {
 229   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S3>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 230 }
 231
 232 void DiameterCalulatorQUAD4S3::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 233 {
 234   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S3>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 235 }
 236
 237 void DiameterCalulatorQUAD4S3::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 238 {
 239   ComputeFor1SGTUMeshFrmt<DiameterCalulatorQUAD4S3>(nbOfCells,connPtr,coordsPtr,resPtr);
 240 }
 241
 242 double DiameterCalulatorTETRA4::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 243 {
 244   if(std::distance(bg,endd)==4)
 245     {
 246       const double *a(coordsPtr+3*bg[0]),*b(coordsPtr+3*bg[1]),*c(coordsPtr+3*bg[2]),*d(coordsPtr+3*bg[3]);
 247       double l0[3],l1[3],l2[3],l3[3],l4[3],l5[3];
 248       DiffV3(a,b,l0);
 249       DiffV3(a,c,l1);
 250       DiffV3(b,c,l2);
 251       DiffV3(a,d,l3);
 252       DiffV3(b,d,l4);
 253       DiffV3(c,d,l5);
 254       double tmp[6];
 255       tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3); tmp[4]=SqNormV3(l4); tmp[5]=SqNormV3(l5);
 256       return std::sqrt(*std::max_element(tmp,tmp+6));
 257     }
 258   else
 259     throw Exception("DiameterCalulatorTETRA4::computeForOneCellInternal : input connectivity must be of size 4 !");
 260 }
 261
 262 void DiameterCalulatorTETRA4::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 263 {
 264   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorTETRA4>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 265 }
 266
 267 void DiameterCalulatorTETRA4::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 268 {
 269   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorTETRA4>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 270 }
 271
 272 void DiameterCalulatorTETRA4::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 273 {
 274   ComputeFor1SGTUMeshFrmt<DiameterCalulatorTETRA4>(nbOfCells,connPtr,coordsPtr,resPtr);
 275 }
 276
 277 double DiameterCalulatorHEXA8::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 278 {
 279   if(std::distance(bg,endd)==8)
 280     {
 281       const double *p0(coordsPtr+3*bg[0]),*p1(coordsPtr+3*bg[1]),*p2(coordsPtr+3*bg[2]),*p3(coordsPtr+3*bg[3]),*p4(coordsPtr+3*bg[4]),*p5(coordsPtr+3*bg[5]),*p6(coordsPtr+3*bg[6]),*p7(coordsPtr+3*bg[7]);
 282       double l0[3],l1[3],l2[3],l3[3];
 283       DiffV3(p0,p6,l0);
 284       DiffV3(p1,p7,l1);
 285       DiffV3(p2,p4,l2);
 286       DiffV3(p3,p5,l3);
 287       double tmp[4];
 288       tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3);
 289       return std::sqrt(*std::max_element(tmp,tmp+4));
 290     }
 291   else
 292     throw Exception("DiameterCalulatorHEXA8::computeForOneCellInternal : input connectivity must be of size 8 !");
 293 }
 294
 295 void DiameterCalulatorHEXA8::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 296 {
 297   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorHEXA8>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 298 }
 299
 300 void DiameterCalulatorHEXA8::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 301 {
 302   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorHEXA8>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 303 }
 304
 305 void DiameterCalulatorHEXA8::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 306 {
 307   ComputeFor1SGTUMeshFrmt<DiameterCalulatorHEXA8>(nbOfCells,connPtr,coordsPtr,resPtr);
 308 }
 309
 310 double DiameterCalulatorPENTA6::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 311 {
 312   if(std::distance(bg,endd)==6)
 313     {
 314       const double *p0(coordsPtr+3*bg[0]),*p1(coordsPtr+3*bg[1]),*p2(coordsPtr+3*bg[2]),*p3(coordsPtr+3*bg[3]),*p4(coordsPtr+3*bg[4]),*p5(coordsPtr+3*bg[5]);
 315       double l0[3],l1[3],l2[3],l3[3],l4[3],l5[3];
 316       DiffV3(p0,p4,l0);
 317       DiffV3(p1,p3,l1);
 318       DiffV3(p1,p5,l2);
 319       DiffV3(p2,p4,l3);
 320       DiffV3(p0,p5,l4);
 321       DiffV3(p2,p3,l5);
 322       double tmp[6];
 323       tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3); tmp[4]=SqNormV3(l4); tmp[5]=SqNormV3(l5);
 324       return std::sqrt(*std::max_element(tmp,tmp+6));
 325     }
 326   else
 327     throw Exception("DiameterCalulatorPENTA6::computeForOneCellInternal : input connectivity must be of size 6 !");
 328 }
 329
 330 void DiameterCalulatorPENTA6::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 331 {
 332   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorPENTA6>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 333 }
 334
 335 void DiameterCalulatorPENTA6::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 336 {
 337   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorPENTA6>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 338 }
 339
 340 void DiameterCalulatorPENTA6::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 341 {
 342   ComputeFor1SGTUMeshFrmt<DiameterCalulatorPENTA6>(nbOfCells,connPtr,coordsPtr,resPtr);
 343 }
 344
 345 double DiameterCalulatorPYRA5::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
 346 {
 347   if(std::distance(bg,endd)==5)
 348     {
 349       const double *p0(coordsPtr+3*bg[0]),*p1(coordsPtr+3*bg[1]),*p2(coordsPtr+3*bg[2]),*p3(coordsPtr+3*bg[3]),*p4(coordsPtr+3*bg[4]);
 350       double l0[3],l1[3],l2[3],l3[3],l4[3],l5[3];
 351       DiffV3(p0,p2,l0);
 352       DiffV3(p1,p3,l1);
 353       DiffV3(p0,p4,l2);
 354       DiffV3(p1,p4,l3);
 355       DiffV3(p2,p4,l4);
 356       DiffV3(p3,p4,l5);
 357       double tmp[6];
 358       tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3); tmp[4]=SqNormV3(l4); tmp[5]=SqNormV3(l5);
 359       return std::sqrt(*std::max_element(tmp,tmp+6));
 360     }
 361   else
 362     throw Exception("DiameterCalulatorPYRA5::computeForOneCellInternal : input connectivity must be of size 5 !");
 363 }
 364
 365 void DiameterCalulatorPYRA5::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 366 {
 367   ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorPYRA5>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
 368 }
 369
 370 void DiameterCalulatorPYRA5::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
 371 {
 372   ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorPYRA5>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
 373 }
 374
 375 void DiameterCalulatorPYRA5::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
 376 {
 377   ComputeFor1SGTUMeshFrmt<DiameterCalulatorPYRA5>(nbOfCells,connPtr,coordsPtr,resPtr);
 378 }