// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
#include "MED_GaussUtils.hxx"
#include "MED_Utilities.hxx"
 
#ifdef _DEBUG_
static int MYDEBUG = 0;
static int MYVALUEDEBUG = 0;
#else
// static int MYDEBUG = 0;
// static int MYVALUEDEBUG = 0;
#endif

//#define _DEBUG_REF_COORDS_

namespace MED
{
  //---------------------------------------------------------------
  TGaussCoord
  ::TGaussCoord():
    TModeSwitchInfo(eFULL_INTERLACE),
    myNbElem(0),
    myNbGauss(0),
    myDim(0),
    myGaussStep(0)
  {
  }

  void
  TGaussCoord
  ::Init(TInt theNbElem,
         TInt theNbGauss,
         TInt theDim,
         EModeSwitch theMode)
  {
    myModeSwitch = theMode;

    myNbElem = theNbElem;
    myNbGauss = theNbGauss;
    myDim = theDim;

    myGaussStep = myNbGauss*myDim;

    myGaussCoord.resize(theNbElem*myGaussStep);
  }


  TInt
  TGaussCoord
  ::GetNbElem() const
  { 
    return myNbElem; 
  }
  
  TInt
  TGaussCoord
  ::GetNbGauss() const
  { 
    return myNbGauss; 
  }
  
  TInt
  TGaussCoord
  ::GetDim() const
  { 
    return myDim; 
  }
  
  unsigned char*
  TGaussCoord
  ::GetValuePtr()
  {
    return (unsigned char*)&(myGaussCoord[0]);
  }


  TCCoordSliceArr 
  TGaussCoord
  ::GetCoordSliceArr(TInt theElemId) const
  {
    TCCoordSliceArr aCoordSliceArr(myNbGauss);
    if(GetModeSwitch() == eFULL_INTERLACE){
      TInt anId = theElemId*myGaussStep;
      for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
        aCoordSliceArr[anGaussId] =
          TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
        anId += myDim;
      }
    }
    else{
      for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
        aCoordSliceArr[anGaussId] =
          TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
      }
    }
    return aCoordSliceArr;
  }


  TCoordSliceArr 
  TGaussCoord
  ::GetCoordSliceArr(TInt theElemId)
  {
    TCoordSliceArr aCoordSliceArr(myNbGauss);
    if(GetModeSwitch() == eFULL_INTERLACE){
      TInt anId = theElemId*myGaussStep;
      for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
        aCoordSliceArr[anGaussId] =
          TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
        anId += myDim;
      }
    }
    else{
      for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
        aCoordSliceArr[anGaussId] =
          TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
      }
    }
    return aCoordSliceArr;
  }


  //---------------------------------------------------------------
  inline
  bool 
  IsEqual(TFloat theLeft, TFloat theRight)
  {
    static TFloat EPS = 1.0E-3;
    if(fabs(theLeft) + fabs(theRight) > EPS)
      return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
    return true;
  }


  //---------------------------------------------------------------
  class TShapeFun::TFun
  {
    TFloatVector myFun;
    TInt myNbRef;

  public:

    void
    Init(TInt theNbGauss,
         TInt theNbRef)
    {
      myFun.resize(theNbGauss*theNbRef);
      myNbRef = theNbRef;
    }

    TCFloatVecSlice 
    GetFunSlice(TInt theGaussId) const
    {
      return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
    }

    TFloatVecSlice
    GetFunSlice(TInt theGaussId)
    {
      return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
    }
  };

  //---------------------------------------------------------------

  TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
    myRefCoord(theNbRef*theDim),
    myDim(theDim),
    myNbRef(theNbRef)
  {}

  TCCoordSlice 
  TShapeFun::GetCoord(TInt theRefId) const
  {
    return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
  }

  TCoordSlice
  TShapeFun::GetCoord(TInt theRefId)
  {
    return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
  }

  void 
  TShapeFun::GetFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    TInt aNbRef = theRef.size();
    TInt aNbGauss = theGauss.size();
    theFun.Init(aNbGauss,aNbRef);
  }

  bool 
  TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
  {
    TInt aNbRef = theRefCoord.size();
    TInt aNbRef2 = GetNbRef();
    INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
            "- aNbRef("<<aNbRef<<")"<<
            "; aNbRef2("<<aNbRef2<<")\n");
    bool anIsSatisfy = (aNbRef == aNbRef2);
    if(anIsSatisfy){
      for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
        const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
        TCCoordSlice aCoord = GetCoord(aRefId);
        TInt aDim = aCoord.size();
        bool anIsEqual = false;
        for(TInt anId = 0; anId < aDim; anId++){
          anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
          if(!anIsEqual){
            anIsSatisfy = false;
            break;
          }
        }
        if(!anIsEqual){
#ifdef _DEBUG_
          TCCoordSlice aCoord = GetCoord(aRefId);
          INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
          TInt aDim = aCoord.size();
          for(TInt anId = 0; anId < aDim; anId++)
            ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
          const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
          ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
          for(TInt anId = 0; anId < aDim; anId++)
            ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
          ADDMSG(MYDEBUG,"}\n");
#endif
#ifndef _DEBUG_
          BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
          return anIsSatisfy;
#endif
        }
      }
    }

    BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
    return anIsSatisfy;
  }

  bool
  TShapeFun::Eval(const TCellInfo&       theCellInfo,
                  const TNodeInfo&       theNodeInfo,
                  const TElemNum&        theElemNum,
                  const TCCoordSliceArr& theRef,
                  const TCCoordSliceArr& theGauss,
                  TGaussCoord&           theGaussCoord,
                  EModeSwitch            theMode)
  {
    INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);

    if(IsSatisfy(theRef)){
      const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
      TInt aDim = aMeshInfo->GetDim();
      TInt aNbGauss = theGauss.size();

      bool anIsSubMesh = !theElemNum.empty();
      TInt aNbElem;
      if(anIsSubMesh)
        aNbElem = theElemNum.size();
      else
        aNbElem = theCellInfo.GetNbElem();

      theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);

      TFun aFun;
      InitFun(theRef,theGauss,aFun);
      TInt aConnDim = theCellInfo.GetConnDim();

      INITMSG(MYDEBUG,"aDim = "<<aDim<<
              "; aNbGauss = "<<aNbGauss<<
              "; aNbElem = "<<aNbElem<<
              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
              std::endl);

      for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
        TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
        TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
        TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);

        for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
          TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
          TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);

          for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
            TInt aNodeId = aConnSlice[aConnId] - 1;      
            TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);

            for(TInt aDimId = 0; aDimId < aDim; aDimId++){
              aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
            }
          }
        }
      }

