src/MEDWrapper/Base/MED_GaussUtils.cxx

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