Back to index

salome-med  6.5.0
MED_GaussUtils.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 //  File   : 
00024 //  Author : 
00025 //  Module : 
00026 //  $Header: /home/server/cvs/MED/MED_SRC/src/MEDWrapper/Base/MED_GaussUtils.cxx,v 1.6.2.2.6.3.10.2 2012-05-04 12:25:28 ouv Exp $
00027 //
00028 #include "MED_GaussUtils.hxx"
00029 #include "MED_Utilities.hxx"
00030  
00031 #ifdef _DEBUG_
00032 static int MYDEBUG = 0;
00033 static int MYVALUEDEBUG = 0;
00034 #else
00035 // static int MYDEBUG = 0;
00036 // static int MYVALUEDEBUG = 0;
00037 #endif
00038 
00039 //#define _DEBUG_REF_COORDS_
00040 
00041 namespace MED
00042 {
00043   //---------------------------------------------------------------
00044   TGaussCoord
00045   ::TGaussCoord():
00046     TModeSwitchInfo(eFULL_INTERLACE),
00047     myNbElem(0),
00048     myNbGauss(0),
00049     myDim(0),
00050     myGaussStep(0)
00051   {
00052   }
00053 
00054   void
00055   TGaussCoord
00056   ::Init(TInt theNbElem,
00057          TInt theNbGauss,
00058          TInt theDim,
00059          EModeSwitch theMode)
00060   {
00061     myModeSwitch = theMode;
00062 
00063     myNbElem = theNbElem;
00064     myNbGauss = theNbGauss;
00065     myDim = theDim;
00066 
00067     myGaussStep = myNbGauss*myDim;
00068 
00069     myGaussCoord.resize(theNbElem*myGaussStep);
00070   }
00071 
00072 
00073   TInt
00074   TGaussCoord
00075   ::GetNbElem() const
00076   { 
00077     return myNbElem; 
00078   }
00079   
00080   TInt
00081   TGaussCoord
00082   ::GetNbGauss() const
00083   { 
00084     return myNbGauss; 
00085   }
00086   
00087   TInt
00088   TGaussCoord
00089   ::GetDim() const
00090   { 
00091     return myDim; 
00092   }
00093   
00094   unsigned char*
00095   TGaussCoord
00096   ::GetValuePtr()
00097   {
00098     return (unsigned char*)&(myGaussCoord[0]);
00099   }
00100 
00101 
00102   TCCoordSliceArr 
00103   TGaussCoord
00104   ::GetCoordSliceArr(TInt theElemId) const
00105   {
00106     TCCoordSliceArr aCoordSliceArr(myNbGauss);
00107     if(GetModeSwitch() == eFULL_INTERLACE){
00108       TInt anId = theElemId*myGaussStep;
00109       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
00110         aCoordSliceArr[anGaussId] =
00111           TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
00112         anId += myDim;
00113       }
00114     }
00115     else{
00116       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
00117         aCoordSliceArr[anGaussId] =
00118           TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
00119       }
00120     }
00121     return aCoordSliceArr;
00122   }
00123 
00124 
00125   TCoordSliceArr 
00126   TGaussCoord
00127   ::GetCoordSliceArr(TInt theElemId)
00128   {
00129     TCoordSliceArr aCoordSliceArr(myNbGauss);
00130     if(GetModeSwitch() == eFULL_INTERLACE){
00131       TInt anId = theElemId*myGaussStep;
00132       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
00133         aCoordSliceArr[anGaussId] =
00134           TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
00135         anId += myDim;
00136       }
00137     }
00138     else{
00139       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
00140         aCoordSliceArr[anGaussId] =
00141           TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
00142       }
00143     }
00144     return aCoordSliceArr;
00145   }
00146 
00147 
00148   //---------------------------------------------------------------
00149   inline
00150   bool 
00151   IsEqual(TFloat theLeft, TFloat theRight)
00152   {
00153     static TFloat EPS = 1.0E-3;
00154     if(fabs(theLeft) + fabs(theRight) > EPS)
00155       return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
00156     return true;
00157   }
00158 
00159 
00160   //---------------------------------------------------------------
00161   class TShapeFun::TFun
00162   {
00163     TFloatVector myFun;
00164     TInt myNbRef;
00165 
00166   public:
00167 
00168     void
00169     Init(TInt theNbGauss,
00170          TInt theNbRef)
00171     {
00172       myFun.resize(theNbGauss*theNbRef);
00173       myNbRef = theNbRef;
00174     }
00175 
00176     TCFloatVecSlice 
00177     GetFunSlice(TInt theGaussId) const
00178     {
00179       return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
00180     }
00181 
00182     TFloatVecSlice
00183     GetFunSlice(TInt theGaussId)
00184     {
00185       return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
00186     }
00187   };
00188 
00189   //---------------------------------------------------------------
00190 
00191   TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
00192     myRefCoord(theNbRef*theDim),
00193     myDim(theDim),
00194     myNbRef(theNbRef)
00195   {}
00196 
00197   TCCoordSlice 
00198   TShapeFun::GetCoord(TInt theRefId) const
00199   {
00200     return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
00201   }
00202 
00203   TCoordSlice
00204   TShapeFun::GetCoord(TInt theRefId)
00205   {
00206     return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
00207   }
00208 
00209   void 
00210   TShapeFun::GetFun(const TCCoordSliceArr& theRef,
00211                     const TCCoordSliceArr& theGauss,
00212                     TFun& theFun) const
00213   {
00214     TInt aNbRef = theRef.size();
00215     TInt aNbGauss = theGauss.size();
00216     theFun.Init(aNbGauss,aNbRef);
00217   }
00218 
00219   bool 
00220   TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
00221   {
00222     TInt aNbRef = theRefCoord.size();
00223     TInt aNbRef2 = GetNbRef();
00224     INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
00225             "- aNbRef("<<aNbRef<<")"<<
00226             "; aNbRef2("<<aNbRef2<<")\n");
00227     bool anIsSatisfy = (aNbRef == aNbRef2);
00228     if(anIsSatisfy){
00229       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00230         const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
00231         TCCoordSlice aCoord = GetCoord(aRefId);
00232         TInt aDim = aCoord.size();
00233         bool anIsEqual = false;
00234         for(TInt anId = 0; anId < aDim; anId++){
00235           anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
00236           if(!anIsEqual){
00237             anIsSatisfy = false;
00238             break;
00239           }
00240         }
00241         if(!anIsEqual){
00242 #ifdef _DEBUG_
00243           TCCoordSlice aCoord = GetCoord(aRefId);
00244           INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
00245           TInt aDim = aCoord.size();
00246           for(TInt anId = 0; anId < aDim; anId++)
00247             ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
00248           const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
00249           ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
00250           for(TInt anId = 0; anId < aDim; anId++)
00251             ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
00252           ADDMSG(MYDEBUG,"}\n");
00253 #endif
00254 #ifndef _DEBUG_
00255           BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
00256           return anIsSatisfy;
00257 #endif
00258         }
00259       }
00260     }
00261 
00262     BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
00263     return anIsSatisfy;
00264   }
00265 
00266   bool
00267   TShapeFun::Eval(const TCellInfo&       theCellInfo,
00268                   const TNodeInfo&       theNodeInfo,
00269                   const TElemNum&        theElemNum,
00270                   const TCCoordSliceArr& theRef,
00271                   const TCCoordSliceArr& theGauss,
00272                   TGaussCoord&           theGaussCoord,
00273                   EModeSwitch            theMode)
00274   {
00275     INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);
00276 
00277     if(IsSatisfy(theRef)){
00278       const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
00279       TInt aDim = aMeshInfo->GetDim();
00280       TInt aNbGauss = theGauss.size();
00281 
00282       bool anIsSubMesh = !theElemNum.empty();
00283       TInt aNbElem;
00284       if(anIsSubMesh)
00285         aNbElem = theElemNum.size();
00286       else
00287         aNbElem = theCellInfo.GetNbElem();
00288 
00289       theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
00290 
00291       TFun aFun;
00292       InitFun(theRef,theGauss,aFun);
00293       TInt aConnDim = theCellInfo.GetConnDim();
00294 
00295       INITMSG(MYDEBUG,"aDim = "<<aDim<<
00296               "; aNbGauss = "<<aNbGauss<<
00297               "; aNbElem = "<<aNbElem<<
00298               "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
00299               std::endl);
00300 
00301       for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
00302         TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
00303         TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
00304         TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
00305 
00306         for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00307           TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
00308           TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
00309 
00310           for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
00311             TInt aNodeId = aConnSlice[aConnId] - 1;      
00312             TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
00313 
00314             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
00315               aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
00316             }
00317           }
00318         }
00319       }
00320 
00321 #ifdef _DEBUG_
00322       {
00323         INITMSG(MYVALUEDEBUG,"theGauss: ");
00324         for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00325           TCCoordSlice aCoordSlice = theGauss[aGaussId];
00326           ADDMSG(MYVALUEDEBUG,"{");
00327           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
00328             ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
00329           }
00330           ADDMSG(MYVALUEDEBUG,"} ");
00331         }
00332         ADDMSG(MYVALUEDEBUG,std::endl);
00333       }
00334       for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
00335         TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
00336         INITMSG(MYVALUEDEBUG,"");
00337         for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00338           TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
00339           ADDMSG(MYVALUEDEBUG,"{");
00340           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
00341             ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
00342           }
00343           ADDMSG(MYVALUEDEBUG,"} ");
00344         }
00345         ADDMSG(MYVALUEDEBUG,std::endl);
00346       }
00347 #endif
00348       return true;
00349     }
00350 
00351     return false;
00352   }
00353 
00354 
00355   //---------------------------------------------------------------
00356   TSeg2a::TSeg2a():TShapeFun(1,2)
00357   {
00358     TInt aNbRef = GetNbRef();
00359     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00360       TCoordSlice aCoord = GetCoord(aRefId);
00361       switch(aRefId){
00362       case  0: aCoord[0] = -1.0; break;
00363       case  1: aCoord[0] =  1.0; break;
00364       }
00365     }
00366   }
00367 
00368   void
00369   TSeg2a::InitFun(const TCCoordSliceArr& theRef,
00370                   const TCCoordSliceArr& theGauss,
00371                   TFun& theFun) const
00372   {
00373     GetFun(theRef,theGauss,theFun);
00374 
00375     TInt aNbGauss = theGauss.size();
00376     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00377       const TCCoordSlice& aCoord = theGauss[aGaussId];
00378       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00379 
00380       aSlice[0] = 0.5*(1.0 - aCoord[0]);
00381       aSlice[1] = 0.5*(1.0 + aCoord[0]);
00382     }
00383   }
00384 
00385 
00386   //---------------------------------------------------------------
00387   TSeg3a::TSeg3a():TShapeFun(1,3)
00388   {
00389     TInt aNbRef = GetNbRef();
00390     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00391       TCoordSlice aCoord = GetCoord(aRefId);
00392       switch(aRefId){
00393       case  0: aCoord[0] = -1.0; break;
00394       case  1: aCoord[0] =  1.0; break;
00395       case  2: aCoord[0] =  0.0; break;
00396       }
00397     }
00398   }
00399 
00400   void
00401   TSeg3a::InitFun(const TCCoordSliceArr& theRef,
00402                   const TCCoordSliceArr& theGauss,
00403                   TFun& theFun) const
00404   {
00405     GetFun(theRef,theGauss,theFun);
00406 
00407     TInt aNbGauss = theGauss.size();
00408     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00409       const TCCoordSlice& aCoord = theGauss[aGaussId];
00410       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00411 
00412       aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
00413       aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
00414       aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
00415     }
00416   }
00417 
00418 
00419 
00420   //---------------------------------------------------------------
00421   TTria3a::TTria3a():
00422     TShapeFun(2,3)
00423   {
00424     TInt aNbRef = GetNbRef();
00425     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00426       TCoordSlice aCoord = GetCoord(aRefId);
00427       switch(aRefId){
00428       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00429       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00430       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00431       }
00432     }
00433   }
00434 
00435   void
00436   TTria3a::InitFun(const TCCoordSliceArr& theRef,
00437                    const TCCoordSliceArr& theGauss,
00438                    TFun& theFun) const
00439   {
00440     GetFun(theRef,theGauss,theFun);
00441 
00442     TInt aNbGauss = theGauss.size();
00443     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00444       const TCCoordSlice& aCoord = theGauss[aGaussId];
00445       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00446 
00447       aSlice[0] = 0.5*(1.0 + aCoord[1]);
00448       aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
00449       aSlice[2] = 0.5*(1.0 + aCoord[0]);
00450     }
00451   }
00452 
00453 
00454 
00455   //---------------------------------------------------------------
00456   TTria6a::TTria6a():TShapeFun(2,6)
00457   {
00458     TInt aNbRef = GetNbRef();
00459     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00460       TCoordSlice aCoord = GetCoord(aRefId);
00461       switch(aRefId){
00462       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00463       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00464       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00465 
00466       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00467       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00468       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00469       }
00470     }
00471   }
00472 
00473   void
00474   TTria6a::InitFun(const TCCoordSliceArr& theRef,
00475                    const TCCoordSliceArr& theGauss,
00476                    TFun& theFun) const
00477   {
00478     GetFun(theRef,theGauss,theFun);
00479 
00480     TInt aNbGauss = theGauss.size();
00481     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00482       const TCCoordSlice& aCoord = theGauss[aGaussId];
00483       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00484 
00485       aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
00486       aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
00487       aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
00488 
00489       aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
00490       aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
00491       aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
00492     }
00493   }
00494 
00495 
00496 
00497   //---------------------------------------------------------------
00498   TTria3b::TTria3b():
00499     TShapeFun(2,3)
00500   {
00501     TInt aNbRef = GetNbRef();
00502     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00503       TCoordSlice aCoord = GetCoord(aRefId);
00504       switch(aRefId){
00505       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00506       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00507       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00508       }
00509     }
00510   }
00511 
00512   void
00513   TTria3b::InitFun(const TCCoordSliceArr& theRef,
00514                    const TCCoordSliceArr& theGauss,
00515                    TFun& theFun) const
00516   {
00517     GetFun(theRef,theGauss,theFun);
00518 
00519     TInt aNbGauss = theGauss.size();
00520     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00521       const TCCoordSlice& aCoord = theGauss[aGaussId];
00522       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00523 
00524       aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
00525       aSlice[1] = aCoord[0];
00526       aSlice[2] = aCoord[1];
00527     }
00528   }
00529 
00530 
00531 
00532   //---------------------------------------------------------------
00533   TTria6b::TTria6b():
00534     TShapeFun(2,6)
00535   {
00536     TInt aNbRef = GetNbRef();
00537     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00538       TCoordSlice aCoord = GetCoord(aRefId);
00539       switch(aRefId){
00540       case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00541       case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00542       case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00543 
00544       case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
00545       case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
00546       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
00547       }
00548     }
00549   }
00550 
00551   void
00552   TTria6b::InitFun(const TCCoordSliceArr& theRef,
00553                    const TCCoordSliceArr& theGauss,
00554                    TFun& theFun) const
00555   {
00556     GetFun(theRef,theGauss,theFun);
00557 
00558     TInt aNbGauss = theGauss.size();
00559     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00560       const TCCoordSlice& aCoord = theGauss[aGaussId];
00561       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00562 
00563       aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
00564       aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
00565       aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
00566 
00567       aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
00568       aSlice[4] = 4.0*aCoord[0]*aCoord[1];
00569       aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
00570     }
00571   }
00572 
00573 
00574 
00575   //---------------------------------------------------------------
00576   TQuad4a::TQuad4a():
00577     TShapeFun(2,4)
00578   {
00579     TInt aNbRef = GetNbRef();
00580     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00581       TCoordSlice aCoord = GetCoord(aRefId);
00582       switch(aRefId){
00583       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00584       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00585       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00586       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00587       }
00588     }
00589   }
00590 
00591   void
00592   TQuad4a::InitFun(const TCCoordSliceArr& theRef,
00593                    const TCCoordSliceArr& theGauss,
00594                    TFun& theFun) const
00595   {
00596     GetFun(theRef,theGauss,theFun);
00597 
00598     TInt aNbGauss = theGauss.size();
00599     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00600       const TCCoordSlice& aCoord = theGauss[aGaussId];
00601       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00602 
00603       aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
00604       aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
00605       aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
00606       aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
00607     }
00608   }
00609 
00610 
00611 
00612   //---------------------------------------------------------------
00613   TQuad8a::TQuad8a():
00614     TShapeFun(2,8)
00615   {
00616     TInt aNbRef = GetNbRef();
00617     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00618       TCoordSlice aCoord = GetCoord(aRefId);
00619       switch(aRefId){
00620       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00621       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00622       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00623       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00624 
00625       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
00626       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00627       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00628       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00629       }
00630     }
00631   }
00632 
00633   void
00634   TQuad8a::InitFun(const TCCoordSliceArr& theRef,
00635                    const TCCoordSliceArr& theGauss,
00636                    TFun& theFun) const
00637   {
00638     GetFun(theRef,theGauss,theFun);
00639 
00640     TInt aNbGauss = theGauss.size();
00641     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00642       const TCCoordSlice& aCoord = theGauss[aGaussId];
00643       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00644 
00645       aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
00646       aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
00647       aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
00648       aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);
00649 
00650       aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
00651       aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
00652       aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
00653       aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
00654     }
00655   }
00656 
00657 
00658 
00659   //---------------------------------------------------------------
00660   TQuad9a::TQuad9a():
00661     TShapeFun(2,9)
00662   {
00663     TInt aNbRef = GetNbRef();
00664     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00665       TCoordSlice aCoord = GetCoord(aRefId);
00666       switch(aRefId){
00667       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00668       case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00669       case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00670       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00671 
00672       case  4: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
00673       case  5: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00674       case  6: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00675       case  7: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00676 
00677       case  8: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00678       }
00679     }
00680   }
00681 
00682   void
00683   TQuad9a::InitFun(const TCCoordSliceArr& theRef,
00684                    const TCCoordSliceArr& theGauss,
00685                    TFun& theFun) const
00686   {
00687     GetFun(theRef,theGauss,theFun);
00688 
00689     TInt aNbGauss = theGauss.size();
00690     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00691       const TCCoordSlice& aCoord = theGauss[aGaussId];
00692       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00693 
00694       aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
00695       aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
00696       aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
00697       aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
00698 
00699       aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
00700       aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
00701       aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
00702       aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
00703 
00704       aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
00705     }
00706   }
00707 
00708 
00709 
00710   //---------------------------------------------------------------
00711   TQuad4b::TQuad4b():
00712     TShapeFun(2,4)
00713   {
00714     TInt aNbRef = GetNbRef();
00715     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00716       TCoordSlice aCoord = GetCoord(aRefId);
00717       switch(aRefId){
00718       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00719       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00720       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00721       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00722       }
00723     }
00724   }
00725 
00726   void
00727   TQuad4b::InitFun(const TCCoordSliceArr& theRef,
00728                    const TCCoordSliceArr& theGauss,
00729                    TFun& theFun) const
00730   {
00731     GetFun(theRef,theGauss,theFun);
00732 
00733     TInt aNbGauss = theGauss.size();
00734     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00735       const TCCoordSlice& aCoord = theGauss[aGaussId];
00736       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00737 
00738       aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
00739       aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
00740       aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
00741       aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
00742     }
00743   }
00744 
00745 
00746 
00747   //---------------------------------------------------------------
00748   TQuad8b::TQuad8b():
00749     TShapeFun(2,8)
00750   {
00751     TInt aNbRef = GetNbRef();
00752     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00753       TCoordSlice aCoord = GetCoord(aRefId);
00754       switch(aRefId){
00755       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00756       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00757       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00758       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00759 
00760       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00761       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00762       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00763       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
00764       }
00765     }
00766   }
00767 
00768   void
00769   TQuad8b::InitFun(const TCCoordSliceArr& theRef,
00770                    const TCCoordSliceArr& theGauss,
00771                    TFun& theFun) const
00772   {
00773     GetFun(theRef,theGauss,theFun);
00774 
00775     TInt aNbGauss = theGauss.size();
00776     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00777       const TCCoordSlice& aCoord = theGauss[aGaussId];
00778       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00779 
00780       aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
00781       aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
00782       aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
00783       aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);
00784 
00785       aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
00786       aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
00787       aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
00788       aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);
00789 
00790       //aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
00791       //aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
00792       //aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
00793       //aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
00794     }
00795   }
00796 
00797 
00798 
00799   //---------------------------------------------------------------
00800   TQuad9b::TQuad9b():
00801     TShapeFun(2,9)
00802   {
00803     TInt aNbRef = GetNbRef();
00804     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00805       TCoordSlice aCoord = GetCoord(aRefId);
00806       switch(aRefId){
00807       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
00808       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
00809       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
00810       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
00811 
00812       case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
00813       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
00814       case  6: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
00815       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0; break;
00816 
00817       case  8: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
00818       }
00819     }
00820   }
00821 
00822   void
00823   TQuad9b::InitFun(const TCCoordSliceArr& theRef,
00824                    const TCCoordSliceArr& theGauss,
00825                    TFun& theFun) const
00826   {
00827     GetFun(theRef,theGauss,theFun);
00828 
00829     TInt aNbGauss = theGauss.size();
00830     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00831       const TCCoordSlice& aCoord = theGauss[aGaussId];
00832       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00833 
00834       aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
00835       aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
00836       aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
00837       aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
00838 
00839       aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
00840       aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
00841       aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
00842       aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
00843 
00844       aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
00845     }
00846   }
00847 
00848 
00849 
00850   //---------------------------------------------------------------
00851   TTetra4a::TTetra4a():
00852     TShapeFun(3,4)
00853   {
00854     TInt aNbRef = GetNbRef();
00855     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00856       TCoordSlice aCoord = GetCoord(aRefId);
00857       switch(aRefId){
00858       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00859       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00860       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00861       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00862       }
00863     }
00864   }
00865 
00866   void
00867   TTetra4a::InitFun(const TCCoordSliceArr& theRef,
00868                     const TCCoordSliceArr& theGauss,
00869                     TFun& theFun) const
00870   {
00871     GetFun(theRef,theGauss,theFun);
00872 
00873     TInt aNbGauss = theGauss.size();
00874     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00875       const TCCoordSlice& aCoord = theGauss[aGaussId];
00876       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00877 
00878       aSlice[0] = aCoord[1];
00879       aSlice[1] = aCoord[2];
00880       aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
00881       aSlice[3] = aCoord[0];
00882     }
00883   }
00884 
00885 
00886 
00887   //---------------------------------------------------------------
00888   TTetra10a::TTetra10a():
00889     TShapeFun(3,10)
00890   {
00891     TInt aNbRef = GetNbRef();
00892     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00893       TCoordSlice aCoord = GetCoord(aRefId);
00894       switch(aRefId){
00895       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00896       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00897       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00898       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00899 
00900       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00901       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00902       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00903 
00904       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00905       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00906       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00907       }
00908     }
00909   }
00910 
00911   void
00912   TTetra10a::InitFun(const TCCoordSliceArr& theRef,
00913                      const TCCoordSliceArr& theGauss,
00914                      TFun& theFun) const
00915   {
00916     GetFun(theRef,theGauss,theFun);
00917 
00918     TInt aNbGauss = theGauss.size();
00919     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00920       const TCCoordSlice& aCoord = theGauss[aGaussId];
00921       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00922 
00923       aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
00924       aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
00925       aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
00926       aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
00927 
00928       aSlice[4] = 4.0*aCoord[1]*aCoord[2];
00929       aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
00930       aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
00931 
00932       aSlice[7] = 4.0*aCoord[0]*aCoord[1];
00933       aSlice[8] = 4.0*aCoord[0]*aCoord[2];
00934       aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
00935     }
00936   }
00937 
00938 
00939 
00940   //---------------------------------------------------------------
00941 
00942 
00943   TTetra4b::TTetra4b():
00944     TShapeFun(3,4)
00945   {
00946     TInt aNbRef = GetNbRef();
00947     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00948       TCoordSlice aCoord = GetCoord(aRefId);
00949       switch(aRefId){
00950       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00951       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00952       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00953       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00954       }
00955     }
00956   }
00957 
00958   void
00959   TTetra4b::InitFun(const TCCoordSliceArr& theRef,
00960                     const TCCoordSliceArr& theGauss,
00961                     TFun& theFun) const
00962   {
00963     GetFun(theRef,theGauss,theFun);
00964 
00965     TInt aNbGauss = theGauss.size();
00966     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
00967       const TCCoordSlice& aCoord = theGauss[aGaussId];
00968       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
00969 
00970       aSlice[0] = aCoord[1];
00971       aSlice[2] = aCoord[2];
00972       aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
00973       aSlice[3] = aCoord[0];
00974     }
00975   }
00976 
00977 
00978 
00979   //---------------------------------------------------------------
00980   TTetra10b::TTetra10b():
00981     TShapeFun(3,10)
00982   {
00983     TInt aNbRef = GetNbRef();
00984     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
00985       TCoordSlice aCoord = GetCoord(aRefId);
00986       switch(aRefId){
00987       case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
00988       case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
00989       case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00990       case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00991 
00992       case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
00993       case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00994       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00995 
00996       case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
00997       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
00998       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
00999       }
01000     }
01001   }
01002 
01003   void
01004   TTetra10b::InitFun(const TCCoordSliceArr& theRef,
01005                      const TCCoordSliceArr& theGauss,
01006                      TFun& theFun) const
01007   {
01008     GetFun(theRef,theGauss,theFun);
01009 
01010     TInt aNbGauss = theGauss.size();
01011     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01012       const TCCoordSlice& aCoord = theGauss[aGaussId];
01013       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01014 
01015       aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
01016       aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
01017       aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
01018       aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
01019 
01020       aSlice[6] = 4.0*aCoord[1]*aCoord[2];
01021       aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
01022       aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
01023 
01024       aSlice[7] = 4.0*aCoord[0]*aCoord[1];
01025       aSlice[9] = 4.0*aCoord[0]*aCoord[2];
01026       aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
01027     }
01028   }
01029 
01030 
01031 
01032   //---------------------------------------------------------------
01033   THexa8a::THexa8a():
01034     TShapeFun(3,8)
01035   {
01036     TInt aNbRef = GetNbRef();
01037     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01038       TCoordSlice aCoord = GetCoord(aRefId);
01039       switch(aRefId){
01040       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01041       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01042       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01043       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01044       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01045       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01046       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01047       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01048       }
01049     }
01050   }
01051 
01052   void
01053   THexa8a::InitFun(const TCCoordSliceArr& theRef,
01054                    const TCCoordSliceArr& theGauss,
01055                    TFun& theFun) const
01056   {
01057     GetFun(theRef,theGauss,theFun);
01058 
01059     TInt aNbGauss = theGauss.size();
01060     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01061       const TCCoordSlice& aCoord = theGauss[aGaussId];
01062       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01063 
01064       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
01065       aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
01066       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
01067       aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
01068 
01069       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
01070       aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
01071       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
01072       aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
01073     }
01074   }
01075 
01076 
01077   //---------------------------------------------------------------
01078   THexa20a::THexa20a(TInt theDim, TInt theNbRef):
01079     TShapeFun(theDim,theNbRef)
01080   {
01081     TInt aNbRef = myRefCoord.size();
01082     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01083       TCoordSlice aCoord = GetCoord(aRefId);
01084       switch(aRefId){
01085       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01086       case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01087       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01088       case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01089       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01090       case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01091       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01092       case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01093 
01094       case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01095       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
01096       case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01097       case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
01098       case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01099       case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01100       case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01101       case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01102       case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01103       case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01104       case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01105       case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01106       }
01107     }
01108   }
01109 
01110   void
01111   THexa20a::InitFun(const TCCoordSliceArr& theRef,
01112                     const TCCoordSliceArr& theGauss,
01113                     TFun& theFun) const
01114   {
01115     GetFun(theRef,theGauss,theFun);
01116 
01117     TInt aNbGauss = theGauss.size();
01118     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01119       const TCCoordSlice& aCoord = theGauss[aGaussId];
01120       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01121 
01122       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
01123         (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
01124       aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
01125         (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
01126       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
01127         (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
01128       aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
01129         (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
01130       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
01131         (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
01132       aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
01133         (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
01134       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
01135         (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
01136       aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
01137         (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
01138 
01139       aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
01140       aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
01141       aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
01142       aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
01143       aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
01144       aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
01145       aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
01146       aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
01147       aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
01148       aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
01149       aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
01150       aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
01151     }
01152   }
01153 
01154 
01155 
01156   //---------------------------------------------------------------
01157   THexa27a::THexa27a():
01158     THexa20a(3,27)
01159   {
01160     TInt aNbRef = myRefCoord.size();
01161     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01162       TCoordSlice aCoord = GetCoord(aRefId);
01163       switch(aRefId){
01164       case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
01165       case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01166       case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01167       case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01168       case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01169       case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01170       case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01171       }
01172     }
01173   }
01174 
01175   void
01176   THexa27a::InitFun(const TCCoordSliceArr& theRef,
01177                     const TCCoordSliceArr& theGauss,
01178                     TFun& theFun) const
01179   {
01180     GetFun(theRef,theGauss,theFun);
01181 
01182     TInt aNbGauss = theGauss.size();
01183     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01184       const TCCoordSlice& aCoord = theGauss[aGaussId];
01185       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01186 
01187       aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
01188       aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
01189       aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
01190       aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
01191       aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
01192       aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
01193       aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
01194       aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
01195 
01196       aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
01197       aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
01198       aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
01199       aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
01200       aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
01201       aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
01202       aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
01203       aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
01204       aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
01205       aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
01206       aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
01207       aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
01208       aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
01209       aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
01210       aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
01211       aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
01212       aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
01213       aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
01214       aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
01215     }
01216   }
01217 
01218 
01219 
01220   //---------------------------------------------------------------
01221   THexa8b::THexa8b():
01222     TShapeFun(3,8)
01223   {
01224     TInt aNbRef = GetNbRef();
01225     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01226       TCoordSlice aCoord = GetCoord(aRefId);
01227       switch(aRefId){
01228       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01229       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01230       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01231       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01232       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01233       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01234       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01235       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01236       }
01237     }
01238   }
01239 
01240   void
01241   THexa8b::InitFun(const TCCoordSliceArr& theRef,
01242                    const TCCoordSliceArr& theGauss,
01243                    TFun& theFun) const
01244   {
01245     GetFun(theRef,theGauss,theFun);
01246 
01247     TInt aNbGauss = theGauss.size();
01248     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01249       const TCCoordSlice& aCoord = theGauss[aGaussId];
01250       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01251 
01252       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
01253       aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
01254       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
01255       aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
01256 
01257       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
01258       aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
01259       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
01260       aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
01261     }
01262   }
01263 
01264 
01265 
01266   //---------------------------------------------------------------
01267   THexa20b::THexa20b(TInt theDim, TInt theNbRef):
01268     TShapeFun(theDim,theNbRef)
01269   {
01270     TInt aNbRef = myRefCoord.size();
01271     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01272       TCoordSlice aCoord = GetCoord(aRefId);
01273       switch(aRefId){
01274       case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01275       case  3: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01276       case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01277       case  1: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01278       case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01279       case  7: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01280       case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01281       case  5: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01282 
01283       case 11: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
01284       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
01285       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
01286       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
01287       case 16: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01288       case 19: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01289       case 18: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01290       case 17: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01291       case 15: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
01292       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01293       case 13: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
01294       case 12: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01295       }
01296     }
01297   }
01298 
01299   void
01300   THexa20b::InitFun(const TCCoordSliceArr& theRef,
01301                     const TCCoordSliceArr& theGauss,
01302                     TFun& theFun) const
01303   {
01304     GetFun(theRef,theGauss,theFun);
01305 
01306     TInt aNbGauss = theGauss.size();
01307     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01308       const TCCoordSlice& aCoord = theGauss[aGaussId];
01309       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01310 
01311       aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
01312         (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
01313       aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
01314         (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
01315       aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
01316         (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
01317       aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
01318         (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
01319       aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
01320         (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
01321       aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
01322         (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
01323       aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
01324         (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
01325       aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
01326         (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
01327 
01328       aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
01329       aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
01330       aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
01331       aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
01332       aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
01333       aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
01334       aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
01335       aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
01336       aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
01337       aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
01338       aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
01339       aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
01340     }
01341   }
01342 
01343 
01344 
01345   //---------------------------------------------------------------
01346   TPenta6a::TPenta6a():
01347     TShapeFun(3,6)
01348   {
01349     TInt aNbRef = myRefCoord.size();
01350     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01351       TCoordSlice aCoord = GetCoord(aRefId);
01352       switch(aRefId){
01353       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01354       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
01355       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01356       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01357       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01358       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01359       }
01360     }
01361   }
01362 
01363   void
01364   TPenta6a::InitFun(const TCCoordSliceArr& theRef,
01365                     const TCCoordSliceArr& theGauss,
01366                     TFun& theFun) const
01367   {
01368     GetFun(theRef,theGauss,theFun);
01369 
01370     TInt aNbGauss = theGauss.size();
01371     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01372       const TCCoordSlice& aCoord = theGauss[aGaussId];
01373       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01374 
01375       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
01376       aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
01377       aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
01378 
01379       aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
01380       aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
01381       aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
01382     }
01383   }
01384 
01385 
01386 
01387   //---------------------------------------------------------------
01388   TPenta6b::TPenta6b():
01389     TShapeFun(3,6)
01390   {
01391     TInt aNbRef = myRefCoord.size();
01392     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01393       TCoordSlice aCoord = GetCoord(aRefId);
01394       switch(aRefId){
01395       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01396       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
01397       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01398       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01399       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01400       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01401       }
01402     }
01403   }
01404 
01405   void
01406   TPenta6b::InitFun(const TCCoordSliceArr& theRef,
01407                     const TCCoordSliceArr& theGauss,
01408                     TFun& theFun) const
01409   {
01410     GetFun(theRef,theGauss,theFun);
01411 
01412     TInt aNbGauss = theGauss.size();
01413     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01414       const TCCoordSlice& aCoord = theGauss[aGaussId];
01415       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01416 
01417       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
01418       aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
01419       aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
01420 
01421       aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
01422       aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
01423       aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
01424     }
01425   }
01426 
01427 
01428 
01429   //---------------------------------------------------------------
01430   TPenta15a::TPenta15a():
01431     TShapeFun(3,15)
01432   {
01433     TInt aNbRef = myRefCoord.size();
01434     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01435       TCoordSlice aCoord = GetCoord(aRefId);
01436       switch(aRefId){
01437       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01438       case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
01439       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01440       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01441       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01442       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01443 
01444       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
01445       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01446       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01447       case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01448       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01449       case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01450       case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
01451       case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01452       case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01453       }
01454     }
01455   }
01456 
01457   void
01458   TPenta15a::InitFun(const TCCoordSliceArr& theRef,
01459                      const TCCoordSliceArr& theGauss,
01460                      TFun& theFun) const
01461   {
01462     GetFun(theRef,theGauss,theFun);
01463 
01464     TInt aNbGauss = theGauss.size();
01465     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01466       const TCCoordSlice& aCoord = theGauss[aGaussId];
01467       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01468 
01469       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
01470       aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
01471       aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
01472 
01473       aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
01474       aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
01475       aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
01476 
01477       aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
01478       aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
01479       aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
01480 
01481       aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
01482       aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
01483       aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
01484 
01485       aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
01486       aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
01487       aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
01488     }
01489   }
01490 
01491 
01492 
01493   //---------------------------------------------------------------
01494   TPenta15b::TPenta15b():
01495     TShapeFun(3,15)
01496   {
01497     TInt aNbRef = myRefCoord.size();
01498     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01499       TCoordSlice aCoord = GetCoord(aRefId);
01500       switch(aRefId){
01501       case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01502       case  2: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
01503       case  1: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01504       case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01505       case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01506       case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01507 
01508       case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
01509       case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01510       case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01511       case 12: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01512       case 14: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01513       case 13: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01514       case 11: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
01515       case 10: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01516       case  9: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01517       }
01518     }
01519   }
01520 
01521   void
01522   TPenta15b::InitFun(const TCCoordSliceArr& theRef,
01523                      const TCCoordSliceArr& theGauss,
01524                      TFun& theFun) const
01525   {
01526     GetFun(theRef,theGauss,theFun);
01527 
01528     TInt aNbGauss = theGauss.size();
01529     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01530       const TCCoordSlice& aCoord = theGauss[aGaussId];
01531       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01532 
01533       aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
01534       aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
01535       aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
01536 
01537       aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
01538       aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
01539       aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
01540 
01541       aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
01542       aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
01543       aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
01544 
01545       aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
01546       aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
01547       aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
01548 
01549       aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
01550       aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
01551       aSlice[9]  = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
01552     }
01553   }
01554 
01555 
01556 
01557   //---------------------------------------------------------------
01558   TPyra5a::TPyra5a():
01559     TShapeFun(3,5)
01560   {
01561     TInt aNbRef = myRefCoord.size();
01562     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01563       TCoordSlice aCoord = GetCoord(aRefId);
01564       switch(aRefId){
01565       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01566       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01567       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01568       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01569       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01570       }
01571     }
01572   }
01573 
01574   void
01575   TPyra5a::InitFun(const TCCoordSliceArr& theRef,
01576                    const TCCoordSliceArr& theGauss,
01577                    TFun& theFun) const
01578   {
01579     GetFun(theRef,theGauss,theFun);
01580 
01581     TInt aNbGauss = theGauss.size();
01582     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01583       const TCCoordSlice& aCoord = theGauss[aGaussId];
01584       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01585       // APO & RNV:
01586       // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
01587       // Seems shape function for ePYRA5 elements is:
01588       // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
01589       // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
01590       // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
01591       // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
01592       // w4 = +Z;
01593       aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01594       aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01595       aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01596       aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01597       aSlice[4] = aCoord[2];
01598     }
01599   }
01600 
01601 
01602 
01603   //---------------------------------------------------------------
01604   TPyra5b::TPyra5b():
01605     TShapeFun(3,5)
01606   {
01607     TInt aNbRef = myRefCoord.size();
01608     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01609       TCoordSlice aCoord = GetCoord(aRefId);
01610       switch(aRefId){        
01611       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01612       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01613       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01614       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01615       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01616       }
01617     }
01618   }
01619 
01620   void
01621   TPyra5b::InitFun(const TCCoordSliceArr& theRef,
01622                    const TCCoordSliceArr& theGauss,
01623                    TFun& theFun) const
01624   {
01625     GetFun(theRef,theGauss,theFun);
01626     
01627     TInt aNbGauss = theGauss.size();
01628     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01629       const TCCoordSlice& aCoord = theGauss[aGaussId];
01630       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01631       // APO & RNV:
01632       // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
01633       // Seems shape function for ePYRA5 elements is:
01634       // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
01635       // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
01636       // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
01637       // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
01638       // w4 = +Z;
01639       aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01640       aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01641       aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01642       aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
01643       aSlice[4] = aCoord[2];
01644     }
01645   }
01646   
01647 
01648 
01649   //---------------------------------------------------------------
01650   TPyra13a::TPyra13a():
01651     TShapeFun(3,13)
01652   {
01653     TInt aNbRef = myRefCoord.size();
01654     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01655       TCoordSlice aCoord = GetCoord(aRefId);
01656       switch(aRefId){
01657       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01658       case  1: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01659       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01660       case  3: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01661       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01662 
01663       case  5: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01664       case  6: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01665       case  7: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
01666       case  8: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
01667       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01668       case 10: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
01669       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01670       case 12: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
01671       }
01672     }
01673   }
01674 
01675   void
01676   TPyra13a::InitFun(const TCCoordSliceArr& theRef,
01677                     const TCCoordSliceArr& theGauss,
01678                     TFun& theFun) const
01679   {
01680     GetFun(theRef,theGauss,theFun);
01681 
01682     TInt aNbGauss = theGauss.size();
01683     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01684       const TCCoordSlice& aCoord = theGauss[aGaussId];
01685       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01686 
01687       aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01688         (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
01689       aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01690         (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
01691       aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01692         (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
01693       aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01694         (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
01695 
01696       aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
01697 
01698       aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01699         (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01700       aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01701         (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01702       aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01703         (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01704       aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01705         (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01706 
01707       aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
01708         (1.0 - aCoord[2]);
01709       aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
01710         (1.0 - aCoord[2]);
01711       aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
01712         (1.0 - aCoord[2]);
01713       aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
01714         (1.0 - aCoord[2]);
01715     }
01716   }
01717 
01718 
01719 
01720   //---------------------------------------------------------------
01721   TPyra13b::TPyra13b():
01722     TShapeFun(3,13)
01723   {
01724     TInt aNbRef = myRefCoord.size();
01725     for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
01726       TCoordSlice aCoord = GetCoord(aRefId);
01727       switch(aRefId){
01728       case  0: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01729       case  3: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
01730       case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
01731       case  1: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
01732       case  4: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
01733 
01734       case  8: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01735       case  7: aCoord[0] = -0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
01736       case  6: aCoord[0] = -0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
01737       case  5: aCoord[0] =  0.5;  aCoord[1] = -0.5;  aCoord[2] =  0.0; break;
01738       case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01739       case 12: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
01740       case 11: aCoord[0] = -0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
01741       case 10: aCoord[0] =  0.0;  aCoord[1] = -0.5;  aCoord[2] =  0.5; break;
01742       }
01743     }
01744   }
01745 
01746   void
01747   TPyra13b::InitFun(const TCCoordSliceArr& theRef,
01748                     const TCCoordSliceArr& theGauss,
01749                     TFun& theFun) const
01750   {
01751     GetFun(theRef,theGauss,theFun);
01752 
01753     TInt aNbGauss = theGauss.size();
01754     for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
01755       const TCCoordSlice& aCoord = theGauss[aGaussId];
01756       TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
01757 
01758       aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01759         (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
01760       aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01761         (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
01762       aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01763         (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
01764       aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01765         (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
01766 
01767       aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
01768 
01769       aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01770         (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01771       aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
01772         (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01773       aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01774         (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01775       aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
01776         (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
01777 
01778       aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
01779         (1.0 - aCoord[2]);
01780       aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
01781         (1.0 - aCoord[2]);
01782       aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
01783         (1.0 - aCoord[2]);
01784       aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
01785         (1.0 - aCoord[2]);
01786     }
01787   }
01788 
01789 
01790 
01791   //---------------------------------------------------------------
01792   bool
01793   GetGaussCoord3D(const TGaussInfo& theGaussInfo, 
01794                   const TCellInfo& theCellInfo,
01795                   const TNodeInfo& theNodeInfo,
01796                   TGaussCoord& theGaussCoord,
01797                   const TElemNum& theElemNum,
01798                   EModeSwitch theMode)
01799   {
01800     INITMSG(MYDEBUG,"GetGaussCoord3D\n");
01801 
01802     if(theGaussInfo.myGeom == theCellInfo.myGeom){
01803       EGeometrieElement aGeom = theGaussInfo.myGeom;
01804 
01805       TInt aNbRef = theGaussInfo.GetNbRef();
01806       TCCoordSliceArr aRefSlice(aNbRef);
01807       for(TInt anId = 0; anId < aNbRef; anId++)
01808         aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
01809 
01810       TInt aNbGauss = theGaussInfo.GetNbGauss();
01811       TCCoordSliceArr aGaussSlice(aNbGauss);
01812       for(TInt anId = 0; anId < aNbGauss; anId++)
01813         aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
01814 
01815       switch(aGeom){
01816       case eSEG2: {
01817         INITMSG(MYDEBUG,"eSEG2"<<std::endl);
01818 
01819         if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01820           return true;
01821 
01822         break;
01823       }
01824       case eSEG3: {
01825         INITMSG(MYDEBUG,"eSEG3"<<std::endl);
01826 
01827         if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01828           return true;
01829 
01830         break;
01831       }
01832       case eTRIA3: {
01833         INITMSG(MYDEBUG,"eTRIA3"<<std::endl);
01834 
01835         if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01836           return true;
01837 
01838         if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01839           return true;
01840 
01841         break;
01842       }
01843       case eTRIA6: {
01844         INITMSG(MYDEBUG,"eTRIA6"<<std::endl);
01845 
01846         if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01847           return true;
01848 
01849         if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01850           return true;
01851 
01852         break;
01853       }
01854       case eQUAD4: {
01855         INITMSG(MYDEBUG,"eQUAD4"<<std::endl);
01856 
01857         if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01858           return true;
01859 
01860         if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01861           return true;
01862 
01863         break;
01864       }
01865       case eQUAD8: {
01866         INITMSG(MYDEBUG,"eQUAD8"<<std::endl);
01867 
01868         if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01869           return true;
01870 
01871         if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01872           return true;
01873 
01874         break;
01875       }
01876       case eQUAD9: {
01877         INITMSG(MYDEBUG,"eQUAD9"<<std::endl);
01878 
01879         if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01880           return true;
01881 
01882         if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01883           return true;
01884 
01885         break;
01886       }
01887       case eTETRA4: {
01888         INITMSG(MYDEBUG,"eTETRA4"<<std::endl);
01889 
01890         if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01891           return true;
01892 
01893         if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01894           return true;
01895 
01896         break;
01897       }
01898       case ePYRA5: {
01899         INITMSG(MYDEBUG,"ePYRA5"<<std::endl);
01900 
01901         if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01902           return true;
01903 
01904         if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01905           return true;
01906 
01907         break;
01908       }
01909       case ePENTA6: {
01910         INITMSG(MYDEBUG,"ePENTA6"<<std::endl);
01911 
01912         if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01913           return true;
01914 
01915         if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01916           return true;
01917 
01918         break;
01919       }
01920       case eHEXA8: {
01921         INITMSG(MYDEBUG,"eHEXA8"<<std::endl);
01922 
01923         if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01924           return true;
01925 
01926         if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01927           return true;
01928 
01929         break;
01930       }
01931       case eTETRA10: {
01932         INITMSG(MYDEBUG,"eTETRA10"<<std::endl);
01933 
01934         if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01935           return true;
01936 
01937         if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01938           return true;
01939 
01940         break;
01941       }
01942       case ePYRA13: {
01943         INITMSG(MYDEBUG,"ePYRA13"<<std::endl);
01944 
01945         if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01946           return true;
01947 
01948         if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01949           return true;
01950 
01951         break;
01952       }
01953       case ePENTA15: {
01954         INITMSG(MYDEBUG,"ePENTA15"<<std::endl);
01955 
01956         if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01957           return true;
01958 
01959         if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01960           return true;
01961 
01962         break;
01963       }
01964       case eHEXA20: {
01965         INITMSG(MYDEBUG,"eHEXA20"<<std::endl);
01966 
01967         if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01968           return true;
01969 
01970         if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
01971           return true;
01972 
01973         break;
01974       }
01975       default: 
01976         INITMSG(MYDEBUG,"eNONE"<<std::endl);
01977         return false;
01978       }
01979     }
01980 
01981     return false;
01982   }
01983 
01984   //---------------------------------------------------------------
01985   bool
01986   GetBaryCenter(const TCellInfo& theCellInfo,
01987                 const TNodeInfo& theNodeInfo,
01988                 TGaussCoord& theGaussCoord,
01989                 const TElemNum& theElemNum,
01990                 EModeSwitch theMode)
01991   {
01992     INITMSG(MYDEBUG,"GetBaryCenter\n");
01993     const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
01994     TInt aDim = aMeshInfo->GetDim();
01995     static TInt aNbGauss = 1;
01996 
01997     bool anIsSubMesh = !theElemNum.empty();
01998     TInt aNbElem;
01999     if(anIsSubMesh)
02000       aNbElem = theElemNum.size();
02001     else
02002       aNbElem = theCellInfo.GetNbElem();
02003 
02004     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
02005 
02006     TInt aConnDim = theCellInfo.GetConnDim();
02007 
02008     INITMSGA(MYDEBUG,0,
02009              "- aDim = "<<aDim<<
02010              "; aNbGauss = "<<aNbGauss<<
02011              "; aNbElem = "<<aNbElem<<
02012              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
02013              std::endl);
02014 
02015     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
02016       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
02017       TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
02018       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
02019 
02020       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
02021         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
02022 
02023         for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
02024           TInt aNodeId = aConnSlice[aConnId] - 1;      
02025           TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
02026 
02027           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02028             aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
02029           }
02030         }
02031 
02032         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02033           aGaussCoordSlice[aDimId] /= aConnDim;
02034         }
02035       }
02036     }
02037 
02038 #ifdef _DEBUG_
02039     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
02040       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
02041       INITMSG(MYVALUEDEBUG,"");
02042       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
02043         TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
02044         ADDMSG(MYVALUEDEBUG,"{");
02045         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02046           ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
02047         }
02048         ADDMSG(MYVALUEDEBUG,"} ");
02049       }
02050       ADDMSG(MYVALUEDEBUG,std::endl);
02051     }
02052 #endif
02053 
02054     return true;
02055   }
02056 
02057 
02058   //---------------------------------------------------------------
02059   bool
02060   GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
02061                 const TNodeInfo& theNodeInfo,
02062                 TGaussCoord& theGaussCoord,
02063                 const TElemNum& theElemNum,
02064                 EModeSwitch theMode)
02065   {
02066     INITMSG(MYDEBUG,"GetBaryCenter\n");
02067     const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
02068     TInt aDim = aMeshInfo->GetDim();
02069     static TInt aNbGauss = 1;
02070 
02071     bool anIsSubMesh = !theElemNum.empty();
02072     TInt aNbElem;
02073     if(anIsSubMesh)
02074       aNbElem = theElemNum.size();
02075     else
02076       aNbElem = thePolygoneInfo.GetNbElem();
02077 
02078     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
02079 
02080     INITMSGA(MYDEBUG,0,
02081              "- aDim = "<<aDim<<
02082              "; aNbGauss = "<<aNbGauss<<
02083              "; aNbElem = "<<aNbElem<<
02084              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
02085              std::endl);
02086 
02087     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
02088       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
02089 
02090       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
02091       TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
02092       TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
02093       TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
02094 
02095       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
02096         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
02097 
02098         for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
02099           TInt aNodeId = aConnSlice[aConnId] - 1;      
02100           TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
02101 
02102           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02103             aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
02104           }
02105         }
02106 
02107         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02108           aGaussCoordSlice[aDimId] /= aNbNodes;
02109         }
02110       }
02111     }
02112 
02113     return true;
02114   }
02115 
02116 
02117   //---------------------------------------------------------------
02118   bool
02119   GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
02120                 const TNodeInfo& theNodeInfo,
02121                 TGaussCoord& theGaussCoord,
02122                 const TElemNum& theElemNum,
02123                 EModeSwitch theMode)
02124   {
02125     INITMSG(MYDEBUG,"GetBaryCenter\n");
02126     const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
02127     TInt aDim = aMeshInfo->GetDim();
02128     static TInt aNbGauss = 1;
02129 
02130     bool anIsSubMesh = !theElemNum.empty();
02131     TInt aNbElem;
02132     if(anIsSubMesh)
02133       aNbElem = theElemNum.size();
02134     else
02135       aNbElem = thePolyedreInfo.GetNbElem();
02136 
02137     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
02138 
02139     INITMSGA(MYDEBUG,0,
02140              "- aDim = "<<aDim<<
02141              "; aNbGauss = "<<aNbGauss<<
02142              "; aNbElem = "<<aNbElem<<
02143              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
02144              std::endl);
02145 
02146     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
02147       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
02148 
02149       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
02150       TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
02151       TInt aNbFaces = aConnSliceArr.size();
02152 
02153       TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
02154 
02155       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
02156         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
02157 
02158         for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
02159           TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
02160           TInt aNbConn = aConnSlice.size();
02161           for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
02162             TInt aNodeId = aConnSlice[aConnId] - 1;      
02163             TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
02164 
02165             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02166               aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
02167             }
02168           }
02169         }
02170         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
02171           aGaussCoordSlice[aDimId] /= aNbNodes;
02172         }
02173       }
02174     }
02175 
02176     return true;
02177   }
02178 }