src/MEDWrapper/MED_GaussUtils.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "MED_GaussUtils.hxx"
  24 #include "MED_Utilities.hxx"
  25
  26 #ifdef _DEBUG_
  27 static int MYDEBUG = 0;
  28 static int MYVALUEDEBUG = 0;
  29 #else
  30 static int MYDEBUG = 0;
  31 static int MYVALUEDEBUG = 0;
  32 #endif
  33
  34 namespace MED
  35 {
  36   //---------------------------------------------------------------
  37   TGaussCoord
  38   ::TGaussCoord():
  39     TModeSwitchInfo(eFULL_INTERLACE),
  40     myNbElem(0),
  41     myNbGauss(0),
  42     myDim(0),
  43     myGaussStep(0)
  44   {
  45   }
  46
  47   void
  48   TGaussCoord
  49   ::Init(TInt theNbElem,
  50          TInt theNbGauss,
  51          TInt theDim,
  52          EModeSwitch theMode)
  53   {
  54     myModeSwitch = theMode;
  55
  56     myNbElem = theNbElem;
  57     myNbGauss = theNbGauss;
  58     myDim = theDim;
  59
  60     myGaussStep = myNbGauss*myDim;
  61
  62     myGaussCoord.resize(theNbElem*myGaussStep);
  63   }
  64
  65   TInt
  66   TGaussCoord
  67   ::GetNbElem() const
  68   {
  69     return myNbElem;
  70   }
  71
  72   TInt
  73   TGaussCoord
  74   ::GetNbGauss() const
  75   {
  76     return myNbGauss;
  77   }
  78
  79   TInt
  80   TGaussCoord
  81   ::GetDim() const
  82   {
  83     return myDim;
  84   }
  85
  86   unsigned char*
  87   TGaussCoord
  88   ::GetValuePtr()
  89   {
  90     return (unsigned char*)&(myGaussCoord[0]);
  91   }
  92
  93   TCCoordSliceArr
  94   TGaussCoord
  95   ::GetCoordSliceArr(TInt theElemId) const
  96   {
  97     TCCoordSliceArr aCoordSliceArr(myNbGauss);
  98     if(GetModeSwitch() == eFULL_INTERLACE){
  99       TInt anId = theElemId*myGaussStep;
 100       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
 101         aCoordSliceArr[anGaussId] =
 102           TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
 103         anId += myDim;
 104       }
 105     }
 106     else{
 107       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
 108         aCoordSliceArr[anGaussId] =
 109           TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
 110       }
 111     }
 112     return aCoordSliceArr;
 113   }
 114
 115   TCoordSliceArr
 116   TGaussCoord
 117   ::GetCoordSliceArr(TInt theElemId)
 118   {
 119     TCoordSliceArr aCoordSliceArr(myNbGauss);
 120     if(GetModeSwitch() == eFULL_INTERLACE){
 121       TInt anId = theElemId*myGaussStep;
 122       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
 123         aCoordSliceArr[anGaussId] =
 124           TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
 125         anId += myDim;
 126       }
 127     }
 128     else{
 129       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
 130         aCoordSliceArr[anGaussId] =
 131           TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
 132       }
 133     }
 134     return aCoordSliceArr;
 135   }
 136
 137   //---------------------------------------------------------------
 138   inline
 139   bool
 140   IsEqual(TFloat theLeft, TFloat theRight)
 141   {
 142     static TFloat EPS = 1.0E-3;
 143     if(fabs(theLeft) + fabs(theRight) > EPS)
 144       return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
 145     return true;
 146   }
 147
 148   //---------------------------------------------------------------
 149   class TShapeFun::TFun
 150   {
 151     TFloatVector myFun;
 152     TInt myNbRef;
 153
 154   public:
 155
 156     void
 157     Init(TInt theNbGauss,
 158          TInt theNbRef)
 159     {
 160       myFun.resize(theNbGauss*theNbRef);
 161       myNbRef = theNbRef;
 162     }
 163
 164     TCFloatVecSlice
 165     GetFunSlice(TInt theGaussId) const
 166     {
 167       return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
 168     }
 169
 170     TFloatVecSlice
 171     GetFunSlice(TInt theGaussId)
 172     {
 173       return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
 174     }
 175   };
 176
 177   //---------------------------------------------------------------
 178
 179   TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
 180     myRefCoord(theNbRef*theDim),
 181     myDim(theDim),
 182     myNbRef(theNbRef)
 183   {}
 184
 185   TCCoordSlice
 186   TShapeFun::GetCoord(TInt theRefId) const
 187   {
 188     return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
 189   }
 190
 191   TCoordSlice
 192   TShapeFun::GetCoord(TInt theRefId)
 193   {
 194     return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
 195   }
 196
 197   void
 198   TShapeFun::GetFun(const TCCoordSliceArr& theRef,
 199                     const TCCoordSliceArr& theGauss,
 200                     TFun& theFun) const
 201   {
 202     TInt aNbRef = theRef.size();
 203     TInt aNbGauss = theGauss.size();
 204     theFun.Init(aNbGauss,aNbRef);
 205   }
 206
 207   bool
 208   TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
 209   {
 210     TInt aNbRef = theRefCoord.size();
 211     TInt aNbRef2 = GetNbRef();
 212     INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
 213             "- aNbRef("<<aNbRef<<")"<<
 214             "; aNbRef2("<<aNbRef2<<")\n");
 215     bool anIsSatisfy = (aNbRef == aNbRef2);
 216     if(anIsSatisfy){
 217       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 218         const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
 219         TCCoordSlice aCoord = GetCoord(aRefId);
 220         TInt aDim = aCoord.size();
 221         bool anIsEqual = false;
 222         for(TInt anId = 0; anId < aDim; anId++){
 223           anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
 224           if(!anIsEqual){
 225             anIsSatisfy = false;
 226             break;
 227           }
 228         }
 229         if(!anIsEqual){
 230 #ifdef _DEBUG_
 231           TCCoordSlice aCoord = GetCoord(aRefId);
 232           INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
 233           TInt aDim = aCoord.size();
 234           for(TInt anId = 0; anId < aDim; anId++)
 235             ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
 236           const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
 237           ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
 238           for(TInt anId = 0; anId < aDim; anId++)
 239             ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
 240           ADDMSG(MYDEBUG,"}\n");
 241 #endif
 242 #ifndef _DEBUG_
 243           BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
 244           return anIsSatisfy;
 245 #endif
 246         }
 247       }
 248     }
 249
 250     BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
 251     return anIsSatisfy;
 252   }
 253
 254   bool
 255   TShapeFun::Eval(const TCellInfo&       theCellInfo,
 256                   const TNodeInfo&       theNodeInfo,
 257                   const TElemNum&        theElemNum,
 258                   const TCCoordSliceArr& theRef,
 259                   const TCCoordSliceArr& theGauss,
 260                   TGaussCoord&           theGaussCoord,
 261                   EModeSwitch            theMode)
 262   {
 263     INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);
 264
 265     if(IsSatisfy(theRef)){
 266       const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
 267       TInt aDim = aMeshInfo->GetDim();
 268       TInt aNbGauss = theGauss.size();
 269
 270       bool anIsSubMesh = !theElemNum.empty();
 271       TInt aNbElem;
 272       if(anIsSubMesh)
 273         aNbElem = theElemNum.size();
 274       else
 275         aNbElem = theCellInfo.GetNbElem();
 276
 277       theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
 278
 279       TFun aFun;
 280       InitFun(theRef,theGauss,aFun);
 281       TInt aConnDim = theCellInfo.GetConnDim();
 282
 283       INITMSG(MYDEBUG,"aDim = "<<aDim<<
 284               "; aNbGauss = "<<aNbGauss<<
 285               "; aNbElem = "<<aNbElem<<
 286               "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
 287               std::endl);
 288
 289       for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
 290         TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
 291         TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
 292         TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
 293
 294         for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 295           TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
 296           TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
 297
 298           for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
 299             TInt aNodeId = aConnSlice[aConnId] - 1;
 300             TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
 301
 302             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
 303               aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
 304             }
 305           }
 306         }
 307       }
 308
 309 #ifdef _DEBUG_
 310       {
 311         INITMSG(MYVALUEDEBUG,"theGauss: ");
 312         for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 313           TCCoordSlice aCoordSlice = theGauss[aGaussId];
 314           ADDMSG(MYVALUEDEBUG,"{");
 315           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
 316             ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
 317           }
 318           ADDMSG(MYVALUEDEBUG,"} ");
 319         }
 320         ADDMSG(MYVALUEDEBUG,std::endl);
 321       }
 322       for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
 323         TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
 324         INITMSG(MYVALUEDEBUG,"");
 325         for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 326           TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
 327           ADDMSG(MYVALUEDEBUG,"{");
 328           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
 329             ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
 330           }
 331           ADDMSG(MYVALUEDEBUG,"} ");
 332         }
 333         ADDMSG(MYVALUEDEBUG,std::endl);
 334       }
 335 #endif
 336       return true;
 337     }
 338
 339     return false;
 340   }
 341
 342   //---------------------------------------------------------------
 343   TSeg2a::TSeg2a():TShapeFun(1,2)
 344   {
 345     TInt aNbRef = GetNbRef();
 346     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 347       TCoordSlice aCoord = GetCoord(aRefId);
 348       switch(aRefId){
 349       case  0: aCoord[0] = -1.0; break;
 350       case  1: aCoord[0] =  1.0; break;
 351       }
 352     }
 353   }
 354
 355   void
 356   TSeg2a::InitFun(const TCCoordSliceArr& theRef,
 357                   const TCCoordSliceArr& theGauss,
 358                   TFun& theFun) const
 359   {
 360     GetFun(theRef,theGauss,theFun);
 361
 362     TInt aNbGauss = theGauss.size();
 363     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 364       const TCCoordSlice& aCoord = theGauss[aGaussId];
 365       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 366
 367       aSlice[0] = 0.5*(1.0 - aCoord[0]);
 368       aSlice[1] = 0.5*(1.0 + aCoord[0]);
 369     }
 370   }
 371
 372   //---------------------------------------------------------------
 373   TSeg3a::TSeg3a():TShapeFun(1,3)
 374   {
 375     TInt aNbRef = GetNbRef();
 376     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 377       TCoordSlice aCoord = GetCoord(aRefId);
 378       switch(aRefId){
 379       case  0: aCoord[0] = -1.0; break;
 380       case  1: aCoord[0] =  1.0; break;
 381       case  2: aCoord[0] =  0.0; break;
 382       }
 383     }
 384   }
 385
 386   void
 387   TSeg3a::InitFun(const TCCoordSliceArr& theRef,
 388                   const TCCoordSliceArr& theGauss,
 389                   TFun& theFun) const
 390   {
 391     GetFun(theRef,theGauss,theFun);
 392
 393     TInt aNbGauss = theGauss.size();
 394     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 395       const TCCoordSlice& aCoord = theGauss[aGaussId];
 396       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 397
 398       aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
 399       aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
 400       aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
 401     }
 402   }
 403
 404   //---------------------------------------------------------------
 405   TTria3a::TTria3a():
 406     TShapeFun(2,3)
 407   {
 408     TInt aNbRef = GetNbRef();
 409     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 410       TCoordSlice aCoord = GetCoord(aRefId);
 411       switch(aRefId){
 412       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 413       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 414       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 415       }
 416     }
 417   }
 418
 419   void
 420   TTria3a::InitFun(const TCCoordSliceArr& theRef,
 421                    const TCCoordSliceArr& theGauss,
 422                    TFun& theFun) const
 423   {
 424     GetFun(theRef,theGauss,theFun);
 425
 426     TInt aNbGauss = theGauss.size();
 427     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 428       const TCCoordSlice& aCoord = theGauss[aGaussId];
 429       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 430
 431       aSlice[0] = 0.5*(1.0 + aCoord[1]);
 432       aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
 433       aSlice[2] = 0.5*(1.0 + aCoord[0]);
 434     }
 435   }
 436
 437   //---------------------------------------------------------------
 438   TTria6a::TTria6a():TShapeFun(2,6)
 439   {
 440     TInt aNbRef = GetNbRef();
 441     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 442       TCoordSlice aCoord = GetCoord(aRefId);
 443       switch(aRefId){
 444       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 445       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 446       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 447
 448       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 449       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 450       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 451       }
 452     }
 453   }
 454
 455   void
 456   TTria6a::InitFun(const TCCoordSliceArr& theRef,
 457                    const TCCoordSliceArr& theGauss,
 458                    TFun& theFun) const
 459   {
 460     GetFun(theRef,theGauss,theFun);
 461
 462     TInt aNbGauss = theGauss.size();
 463     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 464       const TCCoordSlice& aCoord = theGauss[aGaussId];
 465       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 466
 467       aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
 468       aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
 469       aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
 470
 471       aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
 472       aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
 473       aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
 474     }
 475   }
 476
 477   //---------------------------------------------------------------
 478   TTria3b::TTria3b():
 479     TShapeFun(2,3)
 480   {
 481     TInt aNbRef = GetNbRef();
 482     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 483       TCoordSlice aCoord = GetCoord(aRefId);
 484       switch(aRefId){
 485       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 486       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 487       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 488       }
 489     }
 490   }
 491
 492   void
 493   TTria3b::InitFun(const TCCoordSliceArr& theRef,
 494                    const TCCoordSliceArr& theGauss,
 495                    TFun& theFun) const
 496   {
 497     GetFun(theRef,theGauss,theFun);
 498
 499     TInt aNbGauss = theGauss.size();
 500     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 501       const TCCoordSlice& aCoord = theGauss[aGaussId];
 502       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 503
 504       aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
 505       aSlice[1] = aCoord[0];
 506       aSlice[2] = aCoord[1];
 507     }
 508   }
 509
 510   //---------------------------------------------------------------
 511   TTria6b::TTria6b():
 512     TShapeFun(2,6)
 513   {
 514     TInt aNbRef = GetNbRef();
 515     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 516       TCoordSlice aCoord = GetCoord(aRefId);
 517       switch(aRefId){
 518       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 519       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 520       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 521
 522       case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
 523       case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
 524       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
 525       }
 526     }
 527   }
 528
 529   void
 530   TTria6b::InitFun(const TCCoordSliceArr& theRef,
 531                    const TCCoordSliceArr& theGauss,
 532                    TFun& theFun) const
 533   {
 534     GetFun(theRef,theGauss,theFun);
 535
 536     TInt aNbGauss = theGauss.size();
 537     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 538       const TCCoordSlice& aCoord = theGauss[aGaussId];
 539       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 540
 541       aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
 542       aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
 543       aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
 544
 545       aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
 546       aSlice[4] = 4.0*aCoord[0]*aCoord[1];
 547       aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
 548     }
 549   }
 550
 551   //---------------------------------------------------------------
 552   TQuad4a::TQuad4a():
 553     TShapeFun(2,4)
 554   {
 555     TInt aNbRef = GetNbRef();
 556     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 557       TCoordSlice aCoord = GetCoord(aRefId);
 558       switch(aRefId){
 559       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 560       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 561       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 562       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 563       }
 564     }
 565   }
 566
 567   void
 568   TQuad4a::InitFun(const TCCoordSliceArr& theRef,
 569                    const TCCoordSliceArr& theGauss,
 570                    TFun& theFun) const
 571   {
 572     GetFun(theRef,theGauss,theFun);
 573
 574     TInt aNbGauss = theGauss.size();
 575     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 576       const TCCoordSlice& aCoord = theGauss[aGaussId];
 577       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 578
 579       aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
 580       aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
 581       aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
 582       aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
 583     }
 584   }
 585
 586   //---------------------------------------------------------------
 587   TQuad8a::TQuad8a():
 588     TShapeFun(2,8)
 589   {
 590     TInt aNbRef = GetNbRef();
 591     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 592       TCoordSlice aCoord = GetCoord(aRefId);
 593       switch(aRefId){
 594       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 595       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 596       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 597       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 598
 599       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 600       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 601       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 602       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 603       }
 604     }
 605   }
 606
 607   void
 608   TQuad8a::InitFun(const TCCoordSliceArr& theRef,
 609                    const TCCoordSliceArr& theGauss,
 610                    TFun& theFun) const
 611   {
 612     GetFun(theRef,theGauss,theFun);
 613
 614     TInt aNbGauss = theGauss.size();
 615     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 616       const TCCoordSlice& aCoord = theGauss[aGaussId];
 617       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 618
 619       aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
 620       aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
 621       aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
 622       aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);
 623
 624       aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
 625       aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
 626       aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
 627       aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
 628     }
 629   }
 630
 631   //---------------------------------------------------------------
 632   TQuad9a::TQuad9a():
 633     TShapeFun(2,9)
 634   {
 635     TInt aNbRef = GetNbRef();
 636     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 637       TCoordSlice aCoord = GetCoord(aRefId);
 638       switch(aRefId){
 639       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 640       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 641       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 642       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 643
 644       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 645       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 646       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 647       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 648
 649       case  8: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 650       }
 651     }
 652   }
 653
 654   void
 655   TQuad9a::InitFun(const TCCoordSliceArr& theRef,
 656                    const TCCoordSliceArr& theGauss,
 657                    TFun& theFun) const
 658   {
 659     GetFun(theRef,theGauss,theFun);
 660
 661     TInt aNbGauss = theGauss.size();
 662     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 663       const TCCoordSlice& aCoord = theGauss[aGaussId];
 664       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 665
 666       aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
 667       aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
 668       aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
 669       aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
 670
 671       aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
 672       aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
 673       aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
 674       aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
 675
 676       aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
 677     }
 678   }
 679
 680   //---------------------------------------------------------------
 681   TQuad4b::TQuad4b():
 682     TShapeFun(2,4)
 683   {
 684     TInt aNbRef = GetNbRef();
 685     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 686       TCoordSlice aCoord = GetCoord(aRefId);
 687       switch(aRefId){
 688       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 689       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 690       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 691       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 692       }
 693     }
 694   }
 695
 696   void
 697   TQuad4b::InitFun(const TCCoordSliceArr& theRef,
 698                    const TCCoordSliceArr& theGauss,
 699                    TFun& theFun) const
 700   {
 701     GetFun(theRef,theGauss,theFun);
 702
 703     TInt aNbGauss = theGauss.size();
 704     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 705       const TCCoordSlice& aCoord = theGauss[aGaussId];
 706       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 707
 708       aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
 709       aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
 710       aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
 711       aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
 712     }
 713   }
 714
 715   //---------------------------------------------------------------
 716   TQuad8b::TQuad8b():
 717     TShapeFun(2,8)
 718   {
 719     TInt aNbRef = GetNbRef();
 720     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 721       TCoordSlice aCoord = GetCoord(aRefId);
 722       switch(aRefId){
 723       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 724       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 725       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 726       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 727
 728       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 729       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 730       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 731       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 732       }
 733     }
 734   }
 735
 736   void
 737   TQuad8b::InitFun(const TCCoordSliceArr& theRef,
 738                    const TCCoordSliceArr& theGauss,
 739                    TFun& theFun) const
 740   {
 741     GetFun(theRef,theGauss,theFun);
 742
 743     TInt aNbGauss = theGauss.size();
 744     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 745       const TCCoordSlice& aCoord = theGauss[aGaussId];
 746       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 747
 748       aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
 749       aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
 750       aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
 751       aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);
 752
 753       aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
 754       aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
 755       aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
 756       aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);
 757
 758       //aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
 759       //aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
 760       //aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
 761       //aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
 762     }
 763   }
 764
 765   //---------------------------------------------------------------
 766   TQuad9b::TQuad9b():
 767     TShapeFun(2,9)
 768   {
 769     TInt aNbRef = GetNbRef();
 770     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 771       TCoordSlice aCoord = GetCoord(aRefId);
 772       switch(aRefId){
 773       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 774       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 775       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 776       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 777
 778       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 779       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 780       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 781       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
 782
 783       case  8: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 784       }
 785     }
 786   }
 787
 788   void
 789   TQuad9b::InitFun(const TCCoordSliceArr& theRef,
 790                    const TCCoordSliceArr& theGauss,
 791                    TFun& theFun) const
 792   {
 793     GetFun(theRef,theGauss,theFun);
 794
 795     TInt aNbGauss = theGauss.size();
 796     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 797       const TCCoordSlice& aCoord = theGauss[aGaussId];
 798       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 799
 800       aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
 801       aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
 802       aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
 803       aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
 804
 805       aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
 806       aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
 807       aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
 808       aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
 809
 810       aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
 811     }
 812   }
 813
 814   //---------------------------------------------------------------
 815   TTetra4a::TTetra4a():
 816     TShapeFun(3,4)
 817   {
 818     TInt aNbRef = GetNbRef();
 819     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 820       TCoordSlice aCoord = GetCoord(aRefId);
 821       switch(aRefId){
 822       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 823       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 824       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 825       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 826       }
 827     }
 828   }
 829
 830   void
 831   TTetra4a::InitFun(const TCCoordSliceArr& theRef,
 832                     const TCCoordSliceArr& theGauss,
 833                     TFun& theFun) const
 834   {
 835     GetFun(theRef,theGauss,theFun);
 836
 837     TInt aNbGauss = theGauss.size();
 838     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 839       const TCCoordSlice& aCoord = theGauss[aGaussId];
 840       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 841
 842       aSlice[0] = aCoord[1];
 843       aSlice[1] = aCoord[2];
 844       aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
 845       aSlice[3] = aCoord[0];
 846     }
 847   }
 848
 849   //---------------------------------------------------------------
 850   TTetra10a::TTetra10a():
 851     TShapeFun(3,10)
 852   {
 853     TInt aNbRef = GetNbRef();
 854     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 855       TCoordSlice aCoord = GetCoord(aRefId);
 856       switch(aRefId){
 857       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 858       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 859       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 860       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 861
 862       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 863       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 864       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 865
 866       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 867       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 868       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 869       }
 870     }
 871   }
 872
 873   void
 874   TTetra10a::InitFun(const TCCoordSliceArr& theRef,
 875                      const TCCoordSliceArr& theGauss,
 876                      TFun& theFun) const
 877   {
 878     GetFun(theRef,theGauss,theFun);
 879
 880     TInt aNbGauss = theGauss.size();
 881     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 882       const TCCoordSlice& aCoord = theGauss[aGaussId];
 883       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 884
 885       aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
 886       aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
 887       aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
 888       aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
 889
 890       aSlice[4] = 4.0*aCoord[1]*aCoord[2];
 891       aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 892       aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 893
 894       aSlice[7] = 4.0*aCoord[0]*aCoord[1];
 895       aSlice[8] = 4.0*aCoord[0]*aCoord[2];
 896       aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 897     }
 898   }
 899
 900   //---------------------------------------------------------------
 901
 902   TTetra4b::TTetra4b():
 903     TShapeFun(3,4)
 904   {
 905     TInt aNbRef = GetNbRef();
 906     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 907       TCoordSlice aCoord = GetCoord(aRefId);
 908       switch(aRefId){
 909       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 910       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 911       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 912       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 913       }
 914     }
 915   }
 916
 917   void
 918   TTetra4b::InitFun(const TCCoordSliceArr& theRef,
 919                     const TCCoordSliceArr& theGauss,
 920                     TFun& theFun) const
 921   {
 922     GetFun(theRef,theGauss,theFun);
 923
 924     TInt aNbGauss = theGauss.size();
 925     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 926       const TCCoordSlice& aCoord = theGauss[aGaussId];
 927       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 928
 929       aSlice[0] = aCoord[1];
 930       aSlice[2] = aCoord[2];
 931       aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
 932       aSlice[3] = aCoord[0];
 933     }
 934   }
 935
 936   //---------------------------------------------------------------
 937   TTetra10b::TTetra10b():
 938     TShapeFun(3,10)
 939   {
 940     TInt aNbRef = GetNbRef();
 941     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 942       TCoordSlice aCoord = GetCoord(aRefId);
 943       switch(aRefId){
 944       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 945       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 946       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 947       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 948
 949       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 950       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 951       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 952
 953       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 954       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 955       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 956       }
 957     }
 958   }
 959
 960   void
 961   TTetra10b::InitFun(const TCCoordSliceArr& theRef,
 962                      const TCCoordSliceArr& theGauss,
 963                      TFun& theFun) const
 964   {
 965     GetFun(theRef,theGauss,theFun);
 966
 967     TInt aNbGauss = theGauss.size();
 968     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 969       const TCCoordSlice& aCoord = theGauss[aGaussId];
 970       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 971
 972       aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
 973       aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
 974       aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
 975       aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
 976
 977       aSlice[6] = 4.0*aCoord[1]*aCoord[2];
 978       aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 979       aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 980
 981       aSlice[7] = 4.0*aCoord[0]*aCoord[1];
 982       aSlice[9] = 4.0*aCoord[0]*aCoord[2];
 983       aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 984     }
 985   }
 986
 987   //---------------------------------------------------------------
 988   THexa8a::THexa8a():
 989     TShapeFun(3,8)
 990   {
 991     TInt aNbRef = GetNbRef();
 992     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 993       TCoordSlice aCoord = GetCoord(aRefId);
 994       switch(aRefId){
 995       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 996       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 997       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 998       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 999       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1000       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1001       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1002       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1003       }
1004     }
1005   }
1006
1007   void
1008   THexa8a::InitFun(const TCCoordSliceArr& theRef,
1009                    const TCCoordSliceArr& theGauss,
1010                    TFun& theFun) const
1011   {
1012     GetFun(theRef,theGauss,theFun);
1013
1014     TInt aNbGauss = theGauss.size();
1015     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1016       const TCCoordSlice& aCoord = theGauss[aGaussId];
1017       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1018
1019       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1020       aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1021       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1022       aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1023
1024       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1025       aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1026       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1027       aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1028     }
1029   }
1030
1031   //---------------------------------------------------------------
1032   THexa20a::THexa20a(TInt theDim, TInt theNbRef):
1033     TShapeFun(theDim,theNbRef)
1034   {
1035     TInt aNbRef = myRefCoord.size();
1036     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1037       TCoordSlice aCoord = GetCoord(aRefId);
1038       switch(aRefId){
1039       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1040       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1041       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1042       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1043       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1044       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1045       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1046       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1047
1048       case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1049       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
1050       case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1051       case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
1052       case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1053       case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1054       case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1055       case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1056       case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1057       case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1058       case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1059       case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1060       }
1061     }
1062   }
1063
1064   void
1065   THexa20a::InitFun(const TCCoordSliceArr& theRef,
1066                     const TCCoordSliceArr& theGauss,
1067                     TFun& theFun) const
1068   {
1069     GetFun(theRef,theGauss,theFun);
1070
1071     TInt aNbGauss = theGauss.size();
1072     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1073       const TCCoordSlice& aCoord = theGauss[aGaussId];
1074       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1075
1076       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1077         (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1078       aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1079         (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1080       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1081         (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1082       aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1083         (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1084       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1085         (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1086       aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1087         (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1088       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1089         (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1090       aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1091         (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1092
1093       aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1094       aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1095       aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1096       aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1097       aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1098       aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1099       aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1100       aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1101       aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1102       aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1103       aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1104       aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1105     }
1106   }
1107
1108   //---------------------------------------------------------------
1109   THexa27a::THexa27a():
1110     THexa20a(3,27)
1111   {
1112     TInt aNbRef = myRefCoord.size();
1113     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1114       TCoordSlice aCoord = GetCoord(aRefId);
1115       switch(aRefId){
1116       case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
1117       case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1118       case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1119       case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1120       case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1121       case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1122       case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1123       }
1124     }
1125   }
1126
1127   void
1128   THexa27a::InitFun(const TCCoordSliceArr& theRef,
1129                     const TCCoordSliceArr& theGauss,
1130                     TFun& theFun) const
1131   {
1132     GetFun(theRef,theGauss,theFun);
1133
1134     TInt aNbGauss = theGauss.size();
1135     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1136       const TCCoordSlice& aCoord = theGauss[aGaussId];
1137       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1138
1139       aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1140       aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1141       aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1142       aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1143       aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1144       aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1145       aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1146       aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1147
1148       aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1149       aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1150       aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1151       aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1152       aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1153       aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1154       aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1155       aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1156       aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1157       aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1158       aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1159       aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1160       aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1161       aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1162       aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1163       aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1164       aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1165       aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1166       aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
1167     }
1168   }
1169
1170   //---------------------------------------------------------------
1171   THexa8b::THexa8b():
1172     TShapeFun(3,8)
1173   {
1174     TInt aNbRef = GetNbRef();
1175     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1176       TCoordSlice aCoord = GetCoord(aRefId);
1177       switch(aRefId){
1178       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1179       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1180       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1181       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1182       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1183       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1184       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1185       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1186       }
1187     }
1188   }
1189
1190   void
1191   THexa8b::InitFun(const TCCoordSliceArr& theRef,
1192                    const TCCoordSliceArr& theGauss,
1193                    TFun& theFun) const
1194   {
1195     GetFun(theRef,theGauss,theFun);
1196
1197     TInt aNbGauss = theGauss.size();
1198     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1199       const TCCoordSlice& aCoord = theGauss[aGaussId];
1200       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1201
1202       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1203       aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1204       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1205       aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1206
1207       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1208       aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1209       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1210       aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1211     }
1212   }
1213
1214   //---------------------------------------------------------------
1215   THexa20b::THexa20b(TInt theDim, TInt theNbRef):
1216     TShapeFun(theDim,theNbRef)
1217   {
1218     TInt aNbRef = myRefCoord.size();
1219     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1220       TCoordSlice aCoord = GetCoord(aRefId);
1221       switch(aRefId){
1222       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1223       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1224       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1225       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1226       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1227       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1228       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1229       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1230
1231       case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
1232       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
1233       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
1234       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
1235       case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1236       case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1237       case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1238       case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1239       case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
1240       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1241       case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
1242       case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1243       }
1244     }
1245   }
1246
1247   void
1248   THexa20b::InitFun(const TCCoordSliceArr& theRef,
1249                     const TCCoordSliceArr& theGauss,
1250                     TFun& theFun) const
1251   {
1252     GetFun(theRef,theGauss,theFun);
1253
1254     TInt aNbGauss = theGauss.size();
1255     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1256       const TCCoordSlice& aCoord = theGauss[aGaussId];
1257       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1258
1259       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1260         (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1261       aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1262         (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1263       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1264         (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1265       aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1266         (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1267       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1268         (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1269       aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1270         (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1271       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1272         (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1273       aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1274         (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1275
1276       aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1277       aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1278       aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1279       aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1280       aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1281       aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1282       aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1283       aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1284       aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1285       aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1286       aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1287       aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1288     }
1289   }
1290
1291   //---------------------------------------------------------------
1292   TPenta6a::TPenta6a():
1293     TShapeFun(3,6)
1294   {
1295     TInt aNbRef = myRefCoord.size();
1296     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1297       TCoordSlice aCoord = GetCoord(aRefId);
1298       switch(aRefId){
1299       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1300       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
1301       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1302       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1303       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1304       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1305       }
1306     }
1307   }
1308
1309   void
1310   TPenta6a::InitFun(const TCCoordSliceArr& theRef,
1311                     const TCCoordSliceArr& theGauss,
1312                     TFun& theFun) const
1313   {
1314     GetFun(theRef,theGauss,theFun);
1315
1316     TInt aNbGauss = theGauss.size();
1317     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1318       const TCCoordSlice& aCoord = theGauss[aGaussId];
1319       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1320
1321       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1322       aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1323       aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1324
1325       aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1326       aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1327       aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1328     }
1329   }
1330
1331   //---------------------------------------------------------------
1332   TPenta6b::TPenta6b():
1333     TShapeFun(3,6)
1334   {
1335     TInt aNbRef = myRefCoord.size();
1336     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1337       TCoordSlice aCoord = GetCoord(aRefId);
1338       switch(aRefId){
1339       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1340       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
1341       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1342       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1343       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1344       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1345       }
1346     }
1347   }
1348
1349   void
1350   TPenta6b::InitFun(const TCCoordSliceArr& theRef,
1351                     const TCCoordSliceArr& theGauss,
1352                     TFun& theFun) const
1353   {
1354     GetFun(theRef,theGauss,theFun);
1355
1356     TInt aNbGauss = theGauss.size();
1357     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1358       const TCCoordSlice& aCoord = theGauss[aGaussId];
1359       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1360
1361       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1362       aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1363       aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1364
1365       aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1366       aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1367       aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1368     }
1369   }
1370
1371   //---------------------------------------------------------------
1372   TPenta15a::TPenta15a():
1373     TShapeFun(3,15)
1374   {
1375     TInt aNbRef = myRefCoord.size();
1376     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1377       TCoordSlice aCoord = GetCoord(aRefId);
1378       switch(aRefId){
1379       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1380       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
1381       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1382       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1383       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1384       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1385
1386       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
1387       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1388       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1389       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1390       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1391       case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1392       case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
1393       case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1394       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1395       }
1396     }
1397   }
1398
1399   void
1400   TPenta15a::InitFun(const TCCoordSliceArr& theRef,
1401                      const TCCoordSliceArr& theGauss,
1402                      TFun& theFun) const
1403   {
1404     GetFun(theRef,theGauss,theFun);
1405
1406     TInt aNbGauss = theGauss.size();
1407     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1408       const TCCoordSlice& aCoord = theGauss[aGaussId];
1409       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1410
1411       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1412       aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1413       aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1414
1415       aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1416       aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1417       aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1418
1419       aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1420       aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1421       aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1422
1423       aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1424       aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1425       aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1426
1427       aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1428       aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1429       aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1430     }
1431   }
1432
1433   //---------------------------------------------------------------
1434   TPenta15b::TPenta15b():
1435     TShapeFun(3,15)
1436   {
1437     TInt aNbRef = myRefCoord.size();
1438     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1439       TCoordSlice aCoord = GetCoord(aRefId);
1440       switch(aRefId){
1441       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1442       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
1443       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1444       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1445       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1446       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1447
1448       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
1449       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1450       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1451       case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1452       case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1453       case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1454       case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
1455       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1456       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1457       }
1458     }
1459   }
1460
1461   void
1462   TPenta15b::InitFun(const TCCoordSliceArr& theRef,
1463                      const TCCoordSliceArr& theGauss,
1464                      TFun& theFun) const
1465   {
1466     GetFun(theRef,theGauss,theFun);
1467
1468     TInt aNbGauss = theGauss.size();
1469     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1470       const TCCoordSlice& aCoord = theGauss[aGaussId];
1471       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1472
1473       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1474       aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1475       aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1476
1477       aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1478       aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1479       aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1480
1481       aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1482       aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1483       aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1484
1485       aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1486       aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1487       aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1488
1489       aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1490       aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1491       aSlice[9]  = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1492     }
1493   }
1494
1495   //---------------------------------------------------------------
1496   TPyra5a::TPyra5a():
1497     TShapeFun(3,5)
1498   {
1499     TInt aNbRef = myRefCoord.size();
1500     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1501       TCoordSlice aCoord = GetCoord(aRefId);
1502       switch(aRefId){
1503       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1504       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1505       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1506       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1507       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1508       }
1509     }
1510   }
1511
1512   void
1513   TPyra5a::InitFun(const TCCoordSliceArr& theRef,
1514                    const TCCoordSliceArr& theGauss,
1515                    TFun& theFun) const
1516   {
1517     GetFun(theRef,theGauss,theFun);
1518
1519     TInt aNbGauss = theGauss.size();
1520     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1521       const TCCoordSlice& aCoord = theGauss[aGaussId];
1522       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1523       // APO & RNV:
1524       // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1525       // Seems shape function for ePYRA5 elements is:
1526       // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1527       // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1528       // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1529       // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1530       // w4 = +Z;
1531       aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1532       aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1533       aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1534       aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1535       aSlice[4] = aCoord[2];
1536     }
1537   }
1538
1539   //---------------------------------------------------------------
1540   TPyra5b::TPyra5b():
1541     TShapeFun(3,5)
1542   {
1543     TInt aNbRef = myRefCoord.size();
1544     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1545       TCoordSlice aCoord = GetCoord(aRefId);
1546       switch(aRefId){
1547       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1548       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1549       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1550       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1551       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1552       }
1553     }
1554   }
1555
1556   void
1557   TPyra5b::InitFun(const TCCoordSliceArr& theRef,
1558                    const TCCoordSliceArr& theGauss,
1559                    TFun& theFun) const
1560   {
1561     GetFun(theRef,theGauss,theFun);
1562
1563     TInt aNbGauss = theGauss.size();
1564     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1565       const TCCoordSlice& aCoord = theGauss[aGaussId];
1566       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1567       // APO & RNV:
1568       // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1569       // Seems shape function for ePYRA5 elements is:
1570       // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1571       // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1572       // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1573       // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1574       // w4 = +Z;
1575       aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1576       aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1577       aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1578       aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1579       aSlice[4] = aCoord[2];
1580     }
1581   }
1582
1583   //---------------------------------------------------------------
1584   TPyra13a::TPyra13a():
1585     TShapeFun(3,13)
1586   {
1587     TInt aNbRef = myRefCoord.size();
1588     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1589       TCoordSlice aCoord = GetCoord(aRefId);
1590       switch(aRefId){
1591       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1592       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1593       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1594       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1595       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1596
1597       case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1598       case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1599       case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
1600       case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
1601       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1602       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
1603       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1604       case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
1605       }
1606     }
1607   }
1608
1609   void
1610   TPyra13a::InitFun(const TCCoordSliceArr& theRef,
1611                     const TCCoordSliceArr& theGauss,
1612                     TFun& theFun) const
1613   {
1614     GetFun(theRef,theGauss,theFun);
1615
1616     TInt aNbGauss = theGauss.size();
1617     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1618       const TCCoordSlice& aCoord = theGauss[aGaussId];
1619       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1620
1621       aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1622         (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1623       aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1624         (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1625       aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1626         (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1627       aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1628         (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1629
1630       aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1631
1632       aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1633         (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1634       aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1635         (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1636       aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1637         (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1638       aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1639         (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1640
1641       aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1642         (1.0 - aCoord[2]);
1643       aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1644         (1.0 - aCoord[2]);
1645       aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1646         (1.0 - aCoord[2]);
1647       aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1648         (1.0 - aCoord[2]);
1649     }
1650   }
1651
1652   //---------------------------------------------------------------
1653   TPyra13b::TPyra13b():
1654     TShapeFun(3,13)
1655   {
1656     TInt aNbRef = myRefCoord.size();
1657     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1658       TCoordSlice aCoord = GetCoord(aRefId);
1659       switch(aRefId){
1660       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1661       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
1662       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
1663       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
1664       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
1665
1666       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1667       case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
1668       case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
1669       case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
1670       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1671       case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
1672       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
1673       case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
1674       }
1675     }
1676   }
1677
1678   void
1679   TPyra13b::InitFun(const TCCoordSliceArr& theRef,
1680                     const TCCoordSliceArr& theGauss,
1681                     TFun& theFun) const
1682   {
1683     GetFun(theRef,theGauss,theFun);
1684
1685     TInt aNbGauss = theGauss.size();
1686     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1687       const TCCoordSlice& aCoord = theGauss[aGaussId];
1688       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1689
1690       aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1691         (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1692       aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1693         (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1694       aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1695         (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1696       aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1697         (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1698
1699       aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1700
1701       aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1702         (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1703       aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1704         (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1705       aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1706         (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1707       aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1708         (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1709
1710       aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1711         (1.0 - aCoord[2]);
1712       aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1713         (1.0 - aCoord[2]);
1714       aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1715         (1.0 - aCoord[2]);
1716       aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1717         (1.0 - aCoord[2]);
1718     }
1719   }
1720
1721   //---------------------------------------------------------------
1722   bool
1723   GetGaussCoord3D(const TGaussInfo& theGaussInfo,
1724                   const TCellInfo& theCellInfo,
1725                   const TNodeInfo& theNodeInfo,
1726                   TGaussCoord& theGaussCoord,
1727                   const TElemNum& theElemNum,
1728                   EModeSwitch theMode)
1729   {
1730     INITMSG(MYDEBUG,"GetGaussCoord3D\n");
1731
1732     if(theGaussInfo.myGeom == theCellInfo.myGeom){
1733       EGeometrieElement aGeom = theGaussInfo.myGeom;
1734
1735       TInt aNbRef = theGaussInfo.GetNbRef();
1736       TCCoordSliceArr aRefSlice(aNbRef);
1737       for(TInt anId = 0; anId < aNbRef; anId++)
1738         aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
1739
1740       TInt aNbGauss = theGaussInfo.GetNbGauss();
1741       TCCoordSliceArr aGaussSlice(aNbGauss);
1742       for(TInt anId = 0; anId < aNbGauss; anId++)
1743         aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
1744
1745       switch(aGeom){
1746       case eSEG2: {
1747         INITMSG(MYDEBUG,"eSEG2"<<std::endl);
1748
1749         if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1750           return true;
1751
1752         break;
1753       }
1754       case eSEG3: {
1755         INITMSG(MYDEBUG,"eSEG3"<<std::endl);
1756
1757         if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1758           return true;
1759
1760         break;
1761       }
1762       case eTRIA3: {
1763         INITMSG(MYDEBUG,"eTRIA3"<<std::endl);
1764
1765         if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1766           return true;
1767
1768         if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1769           return true;
1770
1771         break;
1772       }
1773       case eTRIA6: {
1774         INITMSG(MYDEBUG,"eTRIA6"<<std::endl);
1775
1776         if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1777           return true;
1778
1779         if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1780           return true;
1781
1782         break;
1783       }
1784       case eQUAD4: {
1785         INITMSG(MYDEBUG,"eQUAD4"<<std::endl);
1786
1787         if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1788           return true;
1789
1790         if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1791           return true;
1792
1793         break;
1794       }
1795       case eQUAD8: {
1796         INITMSG(MYDEBUG,"eQUAD8"<<std::endl);
1797
1798         if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1799           return true;
1800
1801         if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1802           return true;
1803
1804         break;
1805       }
1806       case eQUAD9: {
1807         INITMSG(MYDEBUG,"eQUAD9"<<std::endl);
1808
1809         if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1810           return true;
1811
1812         if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1813           return true;
1814
1815         break;
1816       }
1817       case eTETRA4: {
1818         INITMSG(MYDEBUG,"eTETRA4"<<std::endl);
1819
1820         if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1821           return true;
1822
1823         if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1824           return true;
1825
1826         break;
1827       }
1828       case ePYRA5: {
1829         INITMSG(MYDEBUG,"ePYRA5"<<std::endl);
1830
1831         if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1832           return true;
1833
1834         if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1835           return true;
1836
1837         break;
1838       }
1839       case ePENTA6: {
1840         INITMSG(MYDEBUG,"ePENTA6"<<std::endl);
1841
1842         if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1843           return true;
1844
1845         if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1846           return true;
1847
1848         break;
1849       }
1850       case eHEXA8: {
1851         INITMSG(MYDEBUG,"eHEXA8"<<std::endl);
1852
1853         if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1854           return true;
1855
1856         if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1857           return true;
1858
1859         break;
1860       }
1861       case eTETRA10: {
1862         INITMSG(MYDEBUG,"eTETRA10"<<std::endl);
1863
1864         if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1865           return true;
1866
1867         if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1868           return true;
1869
1870         break;
1871       }
1872       case ePYRA13: {
1873         INITMSG(MYDEBUG,"ePYRA13"<<std::endl);
1874
1875         if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1876           return true;
1877
1878         if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1879           return true;
1880
1881         break;
1882       }
1883       case ePENTA15: {
1884         INITMSG(MYDEBUG,"ePENTA15"<<std::endl);
1885
1886         if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1887           return true;
1888
1889         if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1890           return true;
1891
1892         break;
1893       }
1894       case eHEXA20: {
1895         INITMSG(MYDEBUG,"eHEXA20"<<std::endl);
1896
1897         if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1898           return true;
1899
1900         if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1901           return true;
1902
1903         break;
1904       }
1905       default:
1906         INITMSG(MYDEBUG,"eNONE"<<std::endl);
1907         return false;
1908       }
1909     }
1910
1911     return false;
1912   }
1913
1914   //---------------------------------------------------------------
1915   bool
1916   GetBaryCenter(const TCellInfo& theCellInfo,
1917                 const TNodeInfo& theNodeInfo,
1918                 TGaussCoord& theGaussCoord,
1919                 const TElemNum& theElemNum,
1920                 EModeSwitch theMode)
1921   {
1922     INITMSG(MYDEBUG,"GetBaryCenter\n");
1923     const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
1924     TInt aDim = aMeshInfo->GetDim();
1925     static TInt aNbGauss = 1;
1926
1927     bool anIsSubMesh = !theElemNum.empty();
1928     TInt aNbElem;
1929     if(anIsSubMesh)
1930       aNbElem = theElemNum.size();
1931     else
1932       aNbElem = theCellInfo.GetNbElem();
1933
1934     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
1935
1936     TInt aConnDim = theCellInfo.GetConnDim();
1937
1938     INITMSGA(MYDEBUG,0,
1939              "- aDim = "<<aDim<<
1940              "; aNbGauss = "<<aNbGauss<<
1941              "; aNbElem = "<<aNbElem<<
1942              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
1943              std::endl);
1944
1945     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1946       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
1947       TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
1948       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1949
1950       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1951         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
1952
1953         for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
1954           TInt aNodeId = aConnSlice[aConnId] - 1;
1955           TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
1956
1957           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1958             aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
1959           }
1960         }
1961
1962         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1963           aGaussCoordSlice[aDimId] /= aConnDim;
1964         }
1965       }
1966     }
1967
1968 #ifdef _DEBUG_
1969     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1970       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1971       INITMSG(MYVALUEDEBUG,"");
1972       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1973         TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
1974         ADDMSG(MYVALUEDEBUG,"{");
1975         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1976           ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
1977         }
1978         ADDMSG(MYVALUEDEBUG,"} ");
1979       }
1980       ADDMSG(MYVALUEDEBUG,std::endl);
1981     }
1982 #endif
1983
1984     return true;
1985   }
1986
1987   //---------------------------------------------------------------
1988   bool
1989   GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
1990                 const TNodeInfo& theNodeInfo,
1991                 TGaussCoord& theGaussCoord,
1992                 const TElemNum& theElemNum,
1993                 EModeSwitch theMode)
1994   {
1995     INITMSG(MYDEBUG,"GetBaryCenter\n");
1996     const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
1997     TInt aDim = aMeshInfo->GetDim();
1998     static TInt aNbGauss = 1;
1999
2000     bool anIsSubMesh = !theElemNum.empty();
2001     TInt aNbElem;
2002     if(anIsSubMesh)
2003       aNbElem = theElemNum.size();
2004     else
2005       aNbElem = thePolygoneInfo.GetNbElem();
2006
2007     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2008
2009     INITMSGA(MYDEBUG,0,
2010              "- aDim = "<<aDim<<
2011              "; aNbGauss = "<<aNbGauss<<
2012              "; aNbElem = "<<aNbElem<<
2013              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2014              std::endl);
2015
2016     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2017       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2018
2019       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2020       TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
2021       TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
2022       TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
2023
2024       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2025         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2026
2027         for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2028           TInt aNodeId = aConnSlice[aConnId] - 1;
2029           TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2030
2031           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2032             aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2033           }
2034         }
2035
2036         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2037           aGaussCoordSlice[aDimId] /= aNbNodes;
2038         }
2039       }
2040     }
2041
2042     return true;
2043   }
2044
2045   //---------------------------------------------------------------
2046   bool
2047   GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
2048                 const TNodeInfo& theNodeInfo,
2049                 TGaussCoord& theGaussCoord,
2050                 const TElemNum& theElemNum,
2051                 EModeSwitch theMode)
2052   {
2053     INITMSG(MYDEBUG,"GetBaryCenter\n");
2054     const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
2055     TInt aDim = aMeshInfo->GetDim();
2056     static TInt aNbGauss = 1;
2057
2058     bool anIsSubMesh = !theElemNum.empty();
2059     TInt aNbElem;
2060     if(anIsSubMesh)
2061       aNbElem = theElemNum.size();
2062     else
2063       aNbElem = thePolyedreInfo.GetNbElem();
2064
2065     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2066
2067     INITMSGA(MYDEBUG,0,
2068              "- aDim = "<<aDim<<
2069              "; aNbGauss = "<<aNbGauss<<
2070              "; aNbElem = "<<aNbElem<<
2071              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2072              std::endl);
2073
2074     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2075       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2076
2077       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2078       TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
2079       TInt aNbFaces = aConnSliceArr.size();
2080
2081       TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
2082
2083       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2084         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2085
2086         for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
2087           TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
2088           TInt aNbConn = aConnSlice.size();
2089           for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2090             TInt aNodeId = aConnSlice[aConnId] - 1;
2091             TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2092
2093             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2094               aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2095             }
2096           }
2097         }
2098         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2099           aGaussCoordSlice[aDimId] /= aNbNodes;
2100         }
2101       }
2102     }
2103
2104     return true;
2105   }
2106 }