smesh/src/MEDWrapper/MED_GaussUtils.cxx
2019-02-14 15:33:05 +03:00

2107 lines
81 KiB
C++

// Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
#include "MED_GaussUtils.hxx"
#include "MED_Utilities.hxx"
#ifdef _DEBUG_
static int MYDEBUG = 0;
static int MYVALUEDEBUG = 0;
#else
static int MYDEBUG = 0;
static int MYVALUEDEBUG = 0;
#endif
namespace MED
{
//---------------------------------------------------------------
TGaussCoord
::TGaussCoord():
TModeSwitchInfo(eFULL_INTERLACE),
myNbElem(0),
myNbGauss(0),
myDim(0),
myGaussStep(0)
{
}
void
TGaussCoord
::Init(TInt theNbElem,
TInt theNbGauss,
TInt theDim,
EModeSwitch theMode)
{
myModeSwitch = theMode;
myNbElem = theNbElem;
myNbGauss = theNbGauss;
myDim = theDim;
myGaussStep = myNbGauss*myDim;
myGaussCoord.resize(theNbElem*myGaussStep);
}
TInt
TGaussCoord
::GetNbElem() const
{
return myNbElem;
}
TInt
TGaussCoord
::GetNbGauss() const
{
return myNbGauss;
}
TInt
TGaussCoord
::GetDim() const
{
return myDim;
}
unsigned char*
TGaussCoord
::GetValuePtr()
{
return (unsigned char*)&(myGaussCoord[0]);
}
TCCoordSliceArr
TGaussCoord
::GetCoordSliceArr(TInt theElemId) const
{
TCCoordSliceArr aCoordSliceArr(myNbGauss);
if(GetModeSwitch() == eFULL_INTERLACE){
TInt anId = theElemId*myGaussStep;
for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
aCoordSliceArr[anGaussId] =
TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
anId += myDim;
}
}
else{
for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
aCoordSliceArr[anGaussId] =
TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
}
}
return aCoordSliceArr;
}
TCoordSliceArr
TGaussCoord
::GetCoordSliceArr(TInt theElemId)
{
TCoordSliceArr aCoordSliceArr(myNbGauss);
if(GetModeSwitch() == eFULL_INTERLACE){
TInt anId = theElemId*myGaussStep;
for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
aCoordSliceArr[anGaussId] =
TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
anId += myDim;
}
}
else{
for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
aCoordSliceArr[anGaussId] =
TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
}
}
return aCoordSliceArr;
}
//---------------------------------------------------------------
inline
bool
IsEqual(TFloat theLeft, TFloat theRight)
{
static TFloat EPS = 1.0E-3;
if(fabs(theLeft) + fabs(theRight) > EPS)
return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
return true;
}
//---------------------------------------------------------------
class TShapeFun::TFun
{
TFloatVector myFun;
TInt myNbRef;
public:
void
Init(TInt theNbGauss,
TInt theNbRef)
{
myFun.resize(theNbGauss*theNbRef);
myNbRef = theNbRef;
}
TCFloatVecSlice
GetFunSlice(TInt theGaussId) const
{
return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
}
TFloatVecSlice
GetFunSlice(TInt theGaussId)
{
return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
}
};
//---------------------------------------------------------------
TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
myRefCoord(theNbRef*theDim),
myDim(theDim),
myNbRef(theNbRef)
{}
TCCoordSlice
TShapeFun::GetCoord(TInt theRefId) const
{
return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
}
TCoordSlice
TShapeFun::GetCoord(TInt theRefId)
{
return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
}
void
TShapeFun::GetFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
TInt aNbRef = theRef.size();
TInt aNbGauss = theGauss.size();
theFun.Init(aNbGauss,aNbRef);
}
bool
TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
{
TInt aNbRef = theRefCoord.size();
TInt aNbRef2 = GetNbRef();
INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
"- aNbRef("<<aNbRef<<")"<<
"; aNbRef2("<<aNbRef2<<")\n");
bool anIsSatisfy = (aNbRef == aNbRef2);
if(anIsSatisfy){
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
TCCoordSlice aCoord = GetCoord(aRefId);
TInt aDim = aCoord.size();
bool anIsEqual = false;
for(TInt anId = 0; anId < aDim; anId++){
anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
if(!anIsEqual){
anIsSatisfy = false;
break;
}
}
if(!anIsEqual){
#ifdef _DEBUG_
TCCoordSlice aCoord = GetCoord(aRefId);
INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
TInt aDim = aCoord.size();
for(TInt anId = 0; anId < aDim; anId++)
ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
for(TInt anId = 0; anId < aDim; anId++)
ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
ADDMSG(MYDEBUG,"}\n");
#endif
#ifndef _DEBUG_
BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
return anIsSatisfy;
#endif
}
}
}
BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
return anIsSatisfy;
}
bool
TShapeFun::Eval(const TCellInfo& theCellInfo,
const TNodeInfo& theNodeInfo,
const TElemNum& theElemNum,
const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TGaussCoord& theGaussCoord,
EModeSwitch theMode)
{
INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);
if(IsSatisfy(theRef)){
const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
TInt aDim = aMeshInfo->GetDim();
TInt aNbGauss = theGauss.size();
bool anIsSubMesh = !theElemNum.empty();
TInt aNbElem;
if(anIsSubMesh)
aNbElem = theElemNum.size();
else
aNbElem = theCellInfo.GetNbElem();
theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
TFun aFun;
InitFun(theRef,theGauss,aFun);
TInt aConnDim = theCellInfo.GetConnDim();
INITMSG(MYDEBUG,"aDim = "<<aDim<<
"; aNbGauss = "<<aNbGauss<<
"; aNbElem = "<<aNbElem<<
"; aNbNodes = "<<theNodeInfo.GetNbElem()<<
std::endl);
for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
TInt aNodeId = aConnSlice[aConnId] - 1;
TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
}
}
}
}
#ifdef _DEBUG_
{
INITMSG(MYVALUEDEBUG,"theGauss: ");
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCCoordSlice aCoordSlice = theGauss[aGaussId];
ADDMSG(MYVALUEDEBUG,"{");
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
}
ADDMSG(MYVALUEDEBUG,"} ");
}
ADDMSG(MYVALUEDEBUG,std::endl);
}
for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
INITMSG(MYVALUEDEBUG,"");
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
ADDMSG(MYVALUEDEBUG,"{");
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
}
ADDMSG(MYVALUEDEBUG,"} ");
}
ADDMSG(MYVALUEDEBUG,std::endl);
}
#endif
return true;
}
return false;
}
//---------------------------------------------------------------
TSeg2a::TSeg2a():TShapeFun(1,2)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; break;
case 1: aCoord[0] = 1.0; break;
}
}
}
void
TSeg2a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*(1.0 - aCoord[0]);
aSlice[1] = 0.5*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TSeg3a::TSeg3a():TShapeFun(1,3)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; break;
case 1: aCoord[0] = 1.0; break;
case 2: aCoord[0] = 0.0; break;
}
}
}
void
TSeg3a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
}
}
//---------------------------------------------------------------
TTria3a::TTria3a():
TShapeFun(2,3)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
}
}
}
void
TTria3a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*(1.0 + aCoord[1]);
aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
aSlice[2] = 0.5*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TTria6a::TTria6a():TShapeFun(2,6)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
}
}
}
void
TTria6a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
}
}
//---------------------------------------------------------------
TTria3b::TTria3b():
TShapeFun(2,3)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
}
}
}
void
TTria3b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
aSlice[1] = aCoord[0];
aSlice[2] = aCoord[1];
}
}
//---------------------------------------------------------------
TTria6b::TTria6b():
TShapeFun(2,6)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
case 3: aCoord[0] = 0.5; aCoord[1] = 0.0; break;
case 4: aCoord[0] = 0.5; aCoord[1] = 0.5; break;
case 5: aCoord[0] = 0.0; aCoord[1] = 0.5; break;
}
}
}
void
TTria6b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
aSlice[4] = 4.0*aCoord[0]*aCoord[1];
aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
}
}
//---------------------------------------------------------------
TQuad4a::TQuad4a():
TShapeFun(2,4)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
}
}
}
void
TQuad4a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
}
}
//---------------------------------------------------------------
TQuad8a::TQuad8a():
TShapeFun(2,8)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
}
}
}
void
TQuad8a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);
aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TQuad9a::TQuad9a():
TShapeFun(2,9)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
}
}
}
void
TQuad9a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
}
}
//---------------------------------------------------------------
TQuad4b::TQuad4b():
TShapeFun(2,4)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
}
}
}
void
TQuad4b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
}
}
//---------------------------------------------------------------
TQuad8b::TQuad8b():
TShapeFun(2,8)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
}
}
}
void
TQuad8b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);
aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);
//aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
//aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
//aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
//aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
}
}
//---------------------------------------------------------------
TQuad9b::TQuad9b():
TShapeFun(2,9)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
}
}
}
void
TQuad9b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
}
}
//---------------------------------------------------------------
TTetra4a::TTetra4a():
TShapeFun(3,4)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
TTetra4a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = aCoord[1];
aSlice[1] = aCoord[2];
aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
aSlice[3] = aCoord[0];
}
}
//---------------------------------------------------------------
TTetra10a::TTetra10a():
TShapeFun(3,10)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
TTetra10a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
aSlice[4] = 4.0*aCoord[1]*aCoord[2];
aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[7] = 4.0*aCoord[0]*aCoord[1];
aSlice[8] = 4.0*aCoord[0]*aCoord[2];
aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
}
}
//---------------------------------------------------------------
TTetra4b::TTetra4b():
TShapeFun(3,4)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
TTetra4b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = aCoord[1];
aSlice[2] = aCoord[2];
aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
aSlice[3] = aCoord[0];
}
}
//---------------------------------------------------------------
TTetra10b::TTetra10b():
TShapeFun(3,10)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
TTetra10b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
aSlice[6] = 4.0*aCoord[1]*aCoord[2];
aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[7] = 4.0*aCoord[0]*aCoord[1];
aSlice[9] = 4.0*aCoord[0]*aCoord[2];
aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
}
}
//---------------------------------------------------------------
THexa8a::THexa8a():
TShapeFun(3,8)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
}
}
}
void
THexa8a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
}
}
//---------------------------------------------------------------
THexa20a::THexa20a(TInt theDim, TInt theNbRef):
TShapeFun(theDim,theNbRef)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 8: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 9: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
case 10: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 11: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
case 12: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 13: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 14: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 15: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 16: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 17: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 18: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 19: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
}
}
}
void
THexa20a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
(-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
(-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
(-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
(-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
(-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
(-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
(-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
(-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
}
}
//---------------------------------------------------------------
THexa27a::THexa27a():
THexa20a(3,27)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 20: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
case 21: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 22: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 23: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 24: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 25: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 26: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
THexa27a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
}
}
//---------------------------------------------------------------
THexa8b::THexa8b():
TShapeFun(3,8)
{
TInt aNbRef = GetNbRef();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
}
}
}
void
THexa8b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
}
}
//---------------------------------------------------------------
THexa20b::THexa20b(TInt theDim, TInt theNbRef):
TShapeFun(theDim,theNbRef)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 11: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
case 8: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
case 16: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 19: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 18: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 17: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 15: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
case 14: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 13: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
case 12: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
}
}
}
void
THexa20b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
(-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
(-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
(-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
(-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
(-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
(-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
(-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
(-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
}
}
//---------------------------------------------------------------
TPenta6a::TPenta6a():
TShapeFun(3,6)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
TPenta6a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TPenta6b::TPenta6b():
TShapeFun(3,6)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
}
}
}
void
TPenta6b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TPenta15a::TPenta15a():
TShapeFun(3,15)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 10: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 11: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 12: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 13: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 14: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
}
}
}
void
TPenta15a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TPenta15b::TPenta15b():
TShapeFun(3,15)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 12: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 14: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 13: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 11: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 9: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
}
}
}
void
TPenta15b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
aSlice[9] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
}
}
//---------------------------------------------------------------
TPyra5a::TPyra5a():
TShapeFun(3,5)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
}
}
}
void
TPyra5a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
// APO & RNV:
// Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
// Seems shape function for ePYRA5 elements is:
// w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
// w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
// w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
// w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
// w4 = +Z;
aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[4] = aCoord[2];
}
}
//---------------------------------------------------------------
TPyra5b::TPyra5b():
TShapeFun(3,5)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
}
}
}
void
TPyra5b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
// APO & RNV:
// Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
// Seems shape function for ePYRA5 elements is:
// w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
// w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
// w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
// w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
// w4 = +Z;
aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
aSlice[4] = aCoord[2];
}
}
//---------------------------------------------------------------
TPyra13a::TPyra13a():
TShapeFun(3,13)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 5: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 6: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 7: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
case 8: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 10: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 12: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
}
}
}
void
TPyra13a::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[0] - 0.5)/(1.0 - aCoord[2]);
aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[1] - 0.5)/(1.0 - aCoord[2]);
aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
}
}
//---------------------------------------------------------------
TPyra13b::TPyra13b():
TShapeFun(3,13)
{
TInt aNbRef = myRefCoord.size();
for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
TCoordSlice aCoord = GetCoord(aRefId);
switch(aRefId){
case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
case 8: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 7: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
case 6: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
case 5: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 12: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
case 10: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
}
}
}
void
TPyra13b::InitFun(const TCCoordSliceArr& theRef,
const TCCoordSliceArr& theGauss,
TFun& theFun) const
{
GetFun(theRef,theGauss,theFun);
TInt aNbGauss = theGauss.size();
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
const TCCoordSlice& aCoord = theGauss[aGaussId];
TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[0] - 0.5)/(1.0 - aCoord[2]);
aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[1] - 0.5)/(1.0 - aCoord[2]);
aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
(1.0 - aCoord[2]);
}
}
//---------------------------------------------------------------
bool
GetGaussCoord3D(const TGaussInfo& theGaussInfo,
const TCellInfo& theCellInfo,
const TNodeInfo& theNodeInfo,
TGaussCoord& theGaussCoord,
const TElemNum& theElemNum,
EModeSwitch theMode)
{
INITMSG(MYDEBUG,"GetGaussCoord3D\n");
if(theGaussInfo.myGeom == theCellInfo.myGeom){
EGeometrieElement aGeom = theGaussInfo.myGeom;
TInt aNbRef = theGaussInfo.GetNbRef();
TCCoordSliceArr aRefSlice(aNbRef);
for(TInt anId = 0; anId < aNbRef; anId++)
aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
TInt aNbGauss = theGaussInfo.GetNbGauss();
TCCoordSliceArr aGaussSlice(aNbGauss);
for(TInt anId = 0; anId < aNbGauss; anId++)
aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
switch(aGeom){
case eSEG2: {
INITMSG(MYDEBUG,"eSEG2"<<std::endl);
if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eSEG3: {
INITMSG(MYDEBUG,"eSEG3"<<std::endl);
if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eTRIA3: {
INITMSG(MYDEBUG,"eTRIA3"<<std::endl);
if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eTRIA6: {
INITMSG(MYDEBUG,"eTRIA6"<<std::endl);
if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eQUAD4: {
INITMSG(MYDEBUG,"eQUAD4"<<std::endl);
if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eQUAD8: {
INITMSG(MYDEBUG,"eQUAD8"<<std::endl);
if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eQUAD9: {
INITMSG(MYDEBUG,"eQUAD9"<<std::endl);
if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eTETRA4: {
INITMSG(MYDEBUG,"eTETRA4"<<std::endl);
if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case ePYRA5: {
INITMSG(MYDEBUG,"ePYRA5"<<std::endl);
if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case ePENTA6: {
INITMSG(MYDEBUG,"ePENTA6"<<std::endl);
if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eHEXA8: {
INITMSG(MYDEBUG,"eHEXA8"<<std::endl);
if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eTETRA10: {
INITMSG(MYDEBUG,"eTETRA10"<<std::endl);
if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case ePYRA13: {
INITMSG(MYDEBUG,"ePYRA13"<<std::endl);
if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case ePENTA15: {
INITMSG(MYDEBUG,"ePENTA15"<<std::endl);
if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
case eHEXA20: {
INITMSG(MYDEBUG,"eHEXA20"<<std::endl);
if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
return true;
break;
}
default:
INITMSG(MYDEBUG,"eNONE"<<std::endl);
return false;
}
}
return false;
}
//---------------------------------------------------------------
bool
GetBaryCenter(const TCellInfo& theCellInfo,
const TNodeInfo& theNodeInfo,
TGaussCoord& theGaussCoord,
const TElemNum& theElemNum,
EModeSwitch theMode)
{
INITMSG(MYDEBUG,"GetBaryCenter\n");
const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
TInt aDim = aMeshInfo->GetDim();
static TInt aNbGauss = 1;
bool anIsSubMesh = !theElemNum.empty();
TInt aNbElem;
if(anIsSubMesh)
aNbElem = theElemNum.size();
else
aNbElem = theCellInfo.GetNbElem();
theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
TInt aConnDim = theCellInfo.GetConnDim();
INITMSGA(MYDEBUG,0,
"- aDim = "<<aDim<<
"; aNbGauss = "<<aNbGauss<<
"; aNbElem = "<<aNbElem<<
"; aNbNodes = "<<theNodeInfo.GetNbElem()<<
std::endl);
for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
TInt aNodeId = aConnSlice[aConnId] - 1;
TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
}
}
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] /= aConnDim;
}
}
}
#ifdef _DEBUG_
for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
INITMSG(MYVALUEDEBUG,"");
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
ADDMSG(MYVALUEDEBUG,"{");
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
}
ADDMSG(MYVALUEDEBUG,"} ");
}
ADDMSG(MYVALUEDEBUG,std::endl);
}
#endif
return true;
}
//---------------------------------------------------------------
bool
GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
const TNodeInfo& theNodeInfo,
TGaussCoord& theGaussCoord,
const TElemNum& theElemNum,
EModeSwitch theMode)
{
INITMSG(MYDEBUG,"GetBaryCenter\n");
const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
TInt aDim = aMeshInfo->GetDim();
static TInt aNbGauss = 1;
bool anIsSubMesh = !theElemNum.empty();
TInt aNbElem;
if(anIsSubMesh)
aNbElem = theElemNum.size();
else
aNbElem = thePolygoneInfo.GetNbElem();
theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
INITMSGA(MYDEBUG,0,
"- aDim = "<<aDim<<
"; aNbGauss = "<<aNbGauss<<
"; aNbElem = "<<aNbElem<<
"; aNbNodes = "<<theNodeInfo.GetNbElem()<<
std::endl);
for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
TInt aNodeId = aConnSlice[aConnId] - 1;
TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
}
}
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] /= aNbNodes;
}
}
}
return true;
}
//---------------------------------------------------------------
bool
GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
const TNodeInfo& theNodeInfo,
TGaussCoord& theGaussCoord,
const TElemNum& theElemNum,
EModeSwitch theMode)
{
INITMSG(MYDEBUG,"GetBaryCenter\n");
const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
TInt aDim = aMeshInfo->GetDim();
static TInt aNbGauss = 1;
bool anIsSubMesh = !theElemNum.empty();
TInt aNbElem;
if(anIsSubMesh)
aNbElem = theElemNum.size();
else
aNbElem = thePolyedreInfo.GetNbElem();
theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
INITMSGA(MYDEBUG,0,
"- aDim = "<<aDim<<
"; aNbGauss = "<<aNbGauss<<
"; aNbElem = "<<aNbElem<<
"; aNbNodes = "<<theNodeInfo.GetNbElem()<<
std::endl);
for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
TInt aNbFaces = aConnSliceArr.size();
TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
TInt aNbConn = aConnSlice.size();
for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
TInt aNodeId = aConnSlice[aConnId] - 1;
TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
}
}
}
for(TInt aDimId = 0; aDimId < aDim; aDimId++){
aGaussCoordSlice[aDimId] /= aNbNodes;
}
}
}
return true;
}
}