src/MEDWrapper/Base/MED_GaussUtils.cxx

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