1 // Copyright (C) 2007-2014 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
22 #include "MED_GaussUtils.hxx"
23 #include "MED_Utilities.hxx"
26 static int MYDEBUG = 0;
27 static int MYVALUEDEBUG = 0;
29 // static int MYDEBUG = 0;
30 // static int MYVALUEDEBUG = 0;
33 //#define _DEBUG_REF_COORDS_
37 //---------------------------------------------------------------
40 TModeSwitchInfo(eFULL_INTERLACE),
50 ::Init(TInt theNbElem,
55 myModeSwitch = theMode;
58 myNbGauss = theNbGauss;
61 myGaussStep = myNbGauss*myDim;
63 myGaussCoord.resize(theNbElem*myGaussStep);
92 return (unsigned char*)&(myGaussCoord[0]);
98 ::GetCoordSliceArr(TInt theElemId) const
100 TCCoordSliceArr aCoordSliceArr(myNbGauss);
101 if(GetModeSwitch() == eFULL_INTERLACE){
102 TInt anId = theElemId*myGaussStep;
103 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
104 aCoordSliceArr[anGaussId] =
105 TCCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
110 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
111 aCoordSliceArr[anGaussId] =
112 TCCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
115 return aCoordSliceArr;
121 ::GetCoordSliceArr(TInt theElemId)
123 TCoordSliceArr aCoordSliceArr(myNbGauss);
124 if(GetModeSwitch() == eFULL_INTERLACE){
125 TInt anId = theElemId*myGaussStep;
126 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
127 aCoordSliceArr[anGaussId] =
128 TCoordSlice(myGaussCoord,std::slice(anId,myDim,1));
133 for(TInt anGaussId = 0; anGaussId < myNbGauss; anGaussId++){
134 aCoordSliceArr[anGaussId] =
135 TCoordSlice(myGaussCoord,std::slice(theElemId,myDim,myGaussStep));
138 return aCoordSliceArr;
142 //---------------------------------------------------------------
145 IsEqual(TFloat theLeft, TFloat theRight)
147 static TFloat EPS = 1.0E-3;
148 if(fabs(theLeft) + fabs(theRight) > EPS)
149 return fabs(theLeft-theRight)/(fabs(theLeft)+fabs(theRight)) < EPS;
154 //---------------------------------------------------------------
155 class TShapeFun::TFun
163 Init(TInt theNbGauss,
166 myFun.resize(theNbGauss*theNbRef);
171 GetFunSlice(TInt theGaussId) const
173 return TCFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
177 GetFunSlice(TInt theGaussId)
179 return TFloatVecSlice(myFun,std::slice(theGaussId*myNbRef,myNbRef,1));
183 //---------------------------------------------------------------
185 TShapeFun::TShapeFun(TInt theDim, TInt theNbRef):
186 myRefCoord(theNbRef*theDim),
192 TShapeFun::GetCoord(TInt theRefId) const
194 return TCCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
198 TShapeFun::GetCoord(TInt theRefId)
200 return TCoordSlice(myRefCoord,std::slice(theRefId*myDim,myDim,1));
204 TShapeFun::GetFun(const TCCoordSliceArr& theRef,
205 const TCCoordSliceArr& theGauss,
208 TInt aNbRef = theRef.size();
209 TInt aNbGauss = theGauss.size();
210 theFun.Init(aNbGauss,aNbRef);
214 TShapeFun::IsSatisfy(const TCCoordSliceArr& theRefCoord) const
216 TInt aNbRef = theRefCoord.size();
217 TInt aNbRef2 = GetNbRef();
218 INITMSG(MYDEBUG,"TShapeFun::IsSatisfy "<<
219 "- aNbRef("<<aNbRef<<")"<<
220 "; aNbRef2("<<aNbRef2<<")\n");
221 bool anIsSatisfy = (aNbRef == aNbRef2);
223 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
224 const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
225 TCCoordSlice aCoord = GetCoord(aRefId);
226 TInt aDim = aCoord.size();
227 bool anIsEqual = false;
228 for(TInt anId = 0; anId < aDim; anId++){
229 anIsEqual = IsEqual(aCoord[anId],aCoord2[anId]);
237 TCCoordSlice aCoord = GetCoord(aRefId);
238 INITMSG(MYDEBUG,aRefId + 1<<": aCoord = {");
239 TInt aDim = aCoord.size();
240 for(TInt anId = 0; anId < aDim; anId++)
241 ADDMSG(MYDEBUG,"\t"<<aCoord[anId]);
242 const TCCoordSlice& aCoord2 = theRefCoord[aRefId];
243 ADDMSG(MYDEBUG,"}\t!=\taCoord2 = {");
244 for(TInt anId = 0; anId < aDim; anId++)
245 ADDMSG(MYDEBUG,"\t"<<aCoord2[anId]);
246 ADDMSG(MYDEBUG,"}\n");
249 BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
256 BEGMSG(MYDEBUG,"anIsSatisfy = "<<anIsSatisfy<<"\n");
261 TShapeFun::Eval(const TCellInfo& theCellInfo,
262 const TNodeInfo& theNodeInfo,
263 const TElemNum& theElemNum,
264 const TCCoordSliceArr& theRef,
265 const TCCoordSliceArr& theGauss,
266 TGaussCoord& theGaussCoord,
269 INITMSG(MYDEBUG,"TShapeFun::Eval"<<std::endl);
271 if(IsSatisfy(theRef)){
272 const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
273 TInt aDim = aMeshInfo->GetDim();
274 TInt aNbGauss = theGauss.size();
276 bool anIsSubMesh = !theElemNum.empty();
279 aNbElem = theElemNum.size();
281 aNbElem = theCellInfo.GetNbElem();
283 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
286 InitFun(theRef,theGauss,aFun);
287 TInt aConnDim = theCellInfo.GetConnDim();
289 INITMSG(MYDEBUG,"aDim = "<<aDim<<
290 "; aNbGauss = "<<aNbGauss<<
291 "; aNbElem = "<<aNbElem<<
292 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
295 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
296 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
297 TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
298 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
300 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
301 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
302 TCFloatVecSlice aFunSlice = aFun.GetFunSlice(aGaussId);
304 for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
305 TInt aNodeId = aConnSlice[aConnId] - 1;
306 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
308 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
309 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId]*aFunSlice[aConnId];
317 INITMSG(MYVALUEDEBUG,"theGauss: ");
318 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
319 TCCoordSlice aCoordSlice = theGauss[aGaussId];
320 ADDMSG(MYVALUEDEBUG,"{");
321 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
322 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
324 ADDMSG(MYVALUEDEBUG,"} ");
326 ADDMSG(MYVALUEDEBUG,std::endl);
328 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
329 TCCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
330 INITMSG(MYVALUEDEBUG,"");
331 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
332 TCCoordSlice aCoordSlice = aCoordSliceArr[aGaussId];
333 ADDMSG(MYVALUEDEBUG,"{");
334 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
335 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
337 ADDMSG(MYVALUEDEBUG,"} ");
339 ADDMSG(MYVALUEDEBUG,std::endl);
349 //---------------------------------------------------------------
350 TSeg2a::TSeg2a():TShapeFun(1,2)
352 TInt aNbRef = GetNbRef();
353 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
354 TCoordSlice aCoord = GetCoord(aRefId);
356 case 0: aCoord[0] = -1.0; break;
357 case 1: aCoord[0] = 1.0; break;
363 TSeg2a::InitFun(const TCCoordSliceArr& theRef,
364 const TCCoordSliceArr& theGauss,
367 GetFun(theRef,theGauss,theFun);
369 TInt aNbGauss = theGauss.size();
370 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
371 const TCCoordSlice& aCoord = theGauss[aGaussId];
372 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
374 aSlice[0] = 0.5*(1.0 - aCoord[0]);
375 aSlice[1] = 0.5*(1.0 + aCoord[0]);
380 //---------------------------------------------------------------
381 TSeg3a::TSeg3a():TShapeFun(1,3)
383 TInt aNbRef = GetNbRef();
384 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
385 TCoordSlice aCoord = GetCoord(aRefId);
387 case 0: aCoord[0] = -1.0; break;
388 case 1: aCoord[0] = 1.0; break;
389 case 2: aCoord[0] = 0.0; break;
395 TSeg3a::InitFun(const TCCoordSliceArr& theRef,
396 const TCCoordSliceArr& theGauss,
399 GetFun(theRef,theGauss,theFun);
401 TInt aNbGauss = theGauss.size();
402 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
403 const TCCoordSlice& aCoord = theGauss[aGaussId];
404 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
406 aSlice[0] = 0.5*(1.0 - aCoord[0])*aCoord[0];
407 aSlice[1] = 0.5*(1.0 + aCoord[0])*aCoord[0];
408 aSlice[2] = (1.0 + aCoord[0])*(1.0 - aCoord[0]);
414 //---------------------------------------------------------------
418 TInt aNbRef = GetNbRef();
419 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
420 TCoordSlice aCoord = GetCoord(aRefId);
422 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
423 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
424 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
430 TTria3a::InitFun(const TCCoordSliceArr& theRef,
431 const TCCoordSliceArr& theGauss,
434 GetFun(theRef,theGauss,theFun);
436 TInt aNbGauss = theGauss.size();
437 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
438 const TCCoordSlice& aCoord = theGauss[aGaussId];
439 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
441 aSlice[0] = 0.5*(1.0 + aCoord[1]);
442 aSlice[1] = -0.5*(aCoord[0] + aCoord[1]);
443 aSlice[2] = 0.5*(1.0 + aCoord[0]);
449 //---------------------------------------------------------------
450 TTria6a::TTria6a():TShapeFun(2,6)
452 TInt aNbRef = GetNbRef();
453 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
454 TCoordSlice aCoord = GetCoord(aRefId);
456 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
457 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
458 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
460 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
461 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
462 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
468 TTria6a::InitFun(const TCCoordSliceArr& theRef,
469 const TCCoordSliceArr& theGauss,
472 GetFun(theRef,theGauss,theFun);
474 TInt aNbGauss = theGauss.size();
475 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
476 const TCCoordSlice& aCoord = theGauss[aGaussId];
477 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
479 aSlice[0] = 0.5*(1.0 + aCoord[1])*aCoord[1];
480 aSlice[1] = 0.5*(aCoord[0] + aCoord[1])*(aCoord[0] + aCoord[1] + 1);
481 aSlice[2] = 0.5*(1.0 + aCoord[0])*aCoord[0];
483 aSlice[3] = -1.0*(1.0 + aCoord[1])*(aCoord[0] + aCoord[1]);
484 aSlice[4] = -1.0*(1.0 + aCoord[0])*(aCoord[0] + aCoord[1]);
485 aSlice[5] = (1.0 + aCoord[1])*(1.0 + aCoord[1]);
491 //---------------------------------------------------------------
495 TInt aNbRef = GetNbRef();
496 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
497 TCoordSlice aCoord = GetCoord(aRefId);
499 case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
500 case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
501 case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
507 TTria3b::InitFun(const TCCoordSliceArr& theRef,
508 const TCCoordSliceArr& theGauss,
511 GetFun(theRef,theGauss,theFun);
513 TInt aNbGauss = theGauss.size();
514 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
515 const TCCoordSlice& aCoord = theGauss[aGaussId];
516 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
518 aSlice[0] = 1.0 - aCoord[0] - aCoord[1];
519 aSlice[1] = aCoord[0];
520 aSlice[2] = aCoord[1];
526 //---------------------------------------------------------------
530 TInt aNbRef = GetNbRef();
531 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
532 TCoordSlice aCoord = GetCoord(aRefId);
534 case 0: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
535 case 1: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
536 case 2: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
538 case 3: aCoord[0] = 0.5; aCoord[1] = 0.0; break;
539 case 4: aCoord[0] = 0.5; aCoord[1] = 0.5; break;
540 case 5: aCoord[0] = 0.0; aCoord[1] = 0.5; break;
546 TTria6b::InitFun(const TCCoordSliceArr& theRef,
547 const TCCoordSliceArr& theGauss,
550 GetFun(theRef,theGauss,theFun);
552 TInt aNbGauss = theGauss.size();
553 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
554 const TCCoordSlice& aCoord = theGauss[aGaussId];
555 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
557 aSlice[0] = (1.0 - aCoord[0] - aCoord[1])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1]);
558 aSlice[1] = aCoord[0]*(2.0*aCoord[0] - 1.0);
559 aSlice[2] = aCoord[1]*(2.0*aCoord[1] - 1.0);
561 aSlice[3] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1]);
562 aSlice[4] = 4.0*aCoord[0]*aCoord[1];
563 aSlice[5] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1]);
569 //---------------------------------------------------------------
573 TInt aNbRef = GetNbRef();
574 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
575 TCoordSlice aCoord = GetCoord(aRefId);
577 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
578 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
579 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
580 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
586 TQuad4a::InitFun(const TCCoordSliceArr& theRef,
587 const TCCoordSliceArr& theGauss,
590 GetFun(theRef,theGauss,theFun);
592 TInt aNbGauss = theGauss.size();
593 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
594 const TCCoordSlice& aCoord = theGauss[aGaussId];
595 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
597 aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0]);
598 aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0]);
599 aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0]);
600 aSlice[3] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
606 //---------------------------------------------------------------
610 TInt aNbRef = GetNbRef();
611 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
612 TCoordSlice aCoord = GetCoord(aRefId);
614 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
615 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
616 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
617 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
619 case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
620 case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
621 case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
622 case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
628 TQuad8a::InitFun(const TCCoordSliceArr& theRef,
629 const TCCoordSliceArr& theGauss,
632 GetFun(theRef,theGauss,theFun);
634 TInt aNbGauss = theGauss.size();
635 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
636 const TCCoordSlice& aCoord = theGauss[aGaussId];
637 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
639 aSlice[0] = 0.25*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(aCoord[1] - aCoord[0] - 1.0);
640 aSlice[1] = 0.25*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(-aCoord[1] - aCoord[0] - 1.0);
641 aSlice[2] = 0.25*(1.0 - aCoord[1])*(1.0 + aCoord[0])*(-aCoord[1] + aCoord[0] - 1.0);
642 aSlice[3] = 0.25*(1.0 + aCoord[1])*(1.0 + aCoord[0])*(aCoord[1] + aCoord[0] - 1.0);
644 aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
645 aSlice[5] = 0.5*(1.0 - aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
646 aSlice[6] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[1]);
647 aSlice[7] = 0.5*(1.0 + aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[0]);
653 //---------------------------------------------------------------
657 TInt aNbRef = GetNbRef();
658 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
659 TCoordSlice aCoord = GetCoord(aRefId);
661 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
662 case 1: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
663 case 2: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
664 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
666 case 4: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
667 case 5: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
668 case 6: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
669 case 7: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
671 case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
677 TQuad9a::InitFun(const TCCoordSliceArr& theRef,
678 const TCCoordSliceArr& theGauss,
681 GetFun(theRef,theGauss,theFun);
683 TInt aNbGauss = theGauss.size();
684 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
685 const TCCoordSlice& aCoord = theGauss[aGaussId];
686 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
688 aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
689 aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
690 aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
691 aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
693 aSlice[4] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
694 aSlice[5] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
695 aSlice[6] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
696 aSlice[7] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
698 aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
704 //---------------------------------------------------------------
708 TInt aNbRef = GetNbRef();
709 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
710 TCoordSlice aCoord = GetCoord(aRefId);
712 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
713 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
714 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
715 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
721 TQuad4b::InitFun(const TCCoordSliceArr& theRef,
722 const TCCoordSliceArr& theGauss,
725 GetFun(theRef,theGauss,theFun);
727 TInt aNbGauss = theGauss.size();
728 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
729 const TCCoordSlice& aCoord = theGauss[aGaussId];
730 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
732 aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
733 aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
734 aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
735 aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
741 //---------------------------------------------------------------
745 TInt aNbRef = GetNbRef();
746 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
747 TCoordSlice aCoord = GetCoord(aRefId);
749 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
750 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
751 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
752 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
754 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
755 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
756 case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
757 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
763 TQuad8b::InitFun(const TCCoordSliceArr& theRef,
764 const TCCoordSliceArr& theGauss,
767 GetFun(theRef,theGauss,theFun);
769 TInt aNbGauss = theGauss.size();
770 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
771 const TCCoordSlice& aCoord = theGauss[aGaussId];
772 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
774 aSlice[0] = 0.25*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(-1.0 - aCoord[0] - aCoord[1]);
775 aSlice[1] = 0.25*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(-1.0 + aCoord[0] - aCoord[1]);
776 aSlice[2] = 0.25*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(-1.0 + aCoord[0] + aCoord[1]);
777 aSlice[3] = 0.25*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(-1.0 - aCoord[0] + aCoord[1]);
779 aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]);
780 aSlice[5] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0]);
781 aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1]);
782 aSlice[7] = 0.5*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0]);
784 //aSlice[4] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
785 //aSlice[5] = 0.5*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
786 //aSlice[6] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
787 //aSlice[7] = 0.5*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[1]);
793 //---------------------------------------------------------------
797 TInt aNbRef = GetNbRef();
798 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
799 TCoordSlice aCoord = GetCoord(aRefId);
801 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; break;
802 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; break;
803 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; break;
804 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; break;
806 case 4: aCoord[0] = 0.0; aCoord[1] = -1.0; break;
807 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; break;
808 case 6: aCoord[0] = 0.0; aCoord[1] = 1.0; break;
809 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; break;
811 case 8: aCoord[0] = 0.0; aCoord[1] = 0.0; break;
817 TQuad9b::InitFun(const TCCoordSliceArr& theRef,
818 const TCCoordSliceArr& theGauss,
821 GetFun(theRef,theGauss,theFun);
823 TInt aNbGauss = theGauss.size();
824 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
825 const TCCoordSlice& aCoord = theGauss[aGaussId];
826 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
828 aSlice[0] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] - 1.0);
829 aSlice[1] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] - 1.0);
830 aSlice[2] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] + 1.0)*(aCoord[1] + 1.0);
831 aSlice[3] = 0.25*aCoord[0]*aCoord[1]*(aCoord[0] - 1.0)*(aCoord[1] + 1.0);
833 aSlice[4] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0);
834 aSlice[5] = 0.5*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1]);
835 aSlice[6] = 0.5*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0);
836 aSlice[7] = 0.5*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1]);
838 aSlice[8] = (1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1]);
844 //---------------------------------------------------------------
845 TTetra4a::TTetra4a():
848 TInt aNbRef = GetNbRef();
849 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
850 TCoordSlice aCoord = GetCoord(aRefId);
852 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
853 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
854 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
855 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
861 TTetra4a::InitFun(const TCCoordSliceArr& theRef,
862 const TCCoordSliceArr& theGauss,
865 GetFun(theRef,theGauss,theFun);
867 TInt aNbGauss = theGauss.size();
868 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
869 const TCCoordSlice& aCoord = theGauss[aGaussId];
870 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
872 aSlice[0] = aCoord[1];
873 aSlice[1] = aCoord[2];
874 aSlice[2] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
875 aSlice[3] = aCoord[0];
881 //---------------------------------------------------------------
882 TTetra10a::TTetra10a():
885 TInt aNbRef = GetNbRef();
886 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
887 TCoordSlice aCoord = GetCoord(aRefId);
889 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
890 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
891 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
892 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
894 case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
895 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
896 case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
898 case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
899 case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
900 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
906 TTetra10a::InitFun(const TCCoordSliceArr& theRef,
907 const TCCoordSliceArr& theGauss,
910 GetFun(theRef,theGauss,theFun);
912 TInt aNbGauss = theGauss.size();
913 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
914 const TCCoordSlice& aCoord = theGauss[aGaussId];
915 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
917 aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
918 aSlice[1] = aCoord[2]*(2.0*aCoord[2] - 1.0);
919 aSlice[2] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
920 aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
922 aSlice[4] = 4.0*aCoord[1]*aCoord[2];
923 aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
924 aSlice[6] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
926 aSlice[7] = 4.0*aCoord[0]*aCoord[1];
927 aSlice[8] = 4.0*aCoord[0]*aCoord[2];
928 aSlice[9] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
934 //---------------------------------------------------------------
937 TTetra4b::TTetra4b():
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;
953 TTetra4b::InitFun(const TCCoordSliceArr& theRef,
954 const TCCoordSliceArr& theGauss,
957 GetFun(theRef,theGauss,theFun);
959 TInt aNbGauss = theGauss.size();
960 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
961 const TCCoordSlice& aCoord = theGauss[aGaussId];
962 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
964 aSlice[0] = aCoord[1];
965 aSlice[2] = aCoord[2];
966 aSlice[1] = 1.0 - aCoord[0] - aCoord[1] - aCoord[2];
967 aSlice[3] = aCoord[0];
973 //---------------------------------------------------------------
974 TTetra10b::TTetra10b():
977 TInt aNbRef = GetNbRef();
978 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
979 TCoordSlice aCoord = GetCoord(aRefId);
981 case 0: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
982 case 2: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
983 case 1: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
984 case 3: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
986 case 6: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
987 case 5: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
988 case 4: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
990 case 7: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
991 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
992 case 8: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
998 TTetra10b::InitFun(const TCCoordSliceArr& theRef,
999 const TCCoordSliceArr& theGauss,
1002 GetFun(theRef,theGauss,theFun);
1004 TInt aNbGauss = theGauss.size();
1005 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1006 const TCCoordSlice& aCoord = theGauss[aGaussId];
1007 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1009 aSlice[0] = aCoord[1]*(2.0*aCoord[1] - 1.0);
1010 aSlice[2] = aCoord[2]*(2.0*aCoord[2] - 1.0);
1011 aSlice[1] = (1.0 - aCoord[0] - aCoord[1] - aCoord[2])*(1.0 - 2.0*aCoord[0] - 2.0*aCoord[1] - 2.0*aCoord[2]);
1012 aSlice[3] = aCoord[0]*(2.0*aCoord[0] - 1.0);
1014 aSlice[6] = 4.0*aCoord[1]*aCoord[2];
1015 aSlice[5] = 4.0*aCoord[2]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1016 aSlice[4] = 4.0*aCoord[1]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1018 aSlice[7] = 4.0*aCoord[0]*aCoord[1];
1019 aSlice[9] = 4.0*aCoord[0]*aCoord[2];
1020 aSlice[8] = 4.0*aCoord[0]*(1.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1026 //---------------------------------------------------------------
1030 TInt aNbRef = GetNbRef();
1031 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1032 TCoordSlice aCoord = GetCoord(aRefId);
1034 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1035 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1036 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1037 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1038 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1039 case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1040 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1041 case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1047 THexa8a::InitFun(const TCCoordSliceArr& theRef,
1048 const TCCoordSliceArr& theGauss,
1051 GetFun(theRef,theGauss,theFun);
1053 TInt aNbGauss = theGauss.size();
1054 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1055 const TCCoordSlice& aCoord = theGauss[aGaussId];
1056 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1058 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1059 aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1060 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1061 aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1063 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1064 aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1065 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1066 aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1071 //---------------------------------------------------------------
1072 THexa20a::THexa20a(TInt theDim, TInt theNbRef):
1073 TShapeFun(theDim,theNbRef)
1075 TInt aNbRef = myRefCoord.size();
1076 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1077 TCoordSlice aCoord = GetCoord(aRefId);
1079 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1080 case 1: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1081 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1082 case 3: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1083 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1084 case 5: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1085 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1086 case 7: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1088 case 8: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1089 case 9: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1090 case 10: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1091 case 11: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1092 case 12: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1093 case 13: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1094 case 14: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1095 case 15: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1096 case 16: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1097 case 17: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1098 case 18: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1099 case 19: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1105 THexa20a::InitFun(const TCCoordSliceArr& theRef,
1106 const TCCoordSliceArr& theGauss,
1109 GetFun(theRef,theGauss,theFun);
1111 TInt aNbGauss = theGauss.size();
1112 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1113 const TCCoordSlice& aCoord = theGauss[aGaussId];
1114 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1116 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1117 (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1118 aSlice[1] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1119 (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1120 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1121 (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1122 aSlice[3] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1123 (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1124 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1125 (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1126 aSlice[5] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1127 (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1128 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1129 (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1130 aSlice[7] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1131 (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1133 aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1134 aSlice[9] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1135 aSlice[10] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1136 aSlice[11] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1137 aSlice[12] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1138 aSlice[13] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1139 aSlice[14] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1140 aSlice[15] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1141 aSlice[16] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1142 aSlice[17] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1143 aSlice[18] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1144 aSlice[19] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1150 //---------------------------------------------------------------
1151 THexa27a::THexa27a():
1154 TInt aNbRef = myRefCoord.size();
1155 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1156 TCoordSlice aCoord = GetCoord(aRefId);
1158 case 20: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1159 case 21: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1160 case 22: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1161 case 23: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1162 case 24: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1163 case 25: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1164 case 26: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1170 THexa27a::InitFun(const TCCoordSliceArr& theRef,
1171 const TCCoordSliceArr& theGauss,
1174 GetFun(theRef,theGauss,theFun);
1176 TInt aNbGauss = theGauss.size();
1177 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1178 const TCCoordSlice& aCoord = theGauss[aGaussId];
1179 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1181 aSlice[0] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1182 aSlice[1] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1183 aSlice[2] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1184 aSlice[3] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1185 aSlice[4] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1186 aSlice[5] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1187 aSlice[6] = 0.125*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1188 aSlice[7] = 0.125*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1190 aSlice[8] = 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] - 1.0);
1191 aSlice[9] = 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1192 aSlice[10]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1193 aSlice[11]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1194 aSlice[12]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1195 aSlice[13]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1196 aSlice[14]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1197 aSlice[15]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1198 aSlice[16]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1199 aSlice[17]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1200 aSlice[18]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*aCoord[2]*(aCoord[2] + 1.0);
1201 aSlice[19]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1202 aSlice[20]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] - 1.0);
1203 aSlice[21]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] - 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1204 aSlice[22]= 0.25*aCoord[0]*(aCoord[0] + 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1205 aSlice[23]= 0.25*(1.0 - aCoord[0]*aCoord[0])*aCoord[1]*(aCoord[1] + 1.0)*(1.0 - aCoord[2]*aCoord[2]);
1206 aSlice[24]= 0.25*aCoord[0]*(aCoord[0] - 1.0)*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[2]*aCoord[2]);
1207 aSlice[25]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*aCoord[2]*(aCoord[2] + 1.0);
1208 aSlice[26]= 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[1]*aCoord[1]);
1214 //---------------------------------------------------------------
1218 TInt aNbRef = GetNbRef();
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;
1235 THexa8b::InitFun(const TCCoordSliceArr& theRef,
1236 const TCCoordSliceArr& theGauss,
1239 GetFun(theRef,theGauss,theFun);
1241 TInt aNbGauss = theGauss.size();
1242 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1243 const TCCoordSlice& aCoord = theGauss[aGaussId];
1244 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1246 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1247 aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1248 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1249 aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1251 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1252 aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1253 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1254 aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1260 //---------------------------------------------------------------
1261 THexa20b::THexa20b(TInt theDim, TInt theNbRef):
1262 TShapeFun(theDim,theNbRef)
1264 TInt aNbRef = myRefCoord.size();
1265 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1266 TCoordSlice aCoord = GetCoord(aRefId);
1268 case 0: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1269 case 3: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1270 case 2: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1271 case 1: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1272 case 4: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1273 case 7: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1274 case 6: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1275 case 5: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1277 case 11: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = -1.0; break;
1278 case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1279 case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = -1.0; break;
1280 case 8: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = -1.0; break;
1281 case 16: aCoord[0] = -1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1282 case 19: aCoord[0] = 1.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1283 case 18: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1284 case 17: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1285 case 15: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 1.0; break;
1286 case 14: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1287 case 13: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 1.0; break;
1288 case 12: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1294 THexa20b::InitFun(const TCCoordSliceArr& theRef,
1295 const TCCoordSliceArr& theGauss,
1298 GetFun(theRef,theGauss,theFun);
1300 TInt aNbGauss = theGauss.size();
1301 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1302 const TCCoordSlice& aCoord = theGauss[aGaussId];
1303 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1305 aSlice[0] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1306 (-2.0 - aCoord[0] - aCoord[1] - aCoord[2]);
1307 aSlice[3] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2])*
1308 (-2.0 + aCoord[0] - aCoord[1] - aCoord[2]);
1309 aSlice[2] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1310 (-2.0 + aCoord[0] + aCoord[1] - aCoord[2]);
1311 aSlice[1] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2])*
1312 (-2.0 - aCoord[0] + aCoord[1] - aCoord[2]);
1313 aSlice[4] = 0.125*(1.0 - aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1314 (-2.0 - aCoord[0] - aCoord[1] + aCoord[2]);
1315 aSlice[7] = 0.125*(1.0 + aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2])*
1316 (-2.0 + aCoord[0] - aCoord[1] + aCoord[2]);
1317 aSlice[6] = 0.125*(1.0 + aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1318 (-2.0 + aCoord[0] + aCoord[1] + aCoord[2]);
1319 aSlice[5] = 0.125*(1.0 - aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2])*
1320 (-2.0 - aCoord[0] + aCoord[1] + aCoord[2]);
1322 aSlice[11] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 - aCoord[2]);
1323 aSlice[10] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 - aCoord[2]);
1324 aSlice[9] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 - aCoord[2]);
1325 aSlice[8] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 - aCoord[2]);
1326 aSlice[16] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 - aCoord[1]);
1327 aSlice[19] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 - aCoord[1]);
1328 aSlice[18] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 + aCoord[0])*(1.0 + aCoord[1]);
1329 aSlice[17] = 0.25*(1.0 - aCoord[2]*aCoord[2])*(1.0 - aCoord[0])*(1.0 + aCoord[1]);
1330 aSlice[15] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 - aCoord[1])*(1.0 + aCoord[2]);
1331 aSlice[14] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 + aCoord[0])*(1.0 + aCoord[2]);
1332 aSlice[13] = 0.25*(1.0 - aCoord[0]*aCoord[0])*(1.0 + aCoord[1])*(1.0 + aCoord[2]);
1333 aSlice[12] = 0.25*(1.0 - aCoord[1]*aCoord[1])*(1.0 - aCoord[0])*(1.0 + aCoord[2]);
1339 //---------------------------------------------------------------
1340 TPenta6a::TPenta6a():
1343 TInt aNbRef = myRefCoord.size();
1344 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1345 TCoordSlice aCoord = GetCoord(aRefId);
1347 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1348 case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1349 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1350 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1351 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1352 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1358 TPenta6a::InitFun(const TCCoordSliceArr& theRef,
1359 const TCCoordSliceArr& theGauss,
1362 GetFun(theRef,theGauss,theFun);
1364 TInt aNbGauss = theGauss.size();
1365 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1366 const TCCoordSlice& aCoord = theGauss[aGaussId];
1367 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1369 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0]);
1370 aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1371 aSlice[2] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1373 aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1374 aSlice[4] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1375 aSlice[5] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1381 //---------------------------------------------------------------
1382 TPenta6b::TPenta6b():
1385 TInt aNbRef = myRefCoord.size();
1386 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1387 TCoordSlice aCoord = GetCoord(aRefId);
1389 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1390 case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1391 case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1392 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1393 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1394 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1400 TPenta6b::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]);
1412 aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0]);
1413 aSlice[1] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1415 aSlice[3] = 0.5*aCoord[1]*(aCoord[0] + 1.0);
1416 aSlice[5] = 0.5*aCoord[2]*(aCoord[0] + 1.0);
1417 aSlice[4] = 0.5*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1423 //---------------------------------------------------------------
1424 TPenta15a::TPenta15a():
1427 TInt aNbRef = myRefCoord.size();
1428 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1429 TCoordSlice aCoord = GetCoord(aRefId);
1431 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1432 case 1: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1433 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1434 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1435 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1436 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1438 case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1439 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1440 case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1441 case 9: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1442 case 10: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1443 case 11: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1444 case 12: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1445 case 13: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1446 case 14: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1452 TPenta15a::InitFun(const TCCoordSliceArr& theRef,
1453 const TCCoordSliceArr& theGauss,
1456 GetFun(theRef,theGauss,theFun);
1458 TInt aNbGauss = theGauss.size();
1459 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1460 const TCCoordSlice& aCoord = theGauss[aGaussId];
1461 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1463 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1464 aSlice[1] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1465 aSlice[2] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1467 aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1468 aSlice[4] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1469 aSlice[5] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1471 aSlice[6] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1472 aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1473 aSlice[8] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1475 aSlice[9] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1476 aSlice[10] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1477 aSlice[11] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1479 aSlice[12] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1480 aSlice[13] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1481 aSlice[14] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1487 //---------------------------------------------------------------
1488 TPenta15b::TPenta15b():
1491 TInt aNbRef = myRefCoord.size();
1492 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1493 TCoordSlice aCoord = GetCoord(aRefId);
1495 case 0: aCoord[0] = -1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1496 case 2: aCoord[0] = -1.0; aCoord[1] = -0.0; aCoord[2] = 1.0; break;
1497 case 1: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1498 case 3: aCoord[0] = 1.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1499 case 5: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1500 case 4: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1502 case 8: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1503 case 7: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1504 case 6: aCoord[0] = -1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1505 case 12: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1506 case 14: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1507 case 13: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1508 case 11: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1509 case 10: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1510 case 9: aCoord[0] = 1.0; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1516 TPenta15b::InitFun(const TCCoordSliceArr& theRef,
1517 const TCCoordSliceArr& theGauss,
1520 GetFun(theRef,theGauss,theFun);
1522 TInt aNbGauss = theGauss.size();
1523 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1524 const TCCoordSlice& aCoord = theGauss[aGaussId];
1525 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1527 aSlice[0] = 0.5*aCoord[1]*(1.0 - aCoord[0])*(2.0*aCoord[1] - 2.0 - aCoord[0]);
1528 aSlice[2] = 0.5*aCoord[2]*(1.0 - aCoord[0])*(2.0*aCoord[2] - 2.0 - aCoord[0]);
1529 aSlice[1] = 0.5*(aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1531 aSlice[3] = 0.5*aCoord[1]*(1.0 + aCoord[0])*(2.0*aCoord[1] - 2.0 + aCoord[0]);
1532 aSlice[5] = 0.5*aCoord[2]*(1.0 + aCoord[0])*(2.0*aCoord[2] - 2.0 + aCoord[0]);
1533 aSlice[4] = 0.5*(-aCoord[0] - 1.0)*(1.0 - aCoord[1] - aCoord[2])*(-aCoord[0] + 2.0*aCoord[1] + 2.0*aCoord[2]);
1535 aSlice[8] = 2.0*aCoord[1]*aCoord[2]*(1.0 - aCoord[0]);
1536 aSlice[7] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1537 aSlice[6] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]);
1539 aSlice[12] = aCoord[1]*(1.0 - aCoord[0]*aCoord[0]);
1540 aSlice[14] = aCoord[2]*(1.0 - aCoord[0]*aCoord[0]);
1541 aSlice[13] = (1.0 - aCoord[1] - aCoord[2])*(1.0 - aCoord[0]*aCoord[0]);
1543 aSlice[11] = 2.0*aCoord[1]*aCoord[2]*(1.0 + aCoord[0]);
1544 aSlice[10] = 2.0*aCoord[2]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1545 aSlice[9] = 2.0*aCoord[1]*(1.0 - aCoord[1] - aCoord[2])*(1.0 + aCoord[0]);
1551 //---------------------------------------------------------------
1555 TInt aNbRef = myRefCoord.size();
1556 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1557 TCoordSlice aCoord = GetCoord(aRefId);
1559 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1560 case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1561 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1562 case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1563 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1569 TPyra5a::InitFun(const TCCoordSliceArr& theRef,
1570 const TCCoordSliceArr& theGauss,
1573 GetFun(theRef,theGauss,theFun);
1575 TInt aNbGauss = theGauss.size();
1576 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1577 const TCCoordSlice& aCoord = theGauss[aGaussId];
1578 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1580 // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1581 // Seems shape function for ePYRA5 elements is:
1582 // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1583 // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1584 // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1585 // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1587 aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1588 aSlice[1] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1589 aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1590 aSlice[3] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1591 aSlice[4] = aCoord[2];
1597 //---------------------------------------------------------------
1601 TInt aNbRef = myRefCoord.size();
1602 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1603 TCoordSlice aCoord = GetCoord(aRefId);
1605 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1606 case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1607 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1608 case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1609 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1615 TPyra5b::InitFun(const TCCoordSliceArr& theRef,
1616 const TCCoordSliceArr& theGauss,
1619 GetFun(theRef,theGauss,theFun);
1621 TInt aNbGauss = theGauss.size();
1622 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1623 const TCCoordSlice& aCoord = theGauss[aGaussId];
1624 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1626 // Fix for 0019920: EDF 788 VISU: Bad localisation of gauss points for pyramids
1627 // Seems shape function for ePYRA5 elements is:
1628 // w0 = 0.25*(-X + Y - 1.0)*(-X - Y - 1.0)*(1.0 - Z);
1629 // w1 = 0.25*(-X - Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1630 // w2 = 0.25*(+X + Y - 1.0)*(+X - Y - 1.0)*(1.0 - Z);
1631 // w3 = 0.25*(+X + Y - 1.0)*(-X + Y - 1.0)*(1.0 - Z);
1633 aSlice[0] = 0.25*(-aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1634 aSlice[3] = 0.25*(-aCoord[0] - aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1635 aSlice[2] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(+aCoord[0] - aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1636 aSlice[1] = 0.25*(+aCoord[0] + aCoord[1] - 1.0)*(-aCoord[0] + aCoord[1] - 1.0)*(1.0 - aCoord[2]);
1637 aSlice[4] = aCoord[2];
1643 //---------------------------------------------------------------
1644 TPyra13a::TPyra13a():
1647 TInt aNbRef = myRefCoord.size();
1648 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1649 TCoordSlice aCoord = GetCoord(aRefId);
1651 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1652 case 1: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1653 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1654 case 3: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1655 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1657 case 5: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1658 case 6: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1659 case 7: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1660 case 8: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1661 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1662 case 10: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1663 case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1664 case 12: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
1670 TPyra13a::InitFun(const TCCoordSliceArr& theRef,
1671 const TCCoordSliceArr& theGauss,
1674 GetFun(theRef,theGauss,theFun);
1676 TInt aNbGauss = theGauss.size();
1677 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1678 const TCCoordSlice& aCoord = theGauss[aGaussId];
1679 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1681 aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1682 (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1683 aSlice[1] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1684 (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1685 aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1686 (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1687 aSlice[3] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1688 (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1690 aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1692 aSlice[5] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1693 (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1694 aSlice[6] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1695 (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1696 aSlice[7] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1697 (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1698 aSlice[8] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1699 (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1701 aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1703 aSlice[10] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1705 aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1707 aSlice[12] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1714 //---------------------------------------------------------------
1715 TPyra13b::TPyra13b():
1718 TInt aNbRef = myRefCoord.size();
1719 for(TInt aRefId = 0; aRefId < aNbRef; aRefId++){
1720 TCoordSlice aCoord = GetCoord(aRefId);
1722 case 0: aCoord[0] = 1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1723 case 3: aCoord[0] = 0.0; aCoord[1] = 1.0; aCoord[2] = 0.0; break;
1724 case 2: aCoord[0] = -1.0; aCoord[1] = 0.0; aCoord[2] = 0.0; break;
1725 case 1: aCoord[0] = 0.0; aCoord[1] = -1.0; aCoord[2] = 0.0; break;
1726 case 4: aCoord[0] = 0.0; aCoord[1] = 0.0; aCoord[2] = 1.0; break;
1728 case 8: aCoord[0] = 0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1729 case 7: aCoord[0] = -0.5; aCoord[1] = 0.5; aCoord[2] = 0.0; break;
1730 case 6: aCoord[0] = -0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1731 case 5: aCoord[0] = 0.5; aCoord[1] = -0.5; aCoord[2] = 0.0; break;
1732 case 9: aCoord[0] = 0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1733 case 12: aCoord[0] = 0.0; aCoord[1] = 0.5; aCoord[2] = 0.5; break;
1734 case 11: aCoord[0] = -0.5; aCoord[1] = 0.0; aCoord[2] = 0.5; break;
1735 case 10: aCoord[0] = 0.0; aCoord[1] = -0.5; aCoord[2] = 0.5; break;
1741 TPyra13b::InitFun(const TCCoordSliceArr& theRef,
1742 const TCCoordSliceArr& theGauss,
1745 GetFun(theRef,theGauss,theFun);
1747 TInt aNbGauss = theGauss.size();
1748 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
1749 const TCCoordSlice& aCoord = theGauss[aGaussId];
1750 TFloatVecSlice aSlice = theFun.GetFunSlice(aGaussId);
1752 aSlice[0] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1753 (aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1754 aSlice[3] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1755 (aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1756 aSlice[2] = 0.5*(+aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1757 (-aCoord[0] - 0.5)/(1.0 - aCoord[2]);
1758 aSlice[1] = 0.5*(+aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1759 (-aCoord[1] - 0.5)/(1.0 - aCoord[2]);
1761 aSlice[4] = 2.0*aCoord[2]*(aCoord[2] - 0.5);
1763 aSlice[8] = 0.5*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1764 (aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1765 aSlice[7] = 0.5*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*
1766 (aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1767 aSlice[6] = 0.5*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1768 (-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1769 aSlice[5] = 0.5*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*
1770 (-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/(1.0 - aCoord[2]);
1772 aSlice[9] = 0.5*aCoord[2]*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1774 aSlice[12] = 0.5*aCoord[2]*(-aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)/
1776 aSlice[11] = 0.5*aCoord[2]*(aCoord[0] - aCoord[1] + aCoord[2] - 1.0)*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1778 aSlice[10] = 0.5*aCoord[2]*(aCoord[0] + aCoord[1] + aCoord[2] - 1.0)*(-aCoord[0] + aCoord[1] + aCoord[2] - 1.0)/
1785 //---------------------------------------------------------------
1787 GetGaussCoord3D(const TGaussInfo& theGaussInfo,
1788 const TCellInfo& theCellInfo,
1789 const TNodeInfo& theNodeInfo,
1790 TGaussCoord& theGaussCoord,
1791 const TElemNum& theElemNum,
1792 EModeSwitch theMode)
1794 INITMSG(MYDEBUG,"GetGaussCoord3D\n");
1796 if(theGaussInfo.myGeom == theCellInfo.myGeom){
1797 EGeometrieElement aGeom = theGaussInfo.myGeom;
1799 TInt aNbRef = theGaussInfo.GetNbRef();
1800 TCCoordSliceArr aRefSlice(aNbRef);
1801 for(TInt anId = 0; anId < aNbRef; anId++)
1802 aRefSlice[anId] = theGaussInfo.GetRefCoordSlice(anId);
1804 TInt aNbGauss = theGaussInfo.GetNbGauss();
1805 TCCoordSliceArr aGaussSlice(aNbGauss);
1806 for(TInt anId = 0; anId < aNbGauss; anId++)
1807 aGaussSlice[anId] = theGaussInfo.GetGaussCoordSlice(anId);
1811 INITMSG(MYDEBUG,"eSEG2"<<std::endl);
1813 if(TSeg2a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1819 INITMSG(MYDEBUG,"eSEG3"<<std::endl);
1821 if(TSeg3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1827 INITMSG(MYDEBUG,"eTRIA3"<<std::endl);
1829 if(TTria3a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1832 if(TTria3b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1838 INITMSG(MYDEBUG,"eTRIA6"<<std::endl);
1840 if(TTria6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1843 if(TTria6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1849 INITMSG(MYDEBUG,"eQUAD4"<<std::endl);
1851 if(TQuad4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1854 if(TQuad4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1860 INITMSG(MYDEBUG,"eQUAD8"<<std::endl);
1862 if(TQuad8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1865 if(TQuad8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1871 INITMSG(MYDEBUG,"eQUAD9"<<std::endl);
1873 if(TQuad9a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1876 if(TQuad9b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1882 INITMSG(MYDEBUG,"eTETRA4"<<std::endl);
1884 if(TTetra4a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1887 if(TTetra4b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1893 INITMSG(MYDEBUG,"ePYRA5"<<std::endl);
1895 if(TPyra5a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1898 if(TPyra5b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1904 INITMSG(MYDEBUG,"ePENTA6"<<std::endl);
1906 if(TPenta6a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1909 if(TPenta6b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1915 INITMSG(MYDEBUG,"eHEXA8"<<std::endl);
1917 if(THexa8a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1920 if(THexa8b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1926 INITMSG(MYDEBUG,"eTETRA10"<<std::endl);
1928 if(TTetra10a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1931 if(TTetra10b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1937 INITMSG(MYDEBUG,"ePYRA13"<<std::endl);
1939 if(TPyra13a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1942 if(TPyra13b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1948 INITMSG(MYDEBUG,"ePENTA15"<<std::endl);
1950 if(TPenta15a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1953 if(TPenta15b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1959 INITMSG(MYDEBUG,"eHEXA20"<<std::endl);
1961 if(THexa20a().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1964 if(THexa20b().Eval(theCellInfo,theNodeInfo,theElemNum,aRefSlice,aGaussSlice,theGaussCoord,theMode))
1970 INITMSG(MYDEBUG,"eNONE"<<std::endl);
1978 //---------------------------------------------------------------
1980 GetBaryCenter(const TCellInfo& theCellInfo,
1981 const TNodeInfo& theNodeInfo,
1982 TGaussCoord& theGaussCoord,
1983 const TElemNum& theElemNum,
1984 EModeSwitch theMode)
1986 INITMSG(MYDEBUG,"GetBaryCenter\n");
1987 const PMeshInfo& aMeshInfo = theCellInfo.GetMeshInfo();
1988 TInt aDim = aMeshInfo->GetDim();
1989 static TInt aNbGauss = 1;
1991 bool anIsSubMesh = !theElemNum.empty();
1994 aNbElem = theElemNum.size();
1996 aNbElem = theCellInfo.GetNbElem();
1998 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2000 TInt aConnDim = theCellInfo.GetConnDim();
2004 "; aNbGauss = "<<aNbGauss<<
2005 "; aNbElem = "<<aNbElem<<
2006 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2009 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2010 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2011 TCConnSlice aConnSlice = theCellInfo.GetConnSlice(aCellId);
2012 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2014 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2015 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2017 for(TInt aConnId = 0; aConnId < aConnDim; aConnId++){
2018 TInt aNodeId = aConnSlice[aConnId] - 1;
2019 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2021 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2022 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2026 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2027 aGaussCoordSlice[aDimId] /= aConnDim;
2033 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2034 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2035 INITMSG(MYVALUEDEBUG,"");
2036 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2037 TCoordSlice& aCoordSlice = aCoordSliceArr[aGaussId];
2038 ADDMSG(MYVALUEDEBUG,"{");
2039 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2040 ADDMSG(MYVALUEDEBUG,aCoordSlice[aDimId]<<" ");
2042 ADDMSG(MYVALUEDEBUG,"} ");
2044 ADDMSG(MYVALUEDEBUG,std::endl);
2052 //---------------------------------------------------------------
2054 GetBaryCenter(const TPolygoneInfo& thePolygoneInfo,
2055 const TNodeInfo& theNodeInfo,
2056 TGaussCoord& theGaussCoord,
2057 const TElemNum& theElemNum,
2058 EModeSwitch theMode)
2060 INITMSG(MYDEBUG,"GetBaryCenter\n");
2061 const PMeshInfo& aMeshInfo = thePolygoneInfo.GetMeshInfo();
2062 TInt aDim = aMeshInfo->GetDim();
2063 static TInt aNbGauss = 1;
2065 bool anIsSubMesh = !theElemNum.empty();
2068 aNbElem = theElemNum.size();
2070 aNbElem = thePolygoneInfo.GetNbElem();
2072 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2076 "; aNbGauss = "<<aNbGauss<<
2077 "; aNbElem = "<<aNbElem<<
2078 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2081 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2082 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2084 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2085 TCConnSlice aConnSlice = thePolygoneInfo.GetConnSlice(aCellId);
2086 TInt aNbConn = thePolygoneInfo.GetNbConn(aCellId);
2087 TInt aNbNodes = thePolygoneInfo.GetNbConn(aCellId);
2089 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2090 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2092 for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2093 TInt aNodeId = aConnSlice[aConnId] - 1;
2094 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2096 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2097 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2101 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2102 aGaussCoordSlice[aDimId] /= aNbNodes;
2111 //---------------------------------------------------------------
2113 GetBaryCenter(const TPolyedreInfo& thePolyedreInfo,
2114 const TNodeInfo& theNodeInfo,
2115 TGaussCoord& theGaussCoord,
2116 const TElemNum& theElemNum,
2117 EModeSwitch theMode)
2119 INITMSG(MYDEBUG,"GetBaryCenter\n");
2120 const PMeshInfo& aMeshInfo = thePolyedreInfo.GetMeshInfo();
2121 TInt aDim = aMeshInfo->GetDim();
2122 static TInt aNbGauss = 1;
2124 bool anIsSubMesh = !theElemNum.empty();
2127 aNbElem = theElemNum.size();
2129 aNbElem = thePolyedreInfo.GetNbElem();
2131 theGaussCoord.Init(aNbElem,aNbGauss,aDim,theMode);
2135 "; aNbGauss = "<<aNbGauss<<
2136 "; aNbElem = "<<aNbElem<<
2137 "; aNbNodes = "<<theNodeInfo.GetNbElem()<<
2140 for(TInt anElemId = 0; anElemId < aNbElem; anElemId++){
2141 TInt aCellId = anIsSubMesh? theElemNum[anElemId]-1: anElemId;
2143 TCoordSliceArr aCoordSliceArr = theGaussCoord.GetCoordSliceArr(anElemId);
2144 TCConnSliceArr aConnSliceArr = thePolyedreInfo.GetConnSliceArr(aCellId);
2145 TInt aNbFaces = aConnSliceArr.size();
2147 TInt aNbNodes = thePolyedreInfo.GetNbNodes(aCellId);
2149 for(TInt aGaussId = 0; aGaussId < aNbGauss; aGaussId++){
2150 TCoordSlice& aGaussCoordSlice = aCoordSliceArr[aGaussId];
2152 for(TInt aFaceId = 0; aFaceId < aNbFaces; aFaceId++){
2153 TCConnSlice aConnSlice = aConnSliceArr[aFaceId];
2154 TInt aNbConn = aConnSlice.size();
2155 for(TInt aConnId = 0; aConnId < aNbConn; aConnId++){
2156 TInt aNodeId = aConnSlice[aConnId] - 1;
2157 TCCoordSlice aNodeCoordSlice = theNodeInfo.GetCoordSlice(aNodeId);
2159 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2160 aGaussCoordSlice[aDimId] += aNodeCoordSlice[aDimId];
2164 for(TInt aDimId = 0; aDimId < aDim; aDimId++){
2165 aGaussCoordSlice[aDimId] /= aNbNodes;