1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "MED_GaussUtils.hxx"
24 #include "MED_Utilities.hxx"
27 static int MYDEBUG = 0;
28 static int MYVALUEDEBUG = 0;
30 // static int MYDEBUG = 0;
31 // static int MYVALUEDEBUG = 0;
34 //#define _DEBUG_REF_COORDS_
38 //---------------------------------------------------------------
41 TModeSwitchInfo(eFULL_INTERLACE),
51 ::Init(TInt theNbElem,
56 myModeSwitch = theMode;
59 myNbGauss = theNbGauss;
62 myGaussStep = myNbGauss*myDim;
64 myGaussCoord.resize(theNbElem*myGaussStep);
93 return (unsigned char*)&(myGaussCoord[0]);
99 ::GetCoordSliceArr(TInt theElemId) const
101 TCCoordSliceArr aCoordSliceArr(myNbGauss);
102 if(GetModeSwitch() == eFULL_INTERLACE){
103 TInt anId = theElemId*myGaussStep;
104 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
105 aCoordSliceArr[anGaussId] =
106 TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
111 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
112 aCoordSliceArr[anGaussId] =
113 TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
116 return aCoordSliceArr;
122 ::GetCoordSliceArr(TInt theElemId)
124 TCoordSliceArr aCoordSliceArr(myNbGauss);
125 if(GetModeSwitch() == eFULL_INTERLACE){
126 TInt anId = theElemId*myGaussStep;
127 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
128 aCoordSliceArr[anGaussId] =
129 TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
134 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
135 aCoordSliceArr[anGaussId] =
136 TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
139 return aCoordSliceArr;
143 //---------------------------------------------------------------
146 IsEqual(TFloat theLeft, TFloat theRight)
148 static TFloat EPS = 1.0E-3;
149 if(fabs(theLeft) + fabs(theRight) > EPS)
150 return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
155 //---------------------------------------------------------------
156 class TShapeFun::TFun
164 Init(TInt theNbGauss,
167 myFun.resize(theNbGauss*theNbRef);
172 GetFunSlice(TInt theGaussId) const
174 return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
178 GetFunSlice(TInt theGaussId)
180 return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
184 //---------------------------------------------------------------
186 TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
187 myRefCoord(theNbRef*theDim),
193 TShapeFun::GetCoord(TInt theRefId) const
195 return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
199 TShapeFun::GetCoord(TInt theRefId)
201 return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
205 TShapeFun::GetFun(const TCCoordSliceArr& theRef,
206 const TCCoordSliceArr& theGauss,
209 TInt aNbRef = theRef.size();
210 TInt aNbGauss = theGauss.size();
211 theFun.Init(aNbGauss,aNbRef);
215 TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
217 TInt aNbRef = theRefCoord.size();
218 TInt aNbRef2 = GetNbRef();
219 INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
220 "- aNbRef("<<aNbRef<<")"<<
221 "; aNbRef2("<<aNbRef2<<")\n");
222 bool anIsSatisfy = (aNbRef == aNbRef2);
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]);
238 TCCoordSlice aCoord = GetCoord(aRefId);
239 INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
240 TInt aDim = aCoord.size();
241 for(TInt anId = 0; anId < aDim; anId++)
242 ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
243 const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
244 ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
245 for(TInt anId = 0; anId < aDim; anId++)
246 ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
247 ADDMSG(MYDEBUG,"}\n");
250 BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
257 BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
262 TShapeFun::Eval(const TCellInfo& theCellInfo,
263 const TNodeInfo& theNodeInfo,
264 const TElemNum& theElemNum,
265 const TCCoordSliceArr& theRef,
266 const TCCoordSliceArr& theGauss,
267 TGaussCoord& theGaussCoord,
270 INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);
272 if(IsSatisfy(theRef)){
273 const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
274 TInt aDim = aMeshInfo->GetDim();
275 TInt aNbGauss = theGauss.size();
277 bool anIsSubMesh = !theElemNum.empty();
280 aNbElem = theElemNum.size();
282 aNbElem = theCellInfo.GetNbElem();
284 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
287 InitFun(theRef,theGauss,aFun);
288 TInt aConnDim = theCellInfo.GetConnDim();
290 INITMSG(MYDEBUG,"aDim = "<<aDim<<
291 "; aNbGauss = "<<aNbGauss<<
292 "; aNbElem = "<<aNbElem<<
293 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
296 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
297 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
298 TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
299 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
301 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
302 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
303 TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
305 for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
306 TInt aNodeId = aConnSlice[aConnId] - 1;
307 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
309 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
310 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
318 INITMSG(MYVALUEDEBUG,"theGauss: ");
319 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
320 TCCoordSlice aCoordSlice = theGauss[aGaussId];
321 ADDMSG(MYVALUEDEBUG,"{");
322 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
323 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
325 ADDMSG(MYVALUEDEBUG,"} ");
327 ADDMSG(MYVALUEDEBUG,std::endl);
329 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
330 TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
331 INITMSG(MYVALUEDEBUG,"");
332 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
333 TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
334 ADDMSG(MYVALUEDEBUG,"{");
335 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
336 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
338 ADDMSG(MYVALUEDEBUG,"} ");
340 ADDMSG(MYVALUEDEBUG,std::endl);
350 //---------------------------------------------------------------
351 TSeg2a::TSeg2a():TShapeFun(1,2)
353 TInt aNbRef = GetNbRef();
354 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
355 TCoordSlice aCoord = GetCoord(aRefId);
357 case 0: aCoord[0] = -1.0; break;
358 case 1: aCoord[0] = 1.0; break;
364 TSeg2a::InitFun(const TCCoordSliceArr& theRef,
365 const TCCoordSliceArr& theGauss,
368 GetFun(theRef,theGauss,theFun);
370 TInt aNbGauss = theGauss.size();
371 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
372 const TCCoordSlice& aCoord = theGauss[aGaussId];
373 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
375 aSlice[0] = 0.5*(1.0 - aCoord[0]);
376 aSlice[1] = 0.5*(1.0 + aCoord[0]);
381 //---------------------------------------------------------------
382 TSeg3a::TSeg3a():TShapeFun(1,3)
384 TInt aNbRef = GetNbRef();
385 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
386 TCoordSlice aCoord = GetCoord(aRefId);
388 case 0: aCoord[0] = -1.0; break;
389 case 1: aCoord[0] = 1.0; break;
390 case 2: aCoord[0] = 0.0; break;
396 TSeg3a::InitFun(const TCCoordSliceArr& theRef,
397 const TCCoordSliceArr& theGauss,
400 GetFun(theRef,theGauss,theFun);
402 TInt aNbGauss = theGauss.size();
403 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
404 const TCCoordSlice& aCoord = theGauss[aGaussId];
405 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
407 aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
408 aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
409 aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
415 //---------------------------------------------------------------
419 TInt aNbRef = GetNbRef();
420 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
421 TCoordSlice aCoord = GetCoord(aRefId);
423 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
424 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
425 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
431 TTria3a::InitFun(const TCCoordSliceArr& theRef,
432 const TCCoordSliceArr& theGauss,
435 GetFun(theRef,theGauss,theFun);
437 TInt aNbGauss = theGauss.size();
438 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
439 const TCCoordSlice& aCoord = theGauss[aGaussId];
440 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
442 aSlice[0] = 0.5*(1.0 + aCoord[1]);
443 aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
444 aSlice[2] = 0.5*(1.0 + aCoord[0]);
450 //---------------------------------------------------------------
451 TTria6a::TTria6a():TShapeFun(2,6)
453 TInt aNbRef = GetNbRef();
454 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
455 TCoordSlice aCoord = GetCoord(aRefId);
457 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
458 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
459 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
461 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
462 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
463 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
469 TTria6a::InitFun(const TCCoordSliceArr& theRef,
470 const TCCoordSliceArr& theGauss,
473 GetFun(theRef,theGauss,theFun);
475 TInt aNbGauss = theGauss.size();
476 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
477 const TCCoordSlice& aCoord = theGauss[aGaussId];
478 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
480 aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
481 aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
482 aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
484 aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
485 aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
486 aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
492 //---------------------------------------------------------------
496 TInt aNbRef = GetNbRef();
497 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
498 TCoordSlice aCoord = GetCoord(aRefId);
500 case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
501 case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
502 case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
508 TTria3b::InitFun(const TCCoordSliceArr& theRef,
509 const TCCoordSliceArr& theGauss,
512 GetFun(theRef,theGauss,theFun);
514 TInt aNbGauss = theGauss.size();
515 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
516 const TCCoordSlice& aCoord = theGauss[aGaussId];
517 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
519 aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
520 aSlice[1] = aCoord[0];
521 aSlice[2] = aCoord[1];
527 //---------------------------------------------------------------
531 TInt aNbRef = GetNbRef();
532 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
533 TCoordSlice aCoord = GetCoord(aRefId);
535 case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
536 case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
537 case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
539 case 3: aCoord[0] = 0.5; aCoord[1] = 0.0; break;
540 case 4: aCoord[0] = 0.5; aCoord[1] = 0.5; break;
541 case 5: aCoord[0] = 0.0; aCoord[1] = 0.5; break;
547 TTria6b::InitFun(const TCCoordSliceArr& theRef,
548 const TCCoordSliceArr& theGauss,
551 GetFun(theRef,theGauss,theFun);
553 TInt aNbGauss = theGauss.size();
554 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
555 const TCCoordSlice& aCoord = theGauss[aGaussId];
556 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
558 aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
559 aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
560 aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
562 aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
563 aSlice[4] = 4.0*aCoord[0]*aCoord[1];
564 aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
570 //---------------------------------------------------------------
574 TInt aNbRef = GetNbRef();
575 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
576 TCoordSlice aCoord = GetCoord(aRefId);
578 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
579 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
580 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
581 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
587 TQuad4a::InitFun(const TCCoordSliceArr& theRef,
588 const TCCoordSliceArr& theGauss,
591 GetFun(theRef,theGauss,theFun);
593 TInt aNbGauss = theGauss.size();
594 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
595 const TCCoordSlice& aCoord = theGauss[aGaussId];
596 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
598 aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
599 aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
600 aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
601 aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
607 //---------------------------------------------------------------
611 TInt aNbRef = GetNbRef();
612 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
613 TCoordSlice aCoord = GetCoord(aRefId);
615 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
616 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
617 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
618 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
620 case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
621 case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
622 case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
623 case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
629 TQuad8a::InitFun(const TCCoordSliceArr& theRef,
630 const TCCoordSliceArr& theGauss,
633 GetFun(theRef,theGauss,theFun);
635 TInt aNbGauss = theGauss.size();
636 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
637 const TCCoordSlice& aCoord = theGauss[aGaussId];
638 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
640 aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
641 aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
642 aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
643 aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);
645 aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
646 aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
647 aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
648 aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
654 //---------------------------------------------------------------
658 TInt aNbRef = GetNbRef();
659 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
660 TCoordSlice aCoord = GetCoord(aRefId);
662 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
663 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
664 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
665 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
667 case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
668 case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
669 case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
670 case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
672 case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
678 TQuad9a::InitFun(const TCCoordSliceArr& theRef,
679 const TCCoordSliceArr& theGauss,
682 GetFun(theRef,theGauss,theFun);
684 TInt aNbGauss = theGauss.size();
685 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
686 const TCCoordSlice& aCoord = theGauss[aGaussId];
687 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
689 aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
690 aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
691 aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
692 aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
694 aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
695 aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
696 aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
697 aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
699 aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
705 //---------------------------------------------------------------
709 TInt aNbRef = GetNbRef();
710 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
711 TCoordSlice aCoord = GetCoord(aRefId);
713 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
714 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
715 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
716 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
722 TQuad4b::InitFun(const TCCoordSliceArr& theRef,
723 const TCCoordSliceArr& theGauss,
726 GetFun(theRef,theGauss,theFun);
728 TInt aNbGauss = theGauss.size();
729 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
730 const TCCoordSlice& aCoord = theGauss[aGaussId];
731 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
733 aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
734 aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
735 aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
736 aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
742 //---------------------------------------------------------------
746 TInt aNbRef = GetNbRef();
747 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
748 TCoordSlice aCoord = GetCoord(aRefId);
750 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
751 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
752 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
753 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
755 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
756 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
757 case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
758 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
764 TQuad8b::InitFun(const TCCoordSliceArr& theRef,
765 const TCCoordSliceArr& theGauss,
768 GetFun(theRef,theGauss,theFun);
770 TInt aNbGauss = theGauss.size();
771 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
772 const TCCoordSlice& aCoord = theGauss[aGaussId];
773 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
775 aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
776 aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
777 aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
778 aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);
780 aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
781 aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
782 aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
783 aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);
785 //aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
786 //aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
787 //aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
788 //aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
794 //---------------------------------------------------------------
798 TInt aNbRef = GetNbRef();
799 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
800 TCoordSlice aCoord = GetCoord(aRefId);
802 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
803 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
804 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
805 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
807 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
808 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
809 case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
810 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
812 case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
818 TQuad9b::InitFun(const TCCoordSliceArr& theRef,
819 const TCCoordSliceArr& theGauss,
822 GetFun(theRef,theGauss,theFun);
824 TInt aNbGauss = theGauss.size();
825 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
826 const TCCoordSlice& aCoord = theGauss[aGaussId];
827 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
829 aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
830 aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
831 aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
832 aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
834 aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
835 aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
836 aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
837 aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
839 aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
845 //---------------------------------------------------------------
846 TTetra4a::TTetra4a():
849 TInt aNbRef = GetNbRef();
850 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
851 TCoordSlice aCoord = GetCoord(aRefId);
853 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
854 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
855 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
856 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
862 TTetra4a::InitFun(const TCCoordSliceArr& theRef,
863 const TCCoordSliceArr& theGauss,
866 GetFun(theRef,theGauss,theFun);
868 TInt aNbGauss = theGauss.size();
869 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
870 const TCCoordSlice& aCoord = theGauss[aGaussId];
871 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
873 aSlice[0] = aCoord[1];
874 aSlice[1] = aCoord[2];
875 aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
876 aSlice[3] = aCoord[0];
882 //---------------------------------------------------------------
883 TTetra10a::TTetra10a():
886 TInt aNbRef = GetNbRef();
887 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
888 TCoordSlice aCoord = GetCoord(aRefId);
890 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
891 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
892 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
893 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
895 case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
896 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
897 case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
899 case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
900 case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
901 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
907 TTetra10a::InitFun(const TCCoordSliceArr& theRef,
908 const TCCoordSliceArr& theGauss,
911 GetFun(theRef,theGauss,theFun);
913 TInt aNbGauss = theGauss.size();
914 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
915 const TCCoordSlice& aCoord = theGauss[aGaussId];
916 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
918 aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
919 aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
920 aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
921 aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
923 aSlice[4] = 4.0*aCoord[1]*aCoord[2];
924 aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
925 aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
927 aSlice[7] = 4.0*aCoord[0]*aCoord[1];
928 aSlice[8] = 4.0*aCoord[0]*aCoord[2];
929 aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
935 //---------------------------------------------------------------
938 TTetra4b::TTetra4b():
941 TInt aNbRef = GetNbRef();
942 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
943 TCoordSlice aCoord = GetCoord(aRefId);
945 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
946 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
947 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
948 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
954 TTetra4b::InitFun(const TCCoordSliceArr& theRef,
955 const TCCoordSliceArr& theGauss,
958 GetFun(theRef,theGauss,theFun);
960 TInt aNbGauss = theGauss.size();
961 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
962 const TCCoordSlice& aCoord = theGauss[aGaussId];
963 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
965 aSlice[0] = aCoord[1];
966 aSlice[2] = aCoord[2];
967 aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
968 aSlice[3] = aCoord[0];
974 //---------------------------------------------------------------
975 TTetra10b::TTetra10b():
978 TInt aNbRef = GetNbRef();
979 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
980 TCoordSlice aCoord = GetCoord(aRefId);
982 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
983 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
984 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
985 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
987 case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
988 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
989 case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
991 case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
992 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
993 case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
999 TTetra10b::InitFun(const TCCoordSliceArr& theRef,
1000 const TCCoordSliceArr& theGauss,
1003 GetFun(theRef,theGauss,theFun);
1005 TInt aNbGauss = theGauss.size();
1006 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1007 const TCCoordSlice& aCoord = theGauss[aGaussId];
1008 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1010 aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
1011 aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
1012 aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
1013 aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
1015 aSlice[6] = 4.0*aCoord[1]*aCoord[2];
1016 aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1017 aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1019 aSlice[7] = 4.0*aCoord[0]*aCoord[1];
1020 aSlice[9] = 4.0*aCoord[0]*aCoord[2];
1021 aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1027 //---------------------------------------------------------------
1031 TInt aNbRef = GetNbRef();
1032 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1033 TCoordSlice aCoord = GetCoord(aRefId);
1035 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1036 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1037 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1038 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1039 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1040 case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1041 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1042 case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1048 THexa8a::InitFun(const TCCoordSliceArr& theRef,
1049 const TCCoordSliceArr& theGauss,
1052 GetFun(theRef,theGauss,theFun);
1054 TInt aNbGauss = theGauss.size();
1055 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1056 const TCCoordSlice& aCoord = theGauss[aGaussId];
1057 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1059 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1060 aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1061 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1062 aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1064 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1065 aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1066 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1067 aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1072 //---------------------------------------------------------------
1073 THexa20a::THexa20a(TInt theDim, TInt theNbRef):
1074 TShapeFun(theDim,theNbRef)
1076 TInt aNbRef = myRefCoord.size();
1077 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1078 TCoordSlice aCoord = GetCoord(aRefId);
1080 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1081 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1082 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1083 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1084 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1085 case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1086 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1087 case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1089 case 8: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1090 case 9: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1091 case 10: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1092 case 11: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1093 case 12: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1094 case 13: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1095 case 14: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1096 case 15: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1097 case 16: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1098 case 17: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1099 case 18: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1100 case 19: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1106 THexa20a::InitFun(const TCCoordSliceArr& theRef,
1107 const TCCoordSliceArr& theGauss,
1110 GetFun(theRef,theGauss,theFun);
1112 TInt aNbGauss = theGauss.size();
1113 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1114 const TCCoordSlice& aCoord = theGauss[aGaussId];
1115 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1117 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1118 (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1119 aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1120 (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1121 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1122 (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1123 aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1124 (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1125 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1126 (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1127 aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1128 (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1129 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1130 (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1131 aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1132 (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1134 aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1135 aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1136 aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1137 aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1138 aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1139 aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1140 aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1141 aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1142 aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1143 aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1144 aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1145 aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1151 //---------------------------------------------------------------
1152 THexa27a::THexa27a():
1155 TInt aNbRef = myRefCoord.size();
1156 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1157 TCoordSlice aCoord = GetCoord(aRefId);
1159 case 20: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1160 case 21: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1161 case 22: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1162 case 23: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1163 case 24: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1164 case 25: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1165 case 26: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1171 THexa27a::InitFun(const TCCoordSliceArr& theRef,
1172 const TCCoordSliceArr& theGauss,
1175 GetFun(theRef,theGauss,theFun);
1177 TInt aNbGauss = theGauss.size();
1178 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1179 const TCCoordSlice& aCoord = theGauss[aGaussId];
1180 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1182 aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1183 aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1184 aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1185 aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1186 aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1187 aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1188 aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1189 aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1191 aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1192 aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1193 aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1194 aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1195 aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1196 aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1197 aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1198 aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1199 aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1200 aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1201 aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1202 aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1203 aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1204 aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1205 aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1206 aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1207 aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1208 aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1209 aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
1215 //---------------------------------------------------------------
1219 TInt aNbRef = GetNbRef();
1220 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1221 TCoordSlice aCoord = GetCoord(aRefId);
1223 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1224 case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1225 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1226 case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1227 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1228 case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1229 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1230 case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1236 THexa8b::InitFun(const TCCoordSliceArr& theRef,
1237 const TCCoordSliceArr& theGauss,
1240 GetFun(theRef,theGauss,theFun);
1242 TInt aNbGauss = theGauss.size();
1243 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1244 const TCCoordSlice& aCoord = theGauss[aGaussId];
1245 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1247 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1248 aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1249 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1250 aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1252 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1253 aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1254 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1255 aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1261 //---------------------------------------------------------------
1262 THexa20b::THexa20b(TInt theDim, TInt theNbRef):
1263 TShapeFun(theDim,theNbRef)
1265 TInt aNbRef = myRefCoord.size();
1266 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1267 TCoordSlice aCoord = GetCoord(aRefId);
1269 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1270 case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1271 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1272 case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1273 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1274 case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1275 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1276 case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1278 case 11: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1279 case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1280 case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1281 case 8: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1282 case 16: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1283 case 19: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1284 case 18: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1285 case 17: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1286 case 15: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1287 case 14: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1288 case 13: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1289 case 12: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1295 THexa20b::InitFun(const TCCoordSliceArr& theRef,
1296 const TCCoordSliceArr& theGauss,
1299 GetFun(theRef,theGauss,theFun);
1301 TInt aNbGauss = theGauss.size();
1302 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1303 const TCCoordSlice& aCoord = theGauss[aGaussId];
1304 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1306 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1307 (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1308 aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1309 (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1310 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1311 (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1312 aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1313 (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1314 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1315 (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1316 aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1317 (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1318 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1319 (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1320 aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1321 (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1323 aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1324 aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1325 aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1326 aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1327 aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1328 aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1329 aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1330 aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1331 aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1332 aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1333 aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1334 aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1340 //---------------------------------------------------------------
1341 TPenta6a::TPenta6a():
1344 TInt aNbRef = myRefCoord.size();
1345 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1346 TCoordSlice aCoord = GetCoord(aRefId);
1348 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1349 case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1350 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1351 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1352 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1353 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1359 TPenta6a::InitFun(const TCCoordSliceArr& theRef,
1360 const TCCoordSliceArr& theGauss,
1363 GetFun(theRef,theGauss,theFun);
1365 TInt aNbGauss = theGauss.size();
1366 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1367 const TCCoordSlice& aCoord = theGauss[aGaussId];
1368 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1370 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1371 aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1372 aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1374 aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1375 aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1376 aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1382 //---------------------------------------------------------------
1383 TPenta6b::TPenta6b():
1386 TInt aNbRef = myRefCoord.size();
1387 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1388 TCoordSlice aCoord = GetCoord(aRefId);
1390 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1391 case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1392 case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1393 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1394 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1395 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1401 TPenta6b::InitFun(const TCCoordSliceArr& theRef,
1402 const TCCoordSliceArr& theGauss,
1405 GetFun(theRef,theGauss,theFun);
1407 TInt aNbGauss = theGauss.size();
1408 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1409 const TCCoordSlice& aCoord = theGauss[aGaussId];
1410 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1412 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1413 aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1414 aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1416 aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1417 aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1418 aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1424 //---------------------------------------------------------------
1425 TPenta15a::TPenta15a():
1428 TInt aNbRef = myRefCoord.size();
1429 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1430 TCoordSlice aCoord = GetCoord(aRefId);
1432 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1433 case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1434 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1435 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1436 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1437 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1439 case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1440 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1441 case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1442 case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1443 case 10: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1444 case 11: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1445 case 12: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1446 case 13: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1447 case 14: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1453 TPenta15a::InitFun(const TCCoordSliceArr& theRef,
1454 const TCCoordSliceArr& theGauss,
1457 GetFun(theRef,theGauss,theFun);
1459 TInt aNbGauss = theGauss.size();
1460 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1461 const TCCoordSlice& aCoord = theGauss[aGaussId];
1462 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1464 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1465 aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1466 aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1468 aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1469 aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1470 aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1472 aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1473 aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1474 aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1476 aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1477 aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1478 aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1480 aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1481 aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1482 aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1488 //---------------------------------------------------------------
1489 TPenta15b::TPenta15b():
1492 TInt aNbRef = myRefCoord.size();
1493 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1494 TCoordSlice aCoord = GetCoord(aRefId);
1496 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1497 case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1498 case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1499 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1500 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1501 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1503 case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1504 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1505 case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1506 case 12: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1507 case 14: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1508 case 13: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1509 case 11: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1510 case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1511 case 9: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1517 TPenta15b::InitFun(const TCCoordSliceArr& theRef,
1518 const TCCoordSliceArr& theGauss,
1521 GetFun(theRef,theGauss,theFun);
1523 TInt aNbGauss = theGauss.size();
1524 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1525 const TCCoordSlice& aCoord = theGauss[aGaussId];
1526 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1528 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1529 aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1530 aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1532 aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1533 aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1534 aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1536 aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1537 aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1538 aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1540 aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1541 aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1542 aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1544 aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1545 aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1546 aSlice[9] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1552 //---------------------------------------------------------------
1556 TInt aNbRef = myRefCoord.size();
1557 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1558 TCoordSlice aCoord = GetCoord(aRefId);
1560 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1561 case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1562 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1563 case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1564 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1570 TPyra5a::InitFun(const TCCoordSliceArr& theRef,
1571 const TCCoordSliceArr& theGauss,
1574 GetFun(theRef,theGauss,theFun);
1576 TInt aNbGauss = theGauss.size();
1577 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1578 const TCCoordSlice& aCoord = theGauss[aGaussId];
1579 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1581 // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1582 // Seems shape function for ePYRA5 elements is:
1583 // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1584 // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1585 // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1586 // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1588 aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1589 aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1590 aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1591 aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1592 aSlice[4] = aCoord[2];
1598 //---------------------------------------------------------------
1602 TInt aNbRef = myRefCoord.size();
1603 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1604 TCoordSlice aCoord = GetCoord(aRefId);
1606 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1607 case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1608 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1609 case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1610 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1616 TPyra5b::InitFun(const TCCoordSliceArr& theRef,
1617 const TCCoordSliceArr& theGauss,
1620 GetFun(theRef,theGauss,theFun);
1622 TInt aNbGauss = theGauss.size();
1623 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1624 const TCCoordSlice& aCoord = theGauss[aGaussId];
1625 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1627 // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1628 // Seems shape function for ePYRA5 elements is:
1629 // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1630 // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1631 // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1632 // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1634 aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1635 aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1636 aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1637 aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1638 aSlice[4] = aCoord[2];
1644 //---------------------------------------------------------------
1645 TPyra13a::TPyra13a():
1648 TInt aNbRef = myRefCoord.size();
1649 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1650 TCoordSlice aCoord = GetCoord(aRefId);
1652 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1653 case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1654 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1655 case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1656 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1658 case 5: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1659 case 6: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1660 case 7: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1661 case 8: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1662 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1663 case 10: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1664 case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1665 case 12: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
1671 TPyra13a::InitFun(const TCCoordSliceArr& theRef,
1672 const TCCoordSliceArr& theGauss,
1675 GetFun(theRef,theGauss,theFun);
1677 TInt aNbGauss = theGauss.size();
1678 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1679 const TCCoordSlice& aCoord = theGauss[aGaussId];
1680 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1682 aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1683 (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1684 aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1685 (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1686 aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1687 (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1688 aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1689 (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1691 aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1693 aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1694 (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1695 aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1696 (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1697 aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1698 (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1699 aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1700 (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1702 aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1704 aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1706 aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1708 aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1715 //---------------------------------------------------------------
1716 TPyra13b::TPyra13b():
1719 TInt aNbRef = myRefCoord.size();
1720 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1721 TCoordSlice aCoord = GetCoord(aRefId);
1723 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1724 case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1725 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1726 case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1727 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1729 case 8: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1730 case 7: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1731 case 6: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1732 case 5: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1733 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1734 case 12: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1735 case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1736 case 10: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
1742 TPyra13b::InitFun(const TCCoordSliceArr& theRef,
1743 const TCCoordSliceArr& theGauss,
1746 GetFun(theRef,theGauss,theFun);
1748 TInt aNbGauss = theGauss.size();
1749 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1750 const TCCoordSlice& aCoord = theGauss[aGaussId];
1751 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1753 aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1754 (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1755 aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1756 (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1757 aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1758 (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1759 aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1760 (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1762 aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1764 aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1765 (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1766 aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1767 (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1768 aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1769 (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1770 aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1771 (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1773 aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1775 aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1777 aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1779 aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1786 //---------------------------------------------------------------
1788 GetGaussCoord3D(const TGaussInfo& theGaussInfo,
1789 const TCellInfo& theCellInfo,
1790 const TNodeInfo& theNodeInfo,
1791 TGaussCoord& theGaussCoord,
1792 const TElemNum& theElemNum,
1793 EModeSwitch theMode)
1795 INITMSG(MYDEBUG,"GetGaussCoord3D\n");
1797 if(theGaussInfo.myGeom == theCellInfo.myGeom){
1798 EGeometrieElement aGeom = theGaussInfo.myGeom;
1800 TInt aNbRef = theGaussInfo.GetNbRef();
1801 TCCoordSliceArr aRefSlice(aNbRef);
1802 for(TInt anId = 0; anId < aNbRef; anId++)
1803 aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
1805 TInt aNbGauss = theGaussInfo.GetNbGauss();
1806 TCCoordSliceArr aGaussSlice(aNbGauss);
1807 for(TInt anId = 0; anId < aNbGauss; anId++)
1808 aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
1812 INITMSG(MYDEBUG,"eSEG2"<<std::endl);
1814 if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1820 INITMSG(MYDEBUG,"eSEG3"<<std::endl);
1822 if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1828 INITMSG(MYDEBUG,"eTRIA3"<<std::endl);
1830 if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1833 if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1839 INITMSG(MYDEBUG,"eTRIA6"<<std::endl);
1841 if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1844 if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1850 INITMSG(MYDEBUG,"eQUAD4"<<std::endl);
1852 if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1855 if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1861 INITMSG(MYDEBUG,"eQUAD8"<<std::endl);
1863 if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1866 if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1872 INITMSG(MYDEBUG,"eQUAD9"<<std::endl);
1874 if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1877 if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1883 INITMSG(MYDEBUG,"eTETRA4"<<std::endl);
1885 if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1888 if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1894 INITMSG(MYDEBUG,"ePYRA5"<<std::endl);
1896 if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1899 if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1905 INITMSG(MYDEBUG,"ePENTA6"<<std::endl);
1907 if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1910 if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1916 INITMSG(MYDEBUG,"eHEXA8"<<std::endl);
1918 if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1921 if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1927 INITMSG(MYDEBUG,"eTETRA10"<<std::endl);
1929 if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1932 if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1938 INITMSG(MYDEBUG,"ePYRA13"<<std::endl);
1940 if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1943 if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1949 INITMSG(MYDEBUG,"ePENTA15"<<std::endl);
1951 if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1954 if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1960 INITMSG(MYDEBUG,"eHEXA20"<<std::endl);
1962 if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1965 if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1971 INITMSG(MYDEBUG,"eNONE"<<std::endl);
1979 //---------------------------------------------------------------
1981 GetBaryCenter(const TCellInfo& theCellInfo,
1982 const TNodeInfo& theNodeInfo,
1983 TGaussCoord& theGaussCoord,
1984 const TElemNum& theElemNum,
1985 EModeSwitch theMode)
1987 INITMSG(MYDEBUG,"GetBaryCenter\n");
1988 const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
1989 TInt aDim = aMeshInfo->GetDim();
1990 static TInt aNbGauss = 1;
1992 bool anIsSubMesh = !theElemNum.empty();
1995 aNbElem = theElemNum.size();
1997 aNbElem = theCellInfo.GetNbElem();
1999 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2001 TInt aConnDim = theCellInfo.GetConnDim();
2005 "; aNbGauss = "<<aNbGauss<<
2006 "; aNbElem = "<<aNbElem<<
2007 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2010 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2011 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2012 TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
2013 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2015 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2016 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2018 for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
2019 TInt aNodeId = aConnSlice[aConnId] - 1;
2020 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2022 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2023 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2027 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2028 aGaussCoordSlice[aDimId] /= aConnDim;
2034 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2035 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2036 INITMSG(MYVALUEDEBUG,"");
2037 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2038 TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
2039 ADDMSG(MYVALUEDEBUG,"{");
2040 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2041 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
2043 ADDMSG(MYVALUEDEBUG,"} ");
2045 ADDMSG(MYVALUEDEBUG,std::endl);
2053 //---------------------------------------------------------------
2055 GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
2056 const TNodeInfo& theNodeInfo,
2057 TGaussCoord& theGaussCoord,
2058 const TElemNum& theElemNum,
2059 EModeSwitch theMode)
2061 INITMSG(MYDEBUG,"GetBaryCenter\n");
2062 const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
2063 TInt aDim = aMeshInfo->GetDim();
2064 static TInt aNbGauss = 1;
2066 bool anIsSubMesh = !theElemNum.empty();
2069 aNbElem = theElemNum.size();
2071 aNbElem = thePolygoneInfo.GetNbElem();
2073 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2077 "; aNbGauss = "<<aNbGauss<<
2078 "; aNbElem = "<<aNbElem<<
2079 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2082 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2083 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2085 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2086 TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
2087 TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
2088 TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
2090 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2091 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2093 for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2094 TInt aNodeId = aConnSlice[aConnId] - 1;
2095 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2097 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2098 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2102 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2103 aGaussCoordSlice[aDimId] /= aNbNodes;
2112 //---------------------------------------------------------------
2114 GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
2115 const TNodeInfo& theNodeInfo,
2116 TGaussCoord& theGaussCoord,
2117 const TElemNum& theElemNum,
2118 EModeSwitch theMode)
2120 INITMSG(MYDEBUG,"GetBaryCenter\n");
2121 const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
2122 TInt aDim = aMeshInfo->GetDim();
2123 static TInt aNbGauss = 1;
2125 bool anIsSubMesh = !theElemNum.empty();
2128 aNbElem = theElemNum.size();
2130 aNbElem = thePolyedreInfo.GetNbElem();
2132 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2136 "; aNbGauss = "<<aNbGauss<<
2137 "; aNbElem = "<<aNbElem<<
2138 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2141 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2142 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2144 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2145 TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
2146 TInt aNbFaces = aConnSliceArr.size();
2148 TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
2150 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2151 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2153 for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
2154 TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
2155 TInt aNbConn = aConnSlice.size();
2156 for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2157 TInt aNodeId = aConnSlice[aConnId] - 1;
2158 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2160 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2161 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2165 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2166 aGaussCoordSlice[aDimId] /= aNbNodes;