#ifdef _DEBUG_
      {
        INITMSG(MYVALUEDEBUG,"theGauss: ");
        for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
          TCCoordSlice aCoordSlice = theGauss[aGaussId];
          ADDMSG(MYVALUEDEBUG,"{");
          for(TInt aDimId = 0; aDimId < aDim; aDimId++){
            ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
          }
          ADDMSG(MYVALUEDEBUG,"} ");
        }
        ADDMSG(MYVALUEDEBUG,std::endl);
      }
      for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
        TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
        INITMSG(MYVALUEDEBUG,"");
        for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
          TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
          ADDMSG(MYVALUEDEBUG,"{");
          for(TInt aDimId = 0; aDimId < aDim; aDimId++){
            ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
          }
          ADDMSG(MYVALUEDEBUG,"} ");
        }
        ADDMSG(MYVALUEDEBUG,std::endl);
      }
#endif
      return true;
    }

    return false;
  }


  //---------------------------------------------------------------
  TSeg2a::TSeg2a():TShapeFun(1,2)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0; break;
      case  1: aCoord[0] =  1.0; break;
      }
    }
  }

  void
  TSeg2a::InitFun(const TCCoordSliceArr& theRef,
                  const TCCoordSliceArr& theGauss,
                  TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*(1.0 - aCoord[0]);
      aSlice[1] = 0.5*(1.0 + aCoord[0]);
    }
  }


  //---------------------------------------------------------------
  TSeg3a::TSeg3a():TShapeFun(1,3)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0; break;
      case  1: aCoord[0] =  1.0; break;
      case  2: aCoord[0] =  0.0; break;
      }
    }
  }

  void
  TSeg3a::InitFun(const TCCoordSliceArr& theRef,
                  const TCCoordSliceArr& theGauss,
                  TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
      aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
      aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TTria3a::TTria3a():
    TShapeFun(2,3)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      }
    }
  }

  void
  TTria3a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*(1.0 + aCoord[1]);
      aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
      aSlice[2] = 0.5*(1.0 + aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TTria6a::TTria6a():TShapeFun(2,6)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;

      case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
      case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
      }
    }
  }

  void
  TTria6a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
      aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
      aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];

      aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
      aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
      aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TTria3b::TTria3b():
    TShapeFun(2,3)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
      case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
      }
    }
  }

  void
  TTria3b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
      aSlice[1] = aCoord[0];
      aSlice[2] = aCoord[1];
    }
  }



  //---------------------------------------------------------------
  TTria6b::TTria6b():
    TShapeFun(2,6)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
      case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;

      case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
      case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
      case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
      }
    }
  }

  void
  TTria6b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
      aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
      aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);

      aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
      aSlice[4] = 4.0*aCoord[0]*aCoord[1];
      aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TQuad4a::TQuad4a():
    TShapeFun(2,4)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
      }
    }
  }

  void
  TQuad4a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
      aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
      aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
      aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TQuad8a::TQuad8a():
    TShapeFun(2,8)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;

      case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
      case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
      case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
      case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
      }
    }
  }

  void
  TQuad8a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
      aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
      aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
      aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);

      aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
      aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
      aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
      aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TQuad9a::TQuad9a():
    TShapeFun(2,9)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;

      case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
      case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
      case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
      case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;

      case  8: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
      }
    }
  }

  void
  TQuad9a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
      aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
      aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
      aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);

      aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
      aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
      aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
      aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);

      aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TQuad4b::TQuad4b():
    TShapeFun(2,4)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
      case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
      }
    }
  }

  void
  TQuad4b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
      aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
      aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
      aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TQuad8b::TQuad8b():
    TShapeFun(2,8)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
      case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;

      case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
      case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
      case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
      }
    }
  }

  void
  TQuad8b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
      aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
      aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
      aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);

      aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
      aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
      aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
      aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);

      //aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
      //aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
      //aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
      //aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TQuad9b::TQuad9b():
    TShapeFun(2,9)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
      case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;

      case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
      case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
      case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;

      case  8: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
      }
    }
  }

  void
  TQuad9b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
      aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
      aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
      aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);

      aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
      aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
      aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
      aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);

      aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  TTetra4a::TTetra4a():
    TShapeFun(3,4)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TTetra4a::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = aCoord[1];
      aSlice[1] = aCoord[2];
      aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
      aSlice[3] = aCoord[0];
    }
  }



  //---------------------------------------------------------------
  TTetra10a::TTetra10a():
    TShapeFun(3,10)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;

      case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;

      case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TTetra10a::InitFun(const TCCoordSliceArr& theRef,
                     const TCCoordSliceArr& theGauss,
                     TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
      aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
      aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
      aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);

      aSlice[4] = 4.0*aCoord[1]*aCoord[2];
      aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
      aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);

      aSlice[7] = 4.0*aCoord[0]*aCoord[1];
      aSlice[8] = 4.0*aCoord[0]*aCoord[2];
      aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
    }
  }



  //---------------------------------------------------------------


  TTetra4b::TTetra4b():
    TShapeFun(3,4)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TTetra4b::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = aCoord[1];
      aSlice[2] = aCoord[2];
      aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
      aSlice[3] = aCoord[0];
    }
  }



  //---------------------------------------------------------------
  TTetra10b::TTetra10b():
    TShapeFun(3,10)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;

      case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;

      case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TTetra10b::InitFun(const TCCoordSliceArr& theRef,
                     const TCCoordSliceArr& theGauss,
                     TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
      aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
      aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
      aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);

      aSlice[6] = 4.0*aCoord[1]*aCoord[2];
      aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
      aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);

      aSlice[7] = 4.0*aCoord[0]*aCoord[1];
      aSlice[9] = 4.0*aCoord[0]*aCoord[2];
      aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
    }
  }



  //---------------------------------------------------------------
  THexa8a::THexa8a():
    TShapeFun(3,8)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      }
    }
  }

  void
  THexa8a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
      aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
      aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
      aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);

      aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
      aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
      aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
      aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
    }
  }


  //---------------------------------------------------------------
  THexa20a::THexa20a(TInt theDim, TInt theNbRef):
    TShapeFun(theDim,theNbRef)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;

      case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
      case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
      case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      }
    }
  }

  void
  THexa20a::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
      aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
      aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
      aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
      aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
      aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
      aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
      aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);

      aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
      aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
      aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
      aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
      aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
      aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
      aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
      aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
      aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
      aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
      aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
      aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
    }
  }



  //---------------------------------------------------------------
  THexa27a::THexa27a():
    THexa20a(3,27)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
      case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  THexa27a::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);

      aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
      aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
      aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
      aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
    }
  }



  //---------------------------------------------------------------
  THexa8b::THexa8b():
    TShapeFun(3,8)
  {
    TInt aNbRef = GetNbRef();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      }
    }
  }

  void
  THexa8b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
      aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
      aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
      aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);

      aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
      aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
      aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
      aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
    }
  }



  //---------------------------------------------------------------
  THexa20b::THexa20b(TInt theDim, TInt theNbRef):
    TShapeFun(theDim,theNbRef)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;

      case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
      case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
      case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
      case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
      case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
      case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
      case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      }
    }
  }

  void
  THexa20b::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
      aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
      aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
      aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
        (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
      aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
      aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
      aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
      aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
        (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);

      aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
      aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
      aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
      aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
      aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
      aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
      aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
      aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
      aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
      aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
      aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
      aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
    }
  }



  //---------------------------------------------------------------
  TPenta6a::TPenta6a():
    TShapeFun(3,6)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TPenta6a::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
      aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
      aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);

      aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
      aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
      aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TPenta6b::TPenta6b():
    TShapeFun(3,6)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TPenta6b::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
      aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
      aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);

      aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
      aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
      aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TPenta15a::TPenta15a():
    TShapeFun(3,15)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;

      case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TPenta15a::InitFun(const TCCoordSliceArr& theRef,
                     const TCCoordSliceArr& theGauss,
                     TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
      aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
      aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);

      aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
      aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
      aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);

      aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
      aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
      aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);

      aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
      aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
      aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);

      aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
      aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
      aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TPenta15b::TPenta15b():
    TShapeFun(3,15)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
      case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;

      case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      }
    }
  }

  void
  TPenta15b::InitFun(const TCCoordSliceArr& theRef,
                     const TCCoordSliceArr& theGauss,
                     TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
      aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
      aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);

      aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
      aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
      aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);

      aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
      aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
      aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);

      aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
      aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
      aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);

      aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
      aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
      aSlice[9]  = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
    }
  }



  //---------------------------------------------------------------
  TPyra5a::TPyra5a():
    TShapeFun(3,5)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      }
    }
  }

  void
  TPyra5a::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
      // APO & RNV:
      // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
      // Seems shape function for ePYRA5 elements is:
      // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
      // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
      // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
      // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
      // w4 = +Z;
      aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[4] = aCoord[2];
    }
  }



  //---------------------------------------------------------------
  TPyra5b::TPyra5b():
    TShapeFun(3,5)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){        
      case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
      }
    }
  }

  void
  TPyra5b::InitFun(const TCCoordSliceArr& theRef,
                   const TCCoordSliceArr& theGauss,
                   TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);
    
    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
      // APO & RNV:
      // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
      // Seems shape function for ePYRA5 elements is:
      // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
      // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
      // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
      // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
      // w4 = +Z;
      aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
      aSlice[4] = aCoord[2];
    }
  }
  


  //---------------------------------------------------------------
  TPyra13a::TPyra13a():
    TShapeFun(3,13)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;

      case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
      case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
      case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
      }
    }
  }

  void
  TPyra13a::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
      aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
      aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
      aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);

      aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);

      aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
      aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
      aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
      aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);

      aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
      aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
      aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
      aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
    }
  }



  //---------------------------------------------------------------
  TPyra13b::TPyra13b():
    TShapeFun(3,13)
  {
    TInt aNbRef = myRefCoord.size();
    for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
      TCoordSlice aCoord = GetCoord(aRefId);
      switch(aRefId){
      case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
      case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
      case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
      case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;

      case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
      case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
      case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
      case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
      case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
      case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
      }
    }
  }

  void
  TPyra13b::InitFun(const TCCoordSliceArr& theRef,
                    const TCCoordSliceArr& theGauss,
                    TFun& theFun) const
  {
    GetFun(theRef,theGauss,theFun);

    TInt aNbGauss = theGauss.size();
    for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
      const TCCoordSlice& aCoord = theGauss[aGaussId];
      TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);

      aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
      aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
      aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
      aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);

      aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);

      aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
      aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
        (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
      aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
      aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
        (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);

      aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
      aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
      aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
      aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
        (1.0 - aCoord[2]);
    }
  }



  //---------------------------------------------------------------
  bool
  GetGaussCoord3D(const TGaussInfo& theGaussInfo, 
                  const TCellInfo& theCellInfo,
                  const TNodeInfo& theNodeInfo,
                  TGaussCoord& theGaussCoord,
                  const TElemNum& theElemNum,
                  EModeSwitch theMode)
  {
    INITMSG(MYDEBUG,"GetGaussCoord3D\n");

    if(theGaussInfo.myGeom == theCellInfo.myGeom){
      EGeometrieElement aGeom = theGaussInfo.myGeom;

      TInt aNbRef = theGaussInfo.GetNbRef();
      TCCoordSliceArr aRefSlice(aNbRef);
      for(TInt anId = 0; anId < aNbRef; anId++)
        aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);

      TInt aNbGauss = theGaussInfo.GetNbGauss();
      TCCoordSliceArr aGaussSlice(aNbGauss);
      for(TInt anId = 0; anId < aNbGauss; anId++)
        aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);

      switch(aGeom){
      case eSEG2: {
        INITMSG(MYDEBUG,"eSEG2"<<std::endl);

        if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eSEG3: {
        INITMSG(MYDEBUG,"eSEG3"<<std::endl);

        if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eTRIA3: {
        INITMSG(MYDEBUG,"eTRIA3"<<std::endl);

        if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eTRIA6: {
        INITMSG(MYDEBUG,"eTRIA6"<<std::endl);

        if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eQUAD4: {
        INITMSG(MYDEBUG,"eQUAD4"<<std::endl);

        if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eQUAD8: {
        INITMSG(MYDEBUG,"eQUAD8"<<std::endl);

        if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eQUAD9: {
        INITMSG(MYDEBUG,"eQUAD9"<<std::endl);

        if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eTETRA4: {
        INITMSG(MYDEBUG,"eTETRA4"<<std::endl);

        if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case ePYRA5: {
        INITMSG(MYDEBUG,"ePYRA5"<<std::endl);

        if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case ePENTA6: {
        INITMSG(MYDEBUG,"ePENTA6"<<std::endl);

        if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eHEXA8: {
        INITMSG(MYDEBUG,"eHEXA8"<<std::endl);

        if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eTETRA10: {
        INITMSG(MYDEBUG,"eTETRA10"<<std::endl);

        if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case ePYRA13: {
        INITMSG(MYDEBUG,"ePYRA13"<<std::endl);

        if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case ePENTA15: {
        INITMSG(MYDEBUG,"ePENTA15"<<std::endl);

        if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      case eHEXA20: {
        INITMSG(MYDEBUG,"eHEXA20"<<std::endl);

        if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
          return true;

        break;
      }
      default: 
        INITMSG(MYDEBUG,"eNONE"<<std::endl);
        return false;
      }
    }

    return false;
  }

  //---------------------------------------------------------------
  bool
  GetBaryCenter(const TCellInfo& theCellInfo,
                const TNodeInfo& theNodeInfo,
                TGaussCoord& theGaussCoord,
                const TElemNum& theElemNum,
                EModeSwitch theMode)
  {
    INITMSG(MYDEBUG,"GetBaryCenter\n");
    const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
    TInt aDim = aMeshInfo->GetDim();
    static TInt aNbGauss = 1;

    bool anIsSubMesh = !theElemNum.empty();
    TInt aNbElem;
    if(anIsSubMesh)
      aNbElem = theElemNum.size();
    else
      aNbElem = theCellInfo.GetNbElem();

    theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);

    TInt aConnDim = theCellInfo.GetConnDim();

    INITMSGA(MYDEBUG,0,
             "- aDim = "<<aDim<<
             "; aNbGauss = "<<aNbGauss<<
             "; aNbElem = "<<aNbElem<<
             "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
             std::endl);

    for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
      TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
      TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
      TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);

      for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
        TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];

        for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
          TInt aNodeId = aConnSlice[aConnId] - 1;      
          TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);

          for(TInt aDimId = 0; aDimId < aDim; aDimId++){
            aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
          }
        }

        for(TInt aDimId = 0; aDimId < aDim; aDimId++){
          aGaussCoordSlice[aDimId] /= aConnDim;
        }
      }
    }

#ifdef _DEBUG_
    for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
      TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
      INITMSG(MYVALUEDEBUG,"");
      for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
        TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
        ADDMSG(MYVALUEDEBUG,"{");
        for(TInt aDimId = 0; aDimId < aDim; aDimId++){
          ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
        }
        ADDMSG(MYVALUEDEBUG,"} ");
      }
      ADDMSG(MYVALUEDEBUG,std::endl);
    }
#endif

    return true;
  }


  //---------------------------------------------------------------
  bool
  GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
                const TNodeInfo& theNodeInfo,
                TGaussCoord& theGaussCoord,
                const TElemNum& theElemNum,
                EModeSwitch theMode)
  {
    INITMSG(MYDEBUG,"GetBaryCenter\n");
    const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
    TInt aDim = aMeshInfo->GetDim();
    static TInt aNbGauss = 1;

    bool anIsSubMesh = !theElemNum.empty();
    TInt aNbElem;
    if(anIsSubMesh)
      aNbElem = theElemNum.size();
    else
      aNbElem = thePolygoneInfo.GetNbElem();

    theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);

    INITMSGA(MYDEBUG,0,
             "- aDim = "<<aDim<<
             "; aNbGauss = "<<aNbGauss<<
             "; aNbElem = "<<aNbElem<<
             "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
             std::endl);

    for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
      TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;

      TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
      TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
      TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
      TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);

      for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
        TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];

        for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
          TInt aNodeId = aConnSlice[aConnId] - 1;      
          TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);

          for(TInt aDimId = 0; aDimId < aDim; aDimId++){
            aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
          }
        }

        for(TInt aDimId = 0; aDimId < aDim; aDimId++){
          aGaussCoordSlice[aDimId] /= aNbNodes;
        }
      }
    }

    return true;
  }


  //---------------------------------------------------------------
  bool
  GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
                const TNodeInfo& theNodeInfo,
                TGaussCoord& theGaussCoord,
                const TElemNum& theElemNum,
                EModeSwitch theMode)
  {
    INITMSG(MYDEBUG,"GetBaryCenter\n");
    const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
    TInt aDim = aMeshInfo->GetDim();
    static TInt aNbGauss = 1;

    bool anIsSubMesh = !theElemNum.empty();
    TInt aNbElem;
    if(anIsSubMesh)
      aNbElem = theElemNum.size();
    else
      aNbElem = thePolyedreInfo.GetNbElem();

    theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);

    INITMSGA(MYDEBUG,0,
             "- aDim = "<<aDim<<
             "; aNbGauss = "<<aNbGauss<<
             "; aNbElem = "<<aNbElem<<
             "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
             std::endl);

    for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
      TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;

      TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
      TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
      TInt aNbFaces = aConnSliceArr.size();

      TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);

      for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
        TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];

        for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
          TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
          TInt aNbConn = aConnSlice.size();
          for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
            TInt aNodeId = aConnSlice[aConnId] - 1;      
            TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);

            for(TInt aDimId = 0; aDimId < aDim; aDimId++){
              aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
            }
          }
        }
        for(TInt aDimId = 0; aDimId < aDim; aDimId++){
          aGaussCoordSlice[aDimId] /= aNbNodes;
        }
      }
    }

    return true;
  }
}