src/MEDWrapper/Base/MED_GaussUtils.cxx

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