src/MEDWrapper/Base/MED_GaussUtils.cxx

   1 //
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   :
  25 //  Author :
  26 //  Module :
  27 //  $Header$
  28
  29 #include "MED_GaussUtils.hxx"
  30 #include "MED_Utilities.hxx"
  31
  32 #ifdef _DEBUG_
  33 static int MYDEBUG = 0;
  34 static int MYVALUEDEBUG = 0;
  35 #else
  36 static int MYDEBUG = 0;
  37 static int MYVALUEDEBUG = 0;
  38 #endif
  39
  40 //#define _DEBUG_REF_COORDS_
  41
  42 namespace MED
  43 {
  44   //---------------------------------------------------------------
  45   TGaussCoord
  46   ::TGaussCoord():
  47     TModeSwitchInfo(eFULL_INTERLACE),
  48     myNbElem(0),
  49     myNbGauss(0),
  50     myDim(0),
  51     myGaussStep(0)
  52   {
  53   }
  54
  55   void
  56   TGaussCoord
  57   ::Init(TInt theNbElem,
  58          TInt theNbGauss,
  59          TInt theDim,
  60          EModeSwitch theMode)
  61   {
  62     myModeSwitch = theMode;
  63
  64     myNbElem = theNbElem;
  65     myNbGauss = theNbGauss;
  66     myDim = theDim;
  67
  68     myGaussStep = myNbGauss*myDim;
  69
  70     myGaussCoord.resize(theNbElem*myGaussStep);
  71   }
  72
  73
  74   TCCoordSliceArr
  75   TGaussCoord
  76   ::GetCoordSliceArr(TInt theElemId) const
  77   {
  78     TCCoordSliceArr aCoordSliceArr(myNbGauss);
  79     if(GetModeSwitch() == eFULL_INTERLACE){
  80       TInt anId = theElemId*myGaussStep;
  81       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
  82         aCoordSliceArr[anGaussId] =
  83           TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
  84         anId += myDim;
  85       }
  86     }
  87     else{
  88       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
  89         aCoordSliceArr[anGaussId] =
  90           TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
  91       }
  92     }
  93     return aCoordSliceArr;
  94   }
  95
  96
  97   TCoordSliceArr
  98   TGaussCoord
  99   ::GetCoordSliceArr(TInt theElemId)
 100   {
 101     TCoordSliceArr aCoordSliceArr(myNbGauss);
 102     if(GetModeSwitch() == eFULL_INTERLACE){
 103       TInt anId = theElemId*myGaussStep;
 104       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
 105         aCoordSliceArr[anGaussId] =
 106           TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
 107         anId += myDim;
 108       }
 109     }
 110     else{
 111       for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
 112         aCoordSliceArr[anGaussId] =
 113           TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
 114       }
 115     }
 116     return aCoordSliceArr;
 117   }
 118
 119
 120   //---------------------------------------------------------------
 121   inline
 122   bool
 123   IsEqual(TFloat theLeft, TFloat theRight)
 124   {
 125     static TFloat EPS = 1.0E-3;
 126     if(fabs(theLeft) + fabs(theRight) > EPS)
 127       return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
 128     return true;
 129   }
 130
 131
 132   struct TShapeFun
 133   {
 134     class TFun
 135     {
 136       TFloatVector myFun;
 137       TInt myNbRef;
 138
 139     public:
 140
 141       void
 142       Init(TInt theNbGauss,
 143            TInt theNbRef)
 144       {
 145         myFun.resize(theNbGauss*theNbRef);
 146         myNbRef = theNbRef;
 147       }
 148
 149       TCFloatVecSlice
 150       GetFunSlice(TInt theGaussId) const
 151       {
 152         return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
 153       }
 154
 155       TFloatVecSlice
 156       GetFunSlice(TInt theGaussId)
 157       {
 158         return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
 159       }
 160     };
 161
 162     typedef TVector<TCCoordSlice> TCCoordSliceArr;
 163     typedef TVector<TCoordSlice> TCoordSliceArr;
 164
 165     TFloatVector myRefCoord;
 166     TInt myDim;
 167     TInt myNbRef;
 168
 169     TShapeFun(TInt theDim = 0, TInt theNbRef = 0):
 170       myDim(theDim),
 171       myNbRef(theNbRef),
 172       myRefCoord(theNbRef*theDim)
 173     {}
 174
 175     TInt
 176     GetNbRef() const
 177     {
 178       return myNbRef;
 179     }
 180
 181     TCCoordSlice
 182     GetCoord(TInt theRefId) const
 183     {
 184       return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
 185     }
 186
 187     TCoordSlice
 188     GetCoord(TInt theRefId)
 189     {
 190       return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
 191     }
 192
 193     void
 194     GetFun(const TShapeFun::TCCoordSliceArr& theRef,
 195            const TShapeFun::TCCoordSliceArr& theGauss,
 196            TFun& theFun) const
 197     {
 198       TInt aNbRef = theRef.size();
 199       TInt aNbGauss = theGauss.size();
 200       theFun.Init(aNbGauss,aNbRef);
 201     }
 202
 203     virtual
 204     void
 205     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 206             const TShapeFun::TCCoordSliceArr& theGauss,
 207             TFun& theFun) const = 0;
 208
 209     virtual
 210     bool
 211     IsSatisfy(const TShapeFun::TCCoordSliceArr& theRefCoord) const
 212     {
 213       TInt aNbRef = theRefCoord.size();
 214       TInt aNbRef2 = GetNbRef();
 215       INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
 216               "- aNbRef("<<aNbRef<<")"<<
 217               "; aNbRef2("<<aNbRef2<<")\n");
 218       bool anIsSatisfy = (aNbRef == aNbRef2);
 219 #ifndef _DEBUG_REF_COORDS_
 220       static TInt NB_REF_TO_CHECK = 2;
 221       aNbRef = NB_REF_TO_CHECK;
 222 #endif
 223       if(anIsSatisfy){
 224         for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 225           const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
 226           TCCoordSlice aCoord = GetCoord(aRefId);
 227           TInt aDim = aCoord.size();
 228           bool anIsEqual = false;
 229           for(TInt anId = 0; anId < aDim; anId++){
 230             anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
 231             if(!anIsEqual){
 232               anIsSatisfy = false;
 233               break;
 234             }
 235           }
 236           if(!anIsEqual){
 237             const TCCoordSlice& aCoord = theRefCoord[aRefId];
 238             INITMSG(MYDEBUG,aRefId + 1<<":{");
 239             TInt aDim = aCoord.size();
 240             for(TInt anId = 0; anId < aDim; anId++)
 241               ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
 242             ADDMSG(MYDEBUG,"}\n");
 243 #ifndef _DEBUG_
 244             return anIsSatisfy;
 245 #endif
 246           }
 247         }
 248       }
 249
 250       return anIsSatisfy;
 251     }
 252
 253     bool
 254     Eval(const TCellInfo& theCellInfo,
 255          const TNodeInfo& theNodeInfo,
 256          const TElemNum& theElemNum,
 257          const TShapeFun::TCCoordSliceArr& theRef,
 258          const TShapeFun::TCCoordSliceArr& theGauss,
 259          TGaussCoord& theGaussCoord,
 260          EModeSwitch theMode)
 261     {
 262       INITMSG(MYDEBUG,"TShapeFun::Eval"<<endl);
 263
 264       if(IsSatisfy(theRef)){
 265         const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
 266         TInt aDim = aMeshInfo->GetDim();
 267         TInt aNbGauss = theGauss.size();
 268
 269         bool anIsSubMesh = !theElemNum.empty();
 270         TInt aNbElem;
 271         if(anIsSubMesh)
 272           aNbElem = theElemNum.size();
 273         else
 274           aNbElem = theCellInfo.GetNbElem();
 275
 276         theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
 277
 278         TFun aFun;
 279         InitFun(theRef,theGauss,aFun);
 280         TInt aConnDim = theCellInfo.GetConnDim();
 281
 282         INITMSG(MYDEBUG,"aDim = "<<aDim<<
 283                 "; aNbGauss = "<<aNbGauss<<
 284                 "; aNbElem = "<<aNbElem<<
 285                 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
 286                 endl);
 287
 288         for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
 289           TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
 290           TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
 291           TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
 292
 293           for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 294             TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
 295             TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
 296
 297             for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
 298               TInt aNodeId = aConnSlice[aConnId] - 1;
 299               TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
 300
 301               for(TInt aDimId = 0; aDimId < aDim; aDimId++){
 302                 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
 303               }
 304             }
 305           }
 306         }
 307
 308 #ifdef _DEBUG_
 309         {
 310           INITMSG(MYVALUEDEBUG,"theGauss: ");
 311           for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 312             TCCoordSlice aCoordSlice = theGauss[aGaussId];
 313             ADDMSG(MYVALUEDEBUG,"{");
 314             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
 315               ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
 316             }
 317             ADDMSG(MYVALUEDEBUG,"} ");
 318           }
 319           ADDMSG(MYVALUEDEBUG,endl);
 320         }
 321         for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
 322           TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
 323           INITMSG(MYVALUEDEBUG,"");
 324           for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 325             TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
 326             ADDMSG(MYVALUEDEBUG,"{");
 327             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
 328               ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
 329             }
 330             ADDMSG(MYVALUEDEBUG,"} ");
 331           }
 332           ADDMSG(MYVALUEDEBUG,endl);
 333         }
 334 #endif
 335         return true;
 336       }
 337
 338       return false;
 339     }
 340   };
 341
 342
 343   //---------------------------------------------------------------
 344   struct TSeg2a: TShapeFun
 345   {
 346     TSeg2a():
 347       TShapeFun(1,2)
 348     {
 349       TInt aNbRef = GetNbRef();
 350       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 351         TCoordSlice aCoord = GetCoord(aRefId);
 352         switch(aRefId){
 353         case  0: aCoord[0] = -1.0; break;
 354         case  1: aCoord[0] =  1.0; break;
 355         }
 356       }
 357     }
 358
 359     virtual
 360     void
 361     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 362             const TShapeFun::TCCoordSliceArr& theGauss,
 363             TFun& theFun) const
 364     {
 365       GetFun(theRef,theGauss,theFun);
 366
 367       TInt aNbGauss = theGauss.size();
 368       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 369         const TCCoordSlice& aCoord = theGauss[aGaussId];
 370         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 371
 372         aSlice[0] = 0.5*(1.0 - aCoord[0]);
 373         aSlice[1] = 0.5*(1.0 + aCoord[0]);
 374       }
 375     }
 376   };
 377
 378
 379   //---------------------------------------------------------------
 380   struct TSeg3a: TShapeFun
 381   {
 382     TSeg3a():
 383       TShapeFun(1,3)
 384     {
 385       TInt aNbRef = GetNbRef();
 386       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 387         TCoordSlice aCoord = GetCoord(aRefId);
 388         switch(aRefId){
 389         case  0: aCoord[0] = -1.0; break;
 390         case  1: aCoord[0] =  1.0; break;
 391         case  2: aCoord[0] =  0.0; break;
 392         }
 393       }
 394     }
 395
 396     virtual
 397     void
 398     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 399             const TShapeFun::TCCoordSliceArr& theGauss,
 400             TFun& theFun) const
 401     {
 402       GetFun(theRef,theGauss,theFun);
 403
 404       TInt aNbGauss = theGauss.size();
 405       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 406         const TCCoordSlice& aCoord = theGauss[aGaussId];
 407         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 408
 409         aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
 410         aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
 411         aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
 412       }
 413     }
 414   };
 415
 416
 417   //---------------------------------------------------------------
 418   struct TTria3a: TShapeFun
 419   {
 420     TTria3a():
 421       TShapeFun(2,3)
 422     {
 423       TInt aNbRef = GetNbRef();
 424       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 425         TCoordSlice aCoord = GetCoord(aRefId);
 426         switch(aRefId){
 427         case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 428         case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 429         case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 430         }
 431       }
 432     }
 433
 434     virtual
 435     void
 436     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 437             const TShapeFun::TCCoordSliceArr& theGauss,
 438             TFun& theFun) const
 439     {
 440       GetFun(theRef,theGauss,theFun);
 441
 442       TInt aNbGauss = theGauss.size();
 443       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 444         const TCCoordSlice& aCoord = theGauss[aGaussId];
 445         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 446
 447         aSlice[0] = 0.5*(1.0 + aCoord[1]);
 448         aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
 449         aSlice[2] = 0.5*(1.0 + aCoord[0]);
 450       }
 451     }
 452   };
 453
 454
 455   //---------------------------------------------------------------
 456   struct TTria6a: TShapeFun
 457   {
 458     TTria6a():
 459       TShapeFun(2,6)
 460     {
 461       TInt aNbRef = GetNbRef();
 462       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 463         TCoordSlice aCoord = GetCoord(aRefId);
 464         switch(aRefId){
 465         case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 466         case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 467         case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 468
 469         case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 470         case  4: aCoord[0] =  0.0;  aCoord[1] = -1.0; break;
 471         case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 472         }
 473       }
 474     }
 475
 476     virtual
 477     void
 478     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 479             const TShapeFun::TCCoordSliceArr& theGauss,
 480             TFun& theFun) const
 481     {
 482       GetFun(theRef,theGauss,theFun);
 483
 484       TInt aNbGauss = theGauss.size();
 485       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 486         const TCCoordSlice& aCoord = theGauss[aGaussId];
 487         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 488
 489         aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
 490         aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
 491         aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
 492
 493         aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
 494         aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
 495         aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
 496       }
 497     }
 498   };
 499
 500
 501   //---------------------------------------------------------------
 502   struct TTria3b: TShapeFun
 503   {
 504     TTria3b():
 505       TShapeFun(2,3)
 506     {
 507       TInt aNbRef = GetNbRef();
 508       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 509         TCoordSlice aCoord = GetCoord(aRefId);
 510         switch(aRefId){
 511         case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 512         case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 513         case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 514         }
 515       }
 516     }
 517
 518     virtual
 519     void
 520     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 521             const TShapeFun::TCCoordSliceArr& theGauss,
 522             TFun& theFun) const
 523     {
 524       GetFun(theRef,theGauss,theFun);
 525
 526       TInt aNbGauss = theGauss.size();
 527       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 528         const TCCoordSlice& aCoord = theGauss[aGaussId];
 529         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 530
 531         aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
 532         aSlice[1] = aCoord[0];
 533         aSlice[2] = aCoord[1];
 534       }
 535     }
 536   };
 537
 538
 539   //---------------------------------------------------------------
 540   struct TTria6b: TShapeFun
 541   {
 542     TTria6b():
 543       TShapeFun(2,6)
 544     {
 545       TInt aNbRef = GetNbRef();
 546       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 547         TCoordSlice aCoord = GetCoord(aRefId);
 548         switch(aRefId){
 549         case  0: aCoord[0] =  0.0;  aCoord[1] =  0.0; break;
 550         case  1: aCoord[0] =  1.0;  aCoord[1] =  0.0; break;
 551         case  2: aCoord[0] =  0.0;  aCoord[1] =  1.0; break;
 552
 553         case  3: aCoord[0] =  0.5;  aCoord[1] =  0.0; break;
 554         case  4: aCoord[0] =  0.5;  aCoord[1] =  0.5; break;
 555         case  5: aCoord[0] =  0.0;  aCoord[1] =  0.5; break;
 556         }
 557       }
 558     }
 559
 560     virtual
 561     void
 562     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 563             const TShapeFun::TCCoordSliceArr& theGauss,
 564             TFun& theFun) const
 565     {
 566       GetFun(theRef,theGauss,theFun);
 567
 568       TInt aNbGauss = theGauss.size();
 569       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 570         const TCCoordSlice& aCoord = theGauss[aGaussId];
 571         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 572
 573         aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
 574         aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
 575         aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
 576
 577         aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
 578         aSlice[4] = 4.0*aCoord[0]*aCoord[1];
 579         aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
 580       }
 581     }
 582   };
 583
 584
 585   //---------------------------------------------------------------
 586   struct TQuad4a: TShapeFun
 587   {
 588     TQuad4a():
 589       TShapeFun(2,4)
 590     {
 591       TInt aNbRef = GetNbRef();
 592       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 593         TCoordSlice aCoord = GetCoord(aRefId);
 594         switch(aRefId){
 595         case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0; break;
 596         case  1: aCoord[0] = -1.0;  aCoord[1] = -1.0; break;
 597         case  2: aCoord[0] =  1.0;  aCoord[1] = -1.0; break;
 598         case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0; break;
 599         }
 600       }
 601     }
 602
 603     virtual
 604     void
 605     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 606             const TShapeFun::TCCoordSliceArr& theGauss,
 607             TFun& theFun) const
 608     {
 609       GetFun(theRef,theGauss,theFun);
 610
 611       TInt aNbGauss = theGauss.size();
 612       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 613         const TCCoordSlice& aCoord = theGauss[aGaussId];
 614         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 615
 616         aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
 617         aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
 618         aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
 619         aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
 620       }
 621     }
 622   };
 623
 624
 625   //---------------------------------------------------------------
 626   struct TTetra4a: TShapeFun
 627   {
 628     TTetra4a():
 629       TShapeFun(3,4)
 630     {
 631       TInt aNbRef = GetNbRef();
 632       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 633         TCoordSlice aCoord = GetCoord(aRefId);
 634         switch(aRefId){
 635         case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 636         case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 637         case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 638         case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 639         }
 640       }
 641     }
 642
 643     virtual
 644     void
 645     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 646             const TShapeFun::TCCoordSliceArr& theGauss,
 647             TFun& theFun) const
 648     {
 649       GetFun(theRef,theGauss,theFun);
 650
 651       TInt aNbGauss = theGauss.size();
 652       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 653         const TCCoordSlice& aCoord = theGauss[aGaussId];
 654         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 655
 656         aSlice[0] = aCoord[1];
 657         aSlice[1] = aCoord[2];
 658         aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
 659         aSlice[3] = aCoord[0];
 660       }
 661     }
 662   };
 663
 664
 665   //---------------------------------------------------------------
 666   struct TTetra10a: TShapeFun
 667   {
 668     TTetra10a():
 669       TShapeFun(3,10)
 670     {
 671       TInt aNbRef = GetNbRef();
 672       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 673         TCoordSlice aCoord = GetCoord(aRefId);
 674         switch(aRefId){
 675         case  0: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 676         case  1: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 677         case  2: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 678         case  3: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 679
 680         case  4: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 681         case  5: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 682         case  6: aCoord[0] =  0.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 683
 684         case  7: aCoord[0] =  0.5;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 685         case  8: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 686         case  9: aCoord[0] =  0.5;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 687         }
 688       }
 689     }
 690
 691     virtual
 692     void
 693     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 694             const TShapeFun::TCCoordSliceArr& theGauss,
 695             TFun& theFun) const
 696     {
 697       GetFun(theRef,theGauss,theFun);
 698
 699       TInt aNbGauss = theGauss.size();
 700       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 701         const TCCoordSlice& aCoord = theGauss[aGaussId];
 702         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 703
 704         aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
 705         aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
 706         aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
 707         aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
 708
 709         aSlice[4] = 4.0*aCoord[1]*aCoord[2];
 710         aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 711         aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 712
 713         aSlice[7] = 4.0*aCoord[0]*aCoord[1];
 714         aSlice[8] = 4.0*aCoord[0]*aCoord[2];
 715         aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 716       }
 717     }
 718   };
 719
 720
 721   //---------------------------------------------------------------
 722   struct THexa8a: TShapeFun
 723   {
 724     THexa8a():
 725       TShapeFun(3,8)
 726     {
 727       TInt aNbRef = GetNbRef();
 728       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 729         TCoordSlice aCoord = GetCoord(aRefId);
 730         switch(aRefId){
 731         case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 732         case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 733         case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 734         case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 735         case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 736         case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 737         case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 738         case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 739         }
 740       }
 741     }
 742
 743     virtual
 744     void
 745     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 746             const TShapeFun::TCCoordSliceArr& theGauss,
 747             TFun& theFun) const
 748     {
 749       GetFun(theRef,theGauss,theFun);
 750
 751       TInt aNbGauss = theGauss.size();
 752       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 753         const TCCoordSlice& aCoord = theGauss[aGaussId];
 754         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 755
 756         aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
 757         aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
 758         aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
 759         aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
 760
 761         aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
 762         aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
 763         aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
 764         aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
 765       }
 766     }
 767   };
 768
 769
 770   struct THexa20a: TShapeFun
 771   {
 772     THexa20a(TInt theDim = 3, TInt theNbRef = 20):
 773       TShapeFun(theDim,theNbRef)
 774     {
 775       TInt aNbRef = myRefCoord.size();
 776       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 777         TCoordSlice aCoord = GetCoord(aRefId);
 778         switch(aRefId){
 779         case  0: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 780         case  1: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 781         case  2: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 782         case  3: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 783         case  4: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 784         case  5: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 785         case  6: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 786         case  7: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 787
 788         case  8: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] = -1.0; break;
 789         case  9: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 790         case 10: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] = -1.0; break;
 791         case 11: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 792         case 12: aCoord[0] = -1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 793         case 13: aCoord[0] =  1.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 794         case 14: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 795         case 15: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 796         case 16: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  1.0; break;
 797         case 17: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 798         case 18: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  1.0; break;
 799         case 19: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 800         }
 801       }
 802     }
 803
 804     virtual
 805     void
 806     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 807             const TShapeFun::TCCoordSliceArr& theGauss,
 808             TFun& theFun) const
 809     {
 810       GetFun(theRef,theGauss,theFun);
 811
 812       TInt aNbGauss = theGauss.size();
 813       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 814         const TCCoordSlice& aCoord = theGauss[aGaussId];
 815         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 816
 817         aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
 818           (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
 819         aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
 820           (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
 821         aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
 822           (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
 823         aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
 824           (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
 825         aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
 826           (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
 827         aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
 828           (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
 829         aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
 830           (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
 831         aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
 832           (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
 833
 834         aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
 835         aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
 836         aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
 837         aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
 838         aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
 839         aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
 840         aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
 841         aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
 842         aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
 843         aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
 844         aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
 845         aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
 846       }
 847     }
 848   };
 849
 850
 851   //---------------------------------------------------------------
 852   struct THexa27a: THexa20a
 853   {
 854     THexa27a():
 855       THexa20a(3,27)
 856     {
 857       TInt aNbRef = myRefCoord.size();
 858       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 859         TCoordSlice aCoord = GetCoord(aRefId);
 860         switch(aRefId){
 861         case 20: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] = -1.0; break;
 862         case 21: aCoord[0] =  0.0;  aCoord[1] = -1.0;  aCoord[2] =  0.0; break;
 863         case 22: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 864         case 23: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 865         case 24: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 866         case 25: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 867         case 26: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 868         }
 869       }
 870     }
 871
 872     virtual
 873     void
 874     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 875             const TShapeFun::TCCoordSliceArr& theGauss,
 876             TFun& theFun) const
 877     {
 878       GetFun(theRef,theGauss,theFun);
 879
 880       TInt aNbGauss = theGauss.size();
 881       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 882         const TCCoordSlice& aCoord = theGauss[aGaussId];
 883         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 884
 885         aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
 886         aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
 887         aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
 888         aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
 889         aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
 890         aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
 891         aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
 892         aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
 893
 894         aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
 895         aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
 896         aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
 897         aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
 898         aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
 899         aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
 900         aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
 901         aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
 902         aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
 903         aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
 904         aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
 905         aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
 906         aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
 907         aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
 908         aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
 909         aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
 910         aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
 911         aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
 912         aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
 913       }
 914     }
 915   };
 916
 917
 918   //---------------------------------------------------------------
 919   struct TPenta6a: TShapeFun
 920   {
 921     TPenta6a():
 922       TShapeFun(3,6)
 923     {
 924       TInt aNbRef = myRefCoord.size();
 925       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 926         TCoordSlice aCoord = GetCoord(aRefId);
 927         switch(aRefId){
 928         case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 929         case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 930         case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 931         case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 932         case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 933         case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 934         }
 935       }
 936     }
 937
 938     virtual
 939     void
 940     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 941             const TShapeFun::TCCoordSliceArr& theGauss,
 942             TFun& theFun) const
 943     {
 944       GetFun(theRef,theGauss,theFun);
 945
 946       TInt aNbGauss = theGauss.size();
 947       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
 948         const TCCoordSlice& aCoord = theGauss[aGaussId];
 949         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
 950
 951         aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
 952         aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
 953         aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
 954
 955         aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
 956         aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
 957         aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
 958       }
 959     }
 960   };
 961
 962
 963   //---------------------------------------------------------------
 964   struct TPenta15a: TShapeFun
 965   {
 966     TPenta15a():
 967       TShapeFun(3,15)
 968     {
 969       TInt aNbRef = myRefCoord.size();
 970       for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
 971         TCoordSlice aCoord = GetCoord(aRefId);
 972         switch(aRefId){
 973         case  0: aCoord[0] = -1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 974         case  1: aCoord[0] = -1.0;  aCoord[1] = -0.0;  aCoord[2] =  1.0; break;
 975         case  2: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 976         case  3: aCoord[0] =  1.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 977         case  4: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 978         case  5: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 979
 980         case  6: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 981         case  7: aCoord[0] = -1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 982         case  8: aCoord[0] = -1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 983         case  9: aCoord[0] =  0.0;  aCoord[1] =  1.0;  aCoord[2] =  0.0; break;
 984         case 10: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  1.0; break;
 985         case 11: aCoord[0] =  0.0;  aCoord[1] =  0.0;  aCoord[2] =  0.0; break;
 986         case 12: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.5; break;
 987         case 13: aCoord[0] =  1.0;  aCoord[1] =  0.0;  aCoord[2] =  0.5; break;
 988         case 14: aCoord[0] =  1.0;  aCoord[1] =  0.5;  aCoord[2] =  0.0; break;
 989         }
 990       }
 991     }
 992
 993     virtual
 994     void
 995     InitFun(const TShapeFun::TCCoordSliceArr& theRef,
 996             const TShapeFun::TCCoordSliceArr& theGauss,
 997             TFun& theFun) const
 998     {
 999       GetFun(theRef,theGauss,theFun);
1000
1001       TInt aNbGauss = theGauss.size();
1002       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1003         const TCCoordSlice& aCoord = theGauss[aGaussId];
1004         TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1005
1006         aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1007         aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1008         aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1009
1010         aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1011         aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1012         aSlice[5] = -0.5*(aCoord[0] + 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1013
1014         aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1015         aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1016         aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1017
1018         aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1019         aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1020         aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1021
1022         aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1023         aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1024         aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1025       }
1026     }
1027   };
1028
1029
1030   //---------------------------------------------------------------
1031   bool
1032   GetGaussCoord3D(const TGaussInfo& theGaussInfo,
1033                   const TCellInfo& theCellInfo,
1034                   const TNodeInfo& theNodeInfo,
1035                   TGaussCoord& theGaussCoord,
1036                   const TElemNum& theElemNum,
1037                   EModeSwitch theMode)
1038   {
1039     INITMSG(MYDEBUG,"GetGaussCoord3D\n");
1040
1041     if(theGaussInfo.myGeom == theCellInfo.myGeom){
1042       EGeometrieElement aGeom = theGaussInfo.myGeom;
1043
1044       TInt aNbRef = theGaussInfo.GetNbRef();
1045       TCCoordSliceArr aRefSlice(aNbRef);
1046       for(TInt anId = 0; anId < aNbRef; anId++)
1047         aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
1048
1049       TInt aNbGauss = theGaussInfo.GetNbGauss();
1050       TCCoordSliceArr aGaussSlice(aNbGauss);
1051       for(TInt anId = 0; anId < aNbGauss; anId++)
1052         aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
1053
1054       switch(aGeom){
1055       case eSEG2: {
1056         INITMSG(MYDEBUG,"eSEG2"<<endl);
1057
1058         if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1059           return true;
1060
1061         break;
1062       }
1063       case eSEG3: {
1064         INITMSG(MYDEBUG,"eSEG3"<<endl);
1065
1066         if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1067           return true;
1068
1069         break;
1070       }
1071       case eTRIA3: {
1072         INITMSG(MYDEBUG,"eTRIA3"<<endl);
1073
1074         if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1075           return true;
1076
1077         if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1078           return true;
1079
1080         break;
1081       }
1082       case eTRIA6: {
1083         INITMSG(MYDEBUG,"eTRIA6"<<endl);
1084
1085         if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1086           return true;
1087
1088         if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1089           return true;
1090
1091         break;
1092       }
1093       case eQUAD4: {
1094         INITMSG(MYDEBUG,"eQUAD4"<<endl);
1095
1096         if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1097           return true;
1098
1099         break;
1100       }
1101       case eQUAD8: {
1102         INITMSG(MYDEBUG,"eQUAD8"<<endl);
1103         break;
1104       }
1105       case eTETRA4: {
1106         INITMSG(MYDEBUG,"eTETRA4"<<endl);
1107
1108         if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1109           return true;
1110
1111         break;
1112       }
1113       case ePYRA5: {
1114         INITMSG(MYDEBUG,"ePYRA5"<<endl);
1115         break;
1116       }
1117       case ePENTA6: {
1118         INITMSG(MYDEBUG,"ePENTA6"<<endl);
1119
1120         if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1121           return true;
1122
1123         break;
1124       }
1125       case eHEXA8: {
1126         INITMSG(MYDEBUG,"eHEXA8"<<endl);
1127
1128         if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1129           return true;
1130
1131         break;
1132       }
1133       case eTETRA10: {
1134         INITMSG(MYDEBUG,"eTETRA10"<<endl);
1135
1136         if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1137           return true;
1138
1139         break;
1140       }
1141       case ePYRA13: {
1142         INITMSG(MYDEBUG,"ePYRA13"<<endl);
1143         break;
1144       }
1145       case ePENTA15: {
1146         INITMSG(MYDEBUG,"ePENTA15"<<endl);
1147
1148         if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1149           return true;
1150
1151         break;
1152       }
1153       case eHEXA20: {
1154         INITMSG(MYDEBUG,"eHEXA20"<<endl);
1155
1156         if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1157           return true;
1158
1159         break;
1160       }
1161       default:
1162         INITMSG(MYDEBUG,"eNONE"<<endl);
1163         return false;
1164       }
1165     }
1166
1167     return false;
1168   }
1169
1170   //---------------------------------------------------------------
1171   bool
1172   GetBaryCenter(const TCellInfo& theCellInfo,
1173                 const TNodeInfo& theNodeInfo,
1174                 TGaussCoord& theGaussCoord,
1175                 const TElemNum& theElemNum,
1176                 EModeSwitch theMode)
1177   {
1178     INITMSG(MYDEBUG,"GetBaryCenter\n");
1179     const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
1180     TInt aDim = aMeshInfo->GetDim();
1181     static TInt aNbGauss = 1;
1182
1183     bool anIsSubMesh = !theElemNum.empty();
1184     TInt aNbElem;
1185     if(anIsSubMesh)
1186       aNbElem = theElemNum.size();
1187     else
1188       aNbElem = theCellInfo.GetNbElem();
1189
1190     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
1191
1192     INITMSGA(MYDEBUG,0,
1193              "- aDim = "<<aDim<<
1194              "; aNbGauss = "<<aNbGauss<<
1195              "; aNbElem = "<<aNbElem<<
1196              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
1197              endl);
1198
1199     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1200       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
1201       TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
1202       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1203
1204       TInt aConnDim = aConnSlice.size();
1205
1206       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1207         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
1208
1209         for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
1210           TInt aNodeId = aConnSlice[aConnId] - 1;
1211           TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
1212
1213           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1214             aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
1215           }
1216         }
1217
1218         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1219           aGaussCoordSlice[aDimId] /= aConnDim;
1220         }
1221       }
1222     }
1223
1224 #ifdef _DEBUG_
1225     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1226       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1227       INITMSG(MYVALUEDEBUG,"");
1228       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1229         TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
1230         ADDMSG(MYVALUEDEBUG,"{");
1231         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1232           ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
1233         }
1234         ADDMSG(MYVALUEDEBUG,"} ");
1235       }
1236       ADDMSG(MYVALUEDEBUG,endl);
1237     }
1238 #endif
1239
1240     return true;
1241   }
1242
1243
1244   //---------------------------------------------------------------
1245   bool
1246   GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
1247                 const TNodeInfo& theNodeInfo,
1248                 TGaussCoord& theGaussCoord,
1249                 const TElemNum& theElemNum,
1250                 EModeSwitch theMode)
1251   {
1252     INITMSG(MYDEBUG,"GetBaryCenter\n");
1253     const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
1254     TInt aDim = aMeshInfo->GetDim();
1255     static TInt aNbGauss = 1;
1256
1257     bool anIsSubMesh = !theElemNum.empty();
1258     TInt aNbElem;
1259     if(anIsSubMesh)
1260       aNbElem = theElemNum.size();
1261     else
1262       aNbElem = thePolygoneInfo.GetNbElem();
1263
1264     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
1265
1266     INITMSGA(MYDEBUG,0,
1267              "- aDim = "<<aDim<<
1268              "; aNbGauss = "<<aNbGauss<<
1269              "; aNbElem = "<<aNbElem<<
1270              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
1271              endl);
1272
1273     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1274       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
1275
1276       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1277       TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
1278       TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
1279       TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
1280
1281       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1282         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
1283
1284         for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
1285           TInt aNodeId = aConnSlice[aConnId] - 1;
1286           TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
1287
1288           for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1289             aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
1290           }
1291         }
1292
1293         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1294           aGaussCoordSlice[aDimId] /= aNbNodes;
1295         }
1296       }
1297     }
1298
1299     return true;
1300   }
1301
1302
1303   //---------------------------------------------------------------
1304   bool
1305   GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
1306                 const TNodeInfo& theNodeInfo,
1307                 TGaussCoord& theGaussCoord,
1308                 const TElemNum& theElemNum,
1309                 EModeSwitch theMode)
1310   {
1311     INITMSG(MYDEBUG,"GetBaryCenter\n");
1312     const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
1313     TInt aDim = aMeshInfo->GetDim();
1314     static TInt aNbGauss = 1;
1315
1316     bool anIsSubMesh = !theElemNum.empty();
1317     TInt aNbElem;
1318     if(anIsSubMesh)
1319       aNbElem = theElemNum.size();
1320     else
1321       aNbElem = thePolyedreInfo.GetNbElem();
1322
1323     theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
1324
1325     INITMSGA(MYDEBUG,0,
1326              "- aDim = "<<aDim<<
1327              "; aNbGauss = "<<aNbGauss<<
1328              "; aNbElem = "<<aNbElem<<
1329              "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
1330              endl);
1331
1332     for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1333       TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
1334
1335       TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1336       TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
1337       TInt aNbFaces = aConnSliceArr.size();
1338
1339       TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
1340
1341       for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1342         TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
1343
1344         for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
1345           TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
1346           TInt aNbConn = aConnSlice.size();
1347           for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
1348             TInt aNodeId = aConnSlice[aConnId] - 1;
1349             TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
1350
1351             for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1352               aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
1353             }
1354           }
1355         }
1356         for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1357           aGaussCoordSlice[aDimId] /= aNbNodes;
1358         }
1359       }
1360     }
1361
1362     return true;
1363   }
1364 }