1 // Copyright (C) 2007-2020 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, or (at your option) any later version.
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;
36 //---------------------------------------------------------------
39 TModeSwitchInfo(eFULL_INTERLACE),
49 ::Init(TInt theNbElem,
54 myModeSwitch = theMode;
57 myNbGauss = theNbGauss;
60 myGaussStep = myNbGauss*myDim;
62 myGaussCoord.resize(theNbElem*myGaussStep);
90 return (unsigned char*)&(myGaussCoord[0]);
95 ::GetCoordSliceArr(TInt theElemId) const
97 TCCoordSliceArr aCoordSliceArr(myNbGauss);
98 if(GetModeSwitch() == eFULL_INTERLACE){
99 TInt anId = theElemId*myGaussStep;
100 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
101 aCoordSliceArr[anGaussId] =
102 TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
107 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
108 aCoordSliceArr[anGaussId] =
109 TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
112 return aCoordSliceArr;
117 ::GetCoordSliceArr(TInt theElemId)
119 TCoordSliceArr aCoordSliceArr(myNbGauss);
120 if(GetModeSwitch() == eFULL_INTERLACE){
121 TInt anId = theElemId*myGaussStep;
122 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
123 aCoordSliceArr[anGaussId] =
124 TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
129 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
130 aCoordSliceArr[anGaussId] =
131 TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
134 return aCoordSliceArr;
137 //---------------------------------------------------------------
140 IsEqual(TFloat theLeft, TFloat theRight)
142 static TFloat EPS = 1.0E-3;
143 if(fabs(theLeft) + fabs(theRight) > EPS)
144 return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
148 //---------------------------------------------------------------
149 class TShapeFun::TFun
157 Init(TInt theNbGauss,
160 myFun.resize(theNbGauss*theNbRef);
165 GetFunSlice(TInt theGaussId) const
167 return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
171 GetFunSlice(TInt theGaussId)
173 return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
177 //---------------------------------------------------------------
179 TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
180 myRefCoord(theNbRef*theDim),
186 TShapeFun::GetCoord(TInt theRefId) const
188 return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
192 TShapeFun::GetCoord(TInt theRefId)
194 return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
198 TShapeFun::GetFun(const TCCoordSliceArr& theRef,
199 const TCCoordSliceArr& theGauss,
202 TInt aNbRef = theRef.size();
203 TInt aNbGauss = theGauss.size();
204 theFun.Init(aNbGauss,aNbRef);
208 TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
210 TInt aNbRef = theRefCoord.size();
211 TInt aNbRef2 = GetNbRef();
212 INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
213 "- aNbRef("<<aNbRef<<")"<<
214 "; aNbRef2("<<aNbRef2<<")\n");
215 bool anIsSatisfy = (aNbRef == aNbRef2);
217 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
218 const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
219 TCCoordSlice aCoord = GetCoord(aRefId);
220 TInt aDim = aCoord.size();
221 bool anIsEqual = false;
222 for(TInt anId = 0; anId < aDim; anId++){
223 anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
231 TCCoordSlice aCoord = GetCoord(aRefId);
232 INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
233 TInt aDim = aCoord.size();
234 for(TInt anId = 0; anId < aDim; anId++)
235 ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
236 const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
237 ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
238 for(TInt anId = 0; anId < aDim; anId++)
239 ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
240 ADDMSG(MYDEBUG,"}\n");
243 BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
250 BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
255 TShapeFun::Eval(const TCellInfo& theCellInfo,
256 const TNodeInfo& theNodeInfo,
257 const TElemNum& theElemNum,
258 const TCCoordSliceArr& theRef,
259 const TCCoordSliceArr& theGauss,
260 TGaussCoord& theGaussCoord,
263 INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);
265 if(IsSatisfy(theRef)){
266 const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
267 TInt aDim = aMeshInfo->GetDim();
268 TInt aNbGauss = theGauss.size();
270 bool anIsSubMesh = !theElemNum.empty();
273 aNbElem = theElemNum.size();
275 aNbElem = theCellInfo.GetNbElem();
277 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
280 InitFun(theRef,theGauss,aFun);
281 TInt aConnDim = theCellInfo.GetConnDim();
283 INITMSG(MYDEBUG,"aDim = "<<aDim<<
284 "; aNbGauss = "<<aNbGauss<<
285 "; aNbElem = "<<aNbElem<<
286 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
289 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
290 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
291 TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
292 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
294 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
295 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
296 TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
298 for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
299 TInt aNodeId = aConnSlice[aConnId] - 1;
300 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
302 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
303 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
311 INITMSG(MYVALUEDEBUG,"theGauss: ");
312 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
313 TCCoordSlice aCoordSlice = theGauss[aGaussId];
314 ADDMSG(MYVALUEDEBUG,"{");
315 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
316 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
318 ADDMSG(MYVALUEDEBUG,"} ");
320 ADDMSG(MYVALUEDEBUG,std::endl);
322 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
323 TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
324 INITMSG(MYVALUEDEBUG,"");
325 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
326 TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
327 ADDMSG(MYVALUEDEBUG,"{");
328 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
329 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
331 ADDMSG(MYVALUEDEBUG,"} ");
333 ADDMSG(MYVALUEDEBUG,std::endl);
342 //---------------------------------------------------------------
343 TSeg2a::TSeg2a():TShapeFun(1,2)
345 TInt aNbRef = GetNbRef();
346 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
347 TCoordSlice aCoord = GetCoord(aRefId);
349 case 0: aCoord[0] = -1.0; break;
350 case 1: aCoord[0] = 1.0; break;
356 TSeg2a::InitFun(const TCCoordSliceArr& theRef,
357 const TCCoordSliceArr& theGauss,
360 GetFun(theRef,theGauss,theFun);
362 TInt aNbGauss = theGauss.size();
363 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
364 const TCCoordSlice& aCoord = theGauss[aGaussId];
365 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
367 aSlice[0] = 0.5*(1.0 - aCoord[0]);
368 aSlice[1] = 0.5*(1.0 + aCoord[0]);
372 //---------------------------------------------------------------
373 TSeg3a::TSeg3a():TShapeFun(1,3)
375 TInt aNbRef = GetNbRef();
376 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
377 TCoordSlice aCoord = GetCoord(aRefId);
379 case 0: aCoord[0] = -1.0; break;
380 case 1: aCoord[0] = 1.0; break;
381 case 2: aCoord[0] = 0.0; break;
387 TSeg3a::InitFun(const TCCoordSliceArr& theRef,
388 const TCCoordSliceArr& theGauss,
391 GetFun(theRef,theGauss,theFun);
393 TInt aNbGauss = theGauss.size();
394 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
395 const TCCoordSlice& aCoord = theGauss[aGaussId];
396 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
398 aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
399 aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
400 aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
404 //---------------------------------------------------------------
408 TInt aNbRef = GetNbRef();
409 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
410 TCoordSlice aCoord = GetCoord(aRefId);
412 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
413 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
414 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
420 TTria3a::InitFun(const TCCoordSliceArr& theRef,
421 const TCCoordSliceArr& theGauss,
424 GetFun(theRef,theGauss,theFun);
426 TInt aNbGauss = theGauss.size();
427 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
428 const TCCoordSlice& aCoord = theGauss[aGaussId];
429 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
431 aSlice[0] = 0.5*(1.0 + aCoord[1]);
432 aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
433 aSlice[2] = 0.5*(1.0 + aCoord[0]);
437 //---------------------------------------------------------------
438 TTria6a::TTria6a():TShapeFun(2,6)
440 TInt aNbRef = GetNbRef();
441 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
442 TCoordSlice aCoord = GetCoord(aRefId);
444 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
445 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
446 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
448 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
449 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
450 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
456 TTria6a::InitFun(const TCCoordSliceArr& theRef,
457 const TCCoordSliceArr& theGauss,
460 GetFun(theRef,theGauss,theFun);
462 TInt aNbGauss = theGauss.size();
463 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
464 const TCCoordSlice& aCoord = theGauss[aGaussId];
465 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
467 aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
468 aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
469 aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
471 aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
472 aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
473 aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
477 //---------------------------------------------------------------
481 TInt aNbRef = GetNbRef();
482 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
483 TCoordSlice aCoord = GetCoord(aRefId);
485 case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
486 case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
487 case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
493 TTria3b::InitFun(const TCCoordSliceArr& theRef,
494 const TCCoordSliceArr& theGauss,
497 GetFun(theRef,theGauss,theFun);
499 TInt aNbGauss = theGauss.size();
500 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
501 const TCCoordSlice& aCoord = theGauss[aGaussId];
502 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
504 aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
505 aSlice[1] = aCoord[0];
506 aSlice[2] = aCoord[1];
510 //---------------------------------------------------------------
514 TInt aNbRef = GetNbRef();
515 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
516 TCoordSlice aCoord = GetCoord(aRefId);
518 case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
519 case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
520 case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
522 case 3: aCoord[0] = 0.5; aCoord[1] = 0.0; break;
523 case 4: aCoord[0] = 0.5; aCoord[1] = 0.5; break;
524 case 5: aCoord[0] = 0.0; aCoord[1] = 0.5; break;
530 TTria6b::InitFun(const TCCoordSliceArr& theRef,
531 const TCCoordSliceArr& theGauss,
534 GetFun(theRef,theGauss,theFun);
536 TInt aNbGauss = theGauss.size();
537 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
538 const TCCoordSlice& aCoord = theGauss[aGaussId];
539 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
541 aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
542 aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
543 aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
545 aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
546 aSlice[4] = 4.0*aCoord[0]*aCoord[1];
547 aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
551 //---------------------------------------------------------------
555 TInt aNbRef = GetNbRef();
556 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
557 TCoordSlice aCoord = GetCoord(aRefId);
559 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
560 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
561 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
562 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
568 TQuad4a::InitFun(const TCCoordSliceArr& theRef,
569 const TCCoordSliceArr& theGauss,
572 GetFun(theRef,theGauss,theFun);
574 TInt aNbGauss = theGauss.size();
575 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
576 const TCCoordSlice& aCoord = theGauss[aGaussId];
577 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
579 aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
580 aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
581 aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
582 aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
586 //---------------------------------------------------------------
590 TInt aNbRef = GetNbRef();
591 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
592 TCoordSlice aCoord = GetCoord(aRefId);
594 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
595 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
596 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
597 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
599 case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
600 case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
601 case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
602 case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
608 TQuad8a::InitFun(const TCCoordSliceArr& theRef,
609 const TCCoordSliceArr& theGauss,
612 GetFun(theRef,theGauss,theFun);
614 TInt aNbGauss = theGauss.size();
615 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
616 const TCCoordSlice& aCoord = theGauss[aGaussId];
617 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
619 aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
620 aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
621 aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
622 aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);
624 aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
625 aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
626 aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
627 aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
631 //---------------------------------------------------------------
635 TInt aNbRef = GetNbRef();
636 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
637 TCoordSlice aCoord = GetCoord(aRefId);
639 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
640 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
641 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
642 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
644 case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
645 case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
646 case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
647 case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
649 case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
655 TQuad9a::InitFun(const TCCoordSliceArr& theRef,
656 const TCCoordSliceArr& theGauss,
659 GetFun(theRef,theGauss,theFun);
661 TInt aNbGauss = theGauss.size();
662 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
663 const TCCoordSlice& aCoord = theGauss[aGaussId];
664 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
666 aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
667 aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
668 aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
669 aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
671 aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
672 aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
673 aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
674 aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
676 aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
680 //---------------------------------------------------------------
684 TInt aNbRef = GetNbRef();
685 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
686 TCoordSlice aCoord = GetCoord(aRefId);
688 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
689 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
690 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
691 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
697 TQuad4b::InitFun(const TCCoordSliceArr& theRef,
698 const TCCoordSliceArr& theGauss,
701 GetFun(theRef,theGauss,theFun);
703 TInt aNbGauss = theGauss.size();
704 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
705 const TCCoordSlice& aCoord = theGauss[aGaussId];
706 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
708 aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
709 aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
710 aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
711 aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
715 //---------------------------------------------------------------
719 TInt aNbRef = GetNbRef();
720 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
721 TCoordSlice aCoord = GetCoord(aRefId);
723 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
724 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
725 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
726 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
728 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
729 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
730 case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
731 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
737 TQuad8b::InitFun(const TCCoordSliceArr& theRef,
738 const TCCoordSliceArr& theGauss,
741 GetFun(theRef,theGauss,theFun);
743 TInt aNbGauss = theGauss.size();
744 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
745 const TCCoordSlice& aCoord = theGauss[aGaussId];
746 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
748 aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
749 aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
750 aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
751 aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);
753 aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
754 aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
755 aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
756 aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);
758 //aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
759 //aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
760 //aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
761 //aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
765 //---------------------------------------------------------------
769 TInt aNbRef = GetNbRef();
770 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
771 TCoordSlice aCoord = GetCoord(aRefId);
773 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
774 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
775 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
776 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
778 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
779 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
780 case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
781 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
783 case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
789 TQuad9b::InitFun(const TCCoordSliceArr& theRef,
790 const TCCoordSliceArr& theGauss,
793 GetFun(theRef,theGauss,theFun);
795 TInt aNbGauss = theGauss.size();
796 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
797 const TCCoordSlice& aCoord = theGauss[aGaussId];
798 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
800 aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
801 aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
802 aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
803 aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
805 aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
806 aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
807 aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
808 aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
810 aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
814 //---------------------------------------------------------------
815 TTetra4a::TTetra4a():
818 TInt aNbRef = GetNbRef();
819 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
820 TCoordSlice aCoord = GetCoord(aRefId);
822 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
823 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
824 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
825 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
831 TTetra4a::InitFun(const TCCoordSliceArr& theRef,
832 const TCCoordSliceArr& theGauss,
835 GetFun(theRef,theGauss,theFun);
837 TInt aNbGauss = theGauss.size();
838 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
839 const TCCoordSlice& aCoord = theGauss[aGaussId];
840 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
842 aSlice[0] = aCoord[1];
843 aSlice[1] = aCoord[2];
844 aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
845 aSlice[3] = aCoord[0];
849 //---------------------------------------------------------------
850 TTetra10a::TTetra10a():
853 TInt aNbRef = GetNbRef();
854 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
855 TCoordSlice aCoord = GetCoord(aRefId);
857 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
858 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
859 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
860 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
862 case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
863 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
864 case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
866 case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
867 case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
868 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
874 TTetra10a::InitFun(const TCCoordSliceArr& theRef,
875 const TCCoordSliceArr& theGauss,
878 GetFun(theRef,theGauss,theFun);
880 TInt aNbGauss = theGauss.size();
881 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
882 const TCCoordSlice& aCoord = theGauss[aGaussId];
883 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
885 aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
886 aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
887 aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
888 aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
890 aSlice[4] = 4.0*aCoord[1]*aCoord[2];
891 aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
892 aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
894 aSlice[7] = 4.0*aCoord[0]*aCoord[1];
895 aSlice[8] = 4.0*aCoord[0]*aCoord[2];
896 aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
900 //---------------------------------------------------------------
902 TTetra4b::TTetra4b():
905 TInt aNbRef = GetNbRef();
906 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
907 TCoordSlice aCoord = GetCoord(aRefId);
909 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
910 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
911 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
912 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
918 TTetra4b::InitFun(const TCCoordSliceArr& theRef,
919 const TCCoordSliceArr& theGauss,
922 GetFun(theRef,theGauss,theFun);
924 TInt aNbGauss = theGauss.size();
925 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
926 const TCCoordSlice& aCoord = theGauss[aGaussId];
927 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
929 aSlice[0] = aCoord[1];
930 aSlice[2] = aCoord[2];
931 aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
932 aSlice[3] = aCoord[0];
936 //---------------------------------------------------------------
937 TTetra10b::TTetra10b():
940 TInt aNbRef = GetNbRef();
941 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
942 TCoordSlice aCoord = GetCoord(aRefId);
944 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
945 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
946 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
947 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
949 case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
950 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
951 case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
953 case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
954 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
955 case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
961 TTetra10b::InitFun(const TCCoordSliceArr& theRef,
962 const TCCoordSliceArr& theGauss,
965 GetFun(theRef,theGauss,theFun);
967 TInt aNbGauss = theGauss.size();
968 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
969 const TCCoordSlice& aCoord = theGauss[aGaussId];
970 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
972 aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
973 aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
974 aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
975 aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
977 aSlice[6] = 4.0*aCoord[1]*aCoord[2];
978 aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
979 aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
981 aSlice[7] = 4.0*aCoord[0]*aCoord[1];
982 aSlice[9] = 4.0*aCoord[0]*aCoord[2];
983 aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
987 //---------------------------------------------------------------
991 TInt aNbRef = GetNbRef();
992 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
993 TCoordSlice aCoord = GetCoord(aRefId);
995 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
996 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
997 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
998 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
999 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1000 case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1001 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1002 case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1008 THexa8a::InitFun(const TCCoordSliceArr& theRef,
1009 const TCCoordSliceArr& theGauss,
1012 GetFun(theRef,theGauss,theFun);
1014 TInt aNbGauss = theGauss.size();
1015 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1016 const TCCoordSlice& aCoord = theGauss[aGaussId];
1017 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1019 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1020 aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1021 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1022 aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1024 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1025 aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1026 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1027 aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1031 //---------------------------------------------------------------
1032 THexa20a::THexa20a(TInt theDim, TInt theNbRef):
1033 TShapeFun(theDim,theNbRef)
1035 TInt aNbRef = myRefCoord.size();
1036 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1037 TCoordSlice aCoord = GetCoord(aRefId);
1039 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1040 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1041 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1042 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1043 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1044 case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1045 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1046 case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1048 case 8: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1049 case 9: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1050 case 10: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1051 case 11: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1052 case 12: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1053 case 13: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1054 case 14: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1055 case 15: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1056 case 16: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1057 case 17: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1058 case 18: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1059 case 19: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1065 THexa20a::InitFun(const TCCoordSliceArr& theRef,
1066 const TCCoordSliceArr& theGauss,
1069 GetFun(theRef,theGauss,theFun);
1071 TInt aNbGauss = theGauss.size();
1072 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1073 const TCCoordSlice& aCoord = theGauss[aGaussId];
1074 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1076 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1077 (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1078 aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1079 (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1080 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1081 (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1082 aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1083 (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1084 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1085 (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1086 aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1087 (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1088 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1089 (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1090 aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1091 (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1093 aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1094 aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1095 aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1096 aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1097 aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1098 aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1099 aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1100 aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1101 aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1102 aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1103 aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1104 aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1108 //---------------------------------------------------------------
1109 THexa27a::THexa27a():
1112 TInt aNbRef = myRefCoord.size();
1113 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1114 TCoordSlice aCoord = GetCoord(aRefId);
1116 case 20: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1117 case 21: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1118 case 22: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1119 case 23: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1120 case 24: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1121 case 25: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1122 case 26: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1128 THexa27a::InitFun(const TCCoordSliceArr& theRef,
1129 const TCCoordSliceArr& theGauss,
1132 GetFun(theRef,theGauss,theFun);
1134 TInt aNbGauss = theGauss.size();
1135 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1136 const TCCoordSlice& aCoord = theGauss[aGaussId];
1137 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1139 aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1140 aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1141 aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1142 aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1143 aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1144 aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1145 aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1146 aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1148 aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1149 aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1150 aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1151 aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1152 aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1153 aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1154 aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1155 aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1156 aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1157 aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1158 aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1159 aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1160 aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1161 aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1162 aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1163 aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1164 aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1165 aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1166 aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
1170 //---------------------------------------------------------------
1174 TInt aNbRef = GetNbRef();
1175 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1176 TCoordSlice aCoord = GetCoord(aRefId);
1178 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1179 case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1180 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1181 case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1182 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1183 case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1184 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1185 case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1191 THexa8b::InitFun(const TCCoordSliceArr& theRef,
1192 const TCCoordSliceArr& theGauss,
1195 GetFun(theRef,theGauss,theFun);
1197 TInt aNbGauss = theGauss.size();
1198 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1199 const TCCoordSlice& aCoord = theGauss[aGaussId];
1200 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1202 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1203 aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1204 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1205 aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1207 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1208 aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1209 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1210 aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1214 //---------------------------------------------------------------
1215 THexa20b::THexa20b(TInt theDim, TInt theNbRef):
1216 TShapeFun(theDim,theNbRef)
1218 TInt aNbRef = myRefCoord.size();
1219 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1220 TCoordSlice aCoord = GetCoord(aRefId);
1222 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1223 case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1224 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1225 case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1226 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1227 case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1228 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1229 case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1231 case 11: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1232 case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1233 case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1234 case 8: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1235 case 16: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1236 case 19: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1237 case 18: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1238 case 17: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1239 case 15: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1240 case 14: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1241 case 13: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1242 case 12: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1248 THexa20b::InitFun(const TCCoordSliceArr& theRef,
1249 const TCCoordSliceArr& theGauss,
1252 GetFun(theRef,theGauss,theFun);
1254 TInt aNbGauss = theGauss.size();
1255 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1256 const TCCoordSlice& aCoord = theGauss[aGaussId];
1257 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1259 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1260 (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1261 aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1262 (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1263 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1264 (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1265 aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1266 (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1267 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1268 (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1269 aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1270 (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1271 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1272 (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1273 aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1274 (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1276 aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1277 aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1278 aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1279 aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1280 aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1281 aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1282 aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1283 aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1284 aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1285 aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1286 aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1287 aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1291 //---------------------------------------------------------------
1292 TPenta6a::TPenta6a():
1295 TInt aNbRef = myRefCoord.size();
1296 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1297 TCoordSlice aCoord = GetCoord(aRefId);
1299 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1300 case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1301 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1302 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1303 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1304 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1310 TPenta6a::InitFun(const TCCoordSliceArr& theRef,
1311 const TCCoordSliceArr& theGauss,
1314 GetFun(theRef,theGauss,theFun);
1316 TInt aNbGauss = theGauss.size();
1317 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1318 const TCCoordSlice& aCoord = theGauss[aGaussId];
1319 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1321 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1322 aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1323 aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1325 aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1326 aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1327 aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1331 //---------------------------------------------------------------
1332 TPenta6b::TPenta6b():
1335 TInt aNbRef = myRefCoord.size();
1336 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1337 TCoordSlice aCoord = GetCoord(aRefId);
1339 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1340 case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1341 case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1342 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1343 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1344 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1350 TPenta6b::InitFun(const TCCoordSliceArr& theRef,
1351 const TCCoordSliceArr& theGauss,
1354 GetFun(theRef,theGauss,theFun);
1356 TInt aNbGauss = theGauss.size();
1357 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1358 const TCCoordSlice& aCoord = theGauss[aGaussId];
1359 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1361 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1362 aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1363 aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1365 aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1366 aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1367 aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1371 //---------------------------------------------------------------
1372 TPenta15a::TPenta15a():
1375 TInt aNbRef = myRefCoord.size();
1376 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1377 TCoordSlice aCoord = GetCoord(aRefId);
1379 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1380 case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1381 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1382 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1383 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1384 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1386 case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1387 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1388 case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1389 case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1390 case 10: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1391 case 11: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1392 case 12: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1393 case 13: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1394 case 14: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1400 TPenta15a::InitFun(const TCCoordSliceArr& theRef,
1401 const TCCoordSliceArr& theGauss,
1404 GetFun(theRef,theGauss,theFun);
1406 TInt aNbGauss = theGauss.size();
1407 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1408 const TCCoordSlice& aCoord = theGauss[aGaussId];
1409 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1411 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1412 aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1413 aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1415 aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1416 aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1417 aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1419 aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1420 aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1421 aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1423 aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1424 aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1425 aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1427 aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1428 aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1429 aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1433 //---------------------------------------------------------------
1434 TPenta15b::TPenta15b():
1437 TInt aNbRef = myRefCoord.size();
1438 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1439 TCoordSlice aCoord = GetCoord(aRefId);
1441 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1442 case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1443 case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1444 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1445 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1446 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1448 case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1449 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1450 case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1451 case 12: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1452 case 14: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1453 case 13: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1454 case 11: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1455 case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1456 case 9: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1462 TPenta15b::InitFun(const TCCoordSliceArr& theRef,
1463 const TCCoordSliceArr& theGauss,
1466 GetFun(theRef,theGauss,theFun);
1468 TInt aNbGauss = theGauss.size();
1469 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1470 const TCCoordSlice& aCoord = theGauss[aGaussId];
1471 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1473 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1474 aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1475 aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1477 aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1478 aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1479 aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1481 aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1482 aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1483 aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1485 aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1486 aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1487 aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1489 aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1490 aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1491 aSlice[9] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1495 //---------------------------------------------------------------
1499 TInt aNbRef = myRefCoord.size();
1500 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1501 TCoordSlice aCoord = GetCoord(aRefId);
1503 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1504 case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1505 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1506 case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1507 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1513 TPyra5a::InitFun(const TCCoordSliceArr& theRef,
1514 const TCCoordSliceArr& theGauss,
1517 GetFun(theRef,theGauss,theFun);
1519 TInt aNbGauss = theGauss.size();
1520 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1521 const TCCoordSlice& aCoord = theGauss[aGaussId];
1522 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1524 // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1525 // Seems shape function for ePYRA5 elements is:
1526 // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1527 // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1528 // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1529 // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1531 aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1532 aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1533 aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1534 aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1535 aSlice[4] = aCoord[2];
1539 //---------------------------------------------------------------
1543 TInt aNbRef = myRefCoord.size();
1544 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1545 TCoordSlice aCoord = GetCoord(aRefId);
1547 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1548 case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1549 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1550 case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1551 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1557 TPyra5b::InitFun(const TCCoordSliceArr& theRef,
1558 const TCCoordSliceArr& theGauss,
1561 GetFun(theRef,theGauss,theFun);
1563 TInt aNbGauss = theGauss.size();
1564 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1565 const TCCoordSlice& aCoord = theGauss[aGaussId];
1566 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1568 // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1569 // Seems shape function for ePYRA5 elements is:
1570 // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1571 // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1572 // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1573 // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1575 aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1576 aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1577 aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1578 aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1579 aSlice[4] = aCoord[2];
1583 //---------------------------------------------------------------
1584 TPyra13a::TPyra13a():
1587 TInt aNbRef = myRefCoord.size();
1588 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1589 TCoordSlice aCoord = GetCoord(aRefId);
1591 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1592 case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1593 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1594 case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1595 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1597 case 5: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1598 case 6: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1599 case 7: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1600 case 8: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1601 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1602 case 10: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1603 case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1604 case 12: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
1610 TPyra13a::InitFun(const TCCoordSliceArr& theRef,
1611 const TCCoordSliceArr& theGauss,
1614 GetFun(theRef,theGauss,theFun);
1616 TInt aNbGauss = theGauss.size();
1617 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1618 const TCCoordSlice& aCoord = theGauss[aGaussId];
1619 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1621 aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1622 (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1623 aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1624 (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1625 aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1626 (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1627 aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1628 (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1630 aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1632 aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1633 (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1634 aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1635 (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1636 aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1637 (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1638 aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1639 (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1641 aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1643 aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1645 aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1647 aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1652 //---------------------------------------------------------------
1653 TPyra13b::TPyra13b():
1656 TInt aNbRef = myRefCoord.size();
1657 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1658 TCoordSlice aCoord = GetCoord(aRefId);
1660 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1661 case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1662 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1663 case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1664 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1666 case 8: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1667 case 7: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1668 case 6: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1669 case 5: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1670 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1671 case 12: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1672 case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1673 case 10: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
1679 TPyra13b::InitFun(const TCCoordSliceArr& theRef,
1680 const TCCoordSliceArr& theGauss,
1683 GetFun(theRef,theGauss,theFun);
1685 TInt aNbGauss = theGauss.size();
1686 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1687 const TCCoordSlice& aCoord = theGauss[aGaussId];
1688 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1690 aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1691 (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1692 aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1693 (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1694 aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1695 (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1696 aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1697 (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1699 aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1701 aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1702 (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1703 aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1704 (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1705 aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1706 (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1707 aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1708 (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1710 aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1712 aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1714 aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1716 aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1721 //---------------------------------------------------------------
1723 GetGaussCoord3D(const TGaussInfo& theGaussInfo,
1724 const TCellInfo& theCellInfo,
1725 const TNodeInfo& theNodeInfo,
1726 TGaussCoord& theGaussCoord,
1727 const TElemNum& theElemNum,
1728 EModeSwitch theMode)
1730 INITMSG(MYDEBUG,"GetGaussCoord3D\n");
1732 if(theGaussInfo.myGeom == theCellInfo.myGeom){
1733 EGeometrieElement aGeom = theGaussInfo.myGeom;
1735 TInt aNbRef = theGaussInfo.GetNbRef();
1736 TCCoordSliceArr aRefSlice(aNbRef);
1737 for(TInt anId = 0; anId < aNbRef; anId++)
1738 aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
1740 TInt aNbGauss = theGaussInfo.GetNbGauss();
1741 TCCoordSliceArr aGaussSlice(aNbGauss);
1742 for(TInt anId = 0; anId < aNbGauss; anId++)
1743 aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
1747 INITMSG(MYDEBUG,"eSEG2"<<std::endl);
1749 if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1755 INITMSG(MYDEBUG,"eSEG3"<<std::endl);
1757 if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1763 INITMSG(MYDEBUG,"eTRIA3"<<std::endl);
1765 if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1768 if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1774 INITMSG(MYDEBUG,"eTRIA6"<<std::endl);
1776 if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1779 if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1785 INITMSG(MYDEBUG,"eQUAD4"<<std::endl);
1787 if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1790 if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1796 INITMSG(MYDEBUG,"eQUAD8"<<std::endl);
1798 if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1801 if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1807 INITMSG(MYDEBUG,"eQUAD9"<<std::endl);
1809 if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1812 if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1818 INITMSG(MYDEBUG,"eTETRA4"<<std::endl);
1820 if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1823 if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1829 INITMSG(MYDEBUG,"ePYRA5"<<std::endl);
1831 if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1834 if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1840 INITMSG(MYDEBUG,"ePENTA6"<<std::endl);
1842 if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1845 if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1851 INITMSG(MYDEBUG,"eHEXA8"<<std::endl);
1853 if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1856 if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1862 INITMSG(MYDEBUG,"eTETRA10"<<std::endl);
1864 if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1867 if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1873 INITMSG(MYDEBUG,"ePYRA13"<<std::endl);
1875 if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1878 if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1884 INITMSG(MYDEBUG,"ePENTA15"<<std::endl);
1886 if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1889 if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1895 INITMSG(MYDEBUG,"eHEXA20"<<std::endl);
1897 if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1900 if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1906 INITMSG(MYDEBUG,"eNONE"<<std::endl);
1914 //---------------------------------------------------------------
1916 GetBaryCenter(const TCellInfo& theCellInfo,
1917 const TNodeInfo& theNodeInfo,
1918 TGaussCoord& theGaussCoord,
1919 const TElemNum& theElemNum,
1920 EModeSwitch theMode)
1922 INITMSG(MYDEBUG,"GetBaryCenter\n");
1923 const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
1924 TInt aDim = aMeshInfo->GetDim();
1925 static TInt aNbGauss = 1;
1927 bool anIsSubMesh = !theElemNum.empty();
1930 aNbElem = theElemNum.size();
1932 aNbElem = theCellInfo.GetNbElem();
1934 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
1936 TInt aConnDim = theCellInfo.GetConnDim();
1940 "; aNbGauss = "<<aNbGauss<<
1941 "; aNbElem = "<<aNbElem<<
1942 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
1945 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1946 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
1947 TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
1948 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1950 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1951 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
1953 for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
1954 TInt aNodeId = aConnSlice[aConnId] - 1;
1955 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
1957 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1958 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
1962 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1963 aGaussCoordSlice[aDimId] /= aConnDim;
1969 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
1970 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
1971 INITMSG(MYVALUEDEBUG,"");
1972 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1973 TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
1974 ADDMSG(MYVALUEDEBUG,"{");
1975 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
1976 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
1978 ADDMSG(MYVALUEDEBUG,"} ");
1980 ADDMSG(MYVALUEDEBUG,std::endl);
1987 //---------------------------------------------------------------
1989 GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
1990 const TNodeInfo& theNodeInfo,
1991 TGaussCoord& theGaussCoord,
1992 const TElemNum& theElemNum,
1993 EModeSwitch theMode)
1995 INITMSG(MYDEBUG,"GetBaryCenter\n");
1996 const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
1997 TInt aDim = aMeshInfo->GetDim();
1998 static TInt aNbGauss = 1;
2000 bool anIsSubMesh = !theElemNum.empty();
2003 aNbElem = theElemNum.size();
2005 aNbElem = thePolygoneInfo.GetNbElem();
2007 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2011 "; aNbGauss = "<<aNbGauss<<
2012 "; aNbElem = "<<aNbElem<<
2013 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2016 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2017 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2019 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2020 TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
2021 TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
2022 TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
2024 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2025 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2027 for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2028 TInt aNodeId = aConnSlice[aConnId] - 1;
2029 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2031 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2032 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2036 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2037 aGaussCoordSlice[aDimId] /= aNbNodes;
2045 //---------------------------------------------------------------
2047 GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
2048 const TNodeInfo& theNodeInfo,
2049 TGaussCoord& theGaussCoord,
2050 const TElemNum& theElemNum,
2051 EModeSwitch theMode)
2053 INITMSG(MYDEBUG,"GetBaryCenter\n");
2054 const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
2055 TInt aDim = aMeshInfo->GetDim();
2056 static TInt aNbGauss = 1;
2058 bool anIsSubMesh = !theElemNum.empty();
2061 aNbElem = theElemNum.size();
2063 aNbElem = thePolyedreInfo.GetNbElem();
2065 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2069 "; aNbGauss = "<<aNbGauss<<
2070 "; aNbElem = "<<aNbElem<<
2071 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2074 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2075 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2077 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2078 TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
2079 TInt aNbFaces = aConnSliceArr.size();
2081 TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
2083 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2084 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2086 for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
2087 TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
2088 TInt aNbConn = aConnSlice.size();
2089 for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2090 TInt aNodeId = aConnSlice[aConnId] - 1;
2091 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2093 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2094 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2098 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2099 aGaussCoordSlice[aDimId] /= aNbNodes;