src/MEDWrapper/MED_GaussUtils.cxx

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