49 StdHexExp::StdHexExp()
57 StdExpansion(Ba.GetNumModes()*Bb.GetNumModes()*Bc.GetNumModes(), 3,
141 ASSERTL1(
false,
"input dir is out of range");
190 "Basis[1] is not a general tensor type");
194 "Basis[2] is not a general tensor type");
197 &&
m_base[2]->Collocation())
202 inarray, 1, outarray, 1);
224 inarray,outarray,wsp,
true,
true,
true);
248 bool doCheckCollDir0,
249 bool doCheckCollDir1,
250 bool doCheckCollDir2)
252 int nquad0 =
m_base[0]->GetNumPoints();
253 int nquad1 =
m_base[1]->GetNumPoints();
254 int nquad2 =
m_base[2]->GetNumPoints();
255 int nmodes0 =
m_base[0]->GetNumModes();
256 int nmodes1 =
m_base[1]->GetNumModes();
257 int nmodes2 =
m_base[2]->GetNumModes();
260 bool colldir0 = doCheckCollDir0?(
m_base[0]->Collocation()):
false;
261 bool colldir1 = doCheckCollDir1?(
m_base[1]->Collocation()):
false;
262 bool colldir2 = doCheckCollDir2?(
m_base[2]->Collocation()):
false;
266 if(colldir0 && colldir1 && colldir2)
273 ASSERTL1(wsp.num_elements()>=nquad0*nmodes2*(nmodes1+nquad1),
274 "Workspace size is not sufficient");
280 Blas::Dgemm(
'T',
'T', nmodes1*nmodes2, nquad0, nmodes0,
281 1.0, &inarray[0], nmodes0,
283 0.0, &wsp[0], nmodes1*nmodes2);
284 Blas::Dgemm(
'T',
'T', nquad0*nmodes2, nquad1, nmodes1,
285 1.0, &wsp[0], nmodes1,
287 0.0, &wsp2[0], nquad0*nmodes2);
288 Blas::Dgemm(
'T',
'T', nquad0*nquad1, nquad2, nmodes2,
289 1.0, &wsp2[0], nmodes2,
291 0.0, &outarray[0], nquad0*nquad1);
311 if( (
m_base[0]->Collocation())
312 &&(
m_base[1]->Collocation())
313 &&(
m_base[2]->Collocation()) )
370 if(
m_base[0]->Collocation() &&
371 m_base[1]->Collocation() &&
392 Blas::Dgemv(
'N',
m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
393 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
402 bool multiplybyweights)
404 int nquad0 =
m_base[0]->GetNumPoints();
405 int nquad1 =
m_base[1]->GetNumPoints();
406 int nquad2 =
m_base[2]->GetNumPoints();
407 int order0 =
m_base[0]->GetNumModes();
408 int order1 =
m_base[1]->GetNumModes();
411 order0*order1*nquad2);
413 if(multiplybyweights)
421 tmp,outarray,wsp,
true,
true,
true);
428 inarray,outarray,wsp,
true,
true,
true);
443 bool doCheckCollDir0,
444 bool doCheckCollDir1,
445 bool doCheckCollDir2)
447 int nquad0 =
m_base[0]->GetNumPoints();
448 int nquad1 =
m_base[1]->GetNumPoints();
449 int nquad2 =
m_base[2]->GetNumPoints();
450 int nmodes0 =
m_base[0]->GetNumModes();
451 int nmodes1 =
m_base[1]->GetNumModes();
452 int nmodes2 =
m_base[2]->GetNumModes();
454 bool colldir0 = doCheckCollDir0?(
m_base[0]->Collocation()):
false;
455 bool colldir1 = doCheckCollDir1?(
m_base[1]->Collocation()):
false;
456 bool colldir2 = doCheckCollDir2?(
m_base[2]->Collocation()):
false;
458 if(colldir0 && colldir1 && colldir2)
464 ASSERTL1(wsp.num_elements() >= nmodes0*nquad2*(nquad1+nmodes1),
465 "Insufficient workspace size");
474 for(
int n = 0; n < nmodes0; ++n)
477 tmp0.get()+nquad1*nquad2*n,1);
482 Blas::Dgemm(
'T',
'N', nquad1*nquad2, nmodes0, nquad0,
483 1.0, inarray.get(), nquad0,
485 0.0, tmp0.get(), nquad1*nquad2);
491 for(
int n = 0; n < nmodes1; ++n)
494 tmp1.get()+nquad2*nmodes0*n,1);
499 Blas::Dgemm(
'T',
'N', nquad2*nmodes0, nmodes1, nquad1,
500 1.0, tmp0.get(), nquad1,
502 0.0, tmp1.get(), nquad2*nmodes0);
508 for(
int n = 0; n < nmodes2; ++n)
511 outarray.get()+nmodes0*nmodes1*n,1);
516 Blas::Dgemm(
'T',
'N', nmodes0*nmodes1, nmodes2, nquad2,
517 1.0, tmp1.get(), nquad2,
519 0.0, outarray.get(), nmodes0*nmodes1);
536 ASSERTL0((dir==0)||(dir==1)||(dir==2),
"input dir is out of range");
557 Blas::Dgemv(
'N',
m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
558 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
566 ASSERTL0((dir==0)||(dir==1)||(dir==2),
"input dir is out of range");
568 int nquad1 =
m_base[1]->GetNumPoints();
569 int nquad2 =
m_base[2]->GetNumPoints();
570 int order0 =
m_base[0]->GetNumModes();
571 int order1 =
m_base[1]->GetNumModes();
575 if (outarray.num_elements() < inarray.num_elements())
628 int nquad0 =
m_base[0]->GetNumPoints();
629 int nquad1 =
m_base[1]->GetNumPoints();
630 int nquad2 =
m_base[2]->GetNumPoints();
636 int btmp0 =
m_base[0]->GetNumModes();
637 int btmp1 =
m_base[1]->GetNumModes();
638 int mode2 = mode/(btmp0*btmp1);
639 int mode1 = (mode-mode2*btmp0*btmp1)/btmp0;
640 int mode0 = (mode-mode2*btmp0*btmp1)%btmp0;
642 ASSERTL2(mode2 == (
int)floor((1.0*mode)/(btmp0*btmp1)),
643 "Integer Truncation not Equiv to Floor");
644 ASSERTL2(mode1 == (
int)floor((1.0*mode-mode2*btmp0*btmp1)
646 "Integer Truncation not Equiv to Floor");
648 "calling argument mode is larger than total expansion "
651 for(i = 0; i < nquad1*nquad2; ++i)
654 &outarray[0]+i*nquad0, 1);
657 for(j = 0; j < nquad2; ++j)
659 for(i = 0; i < nquad0; ++i)
662 &outarray[0]+i+j*nquad0*nquad1, nquad0,
663 &outarray[0]+i+j*nquad0*nquad1, nquad0);
667 for(i = 0; i < nquad2; i++)
669 Blas::Dscal(nquad0*nquad1,base2[mode2*nquad2+i],
670 &outarray[0]+i*nquad0*nquad1,1);
703 "BasisType is not a boundary interior form");
706 "BasisType is not a boundary interior form");
709 "BasisType is not a boundary interior form");
711 int nmodes0 =
m_base[0]->GetNumModes();
712 int nmodes1 =
m_base[1]->GetNumModes();
713 int nmodes2 =
m_base[2]->GetNumModes();
715 return ( 2*( nmodes0*nmodes1 + nmodes0*nmodes2
717 - 4*( nmodes0 + nmodes1 + nmodes2 ) + 8 );
724 "BasisType is not a boundary interior form");
727 "BasisType is not a boundary interior form");
730 "BasisType is not a boundary interior form");
732 int nmodes0 =
m_base[0]->GetNumModes();
733 int nmodes1 =
m_base[1]->GetNumModes();
734 int nmodes2 =
m_base[2]->GetNumModes();
736 return 2*( nmodes0*nmodes1 + nmodes0*nmodes2
742 ASSERTL2((i >= 0)&&(i <= 11),
"edge id is out of range");
744 if((i == 0)||(i == 2)||(i == 8)||(i == 10))
748 else if((i == 1)||(i == 3)||(i == 9)||(i == 11))
766 ASSERTL2((i >= 0) && (i <= 5),
"face id is out of range");
767 if((i == 0) || (i == 5))
771 else if((i == 1) || (i == 3))
784 ASSERTL2((i >= 0) && (i <= 5),
"face id is out of range");
785 if((i == 0) || (i == 5))
789 else if((i == 1) || (i == 3))
809 ASSERTL2(i >= 0 && i <= 5,
"face id is out of range");
811 if (i == 0 || i == 5)
813 return m_base[0]->GetNumPoints()*
814 m_base[1]->GetNumPoints();
816 else if (i == 1 || i == 3)
818 return m_base[0]->GetNumPoints()*
819 m_base[2]->GetNumPoints();
823 return m_base[1]->GetNumPoints()*
824 m_base[2]->GetNumPoints();
829 const int i,
const int j)
const
831 ASSERTL2(i >= 0 && i <= 5,
"face id is out of range");
832 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
834 if (i == 0 || i == 5)
836 return m_base[j]->GetPointsKey();
838 else if (i == 1 || i == 3)
840 return m_base[2*j]->GetPointsKey();
844 return m_base[j+1]->GetPointsKey();
850 int nmodes = nummodes[modes_offset]*nummodes[modes_offset+1]*nummodes[modes_offset+2];
858 const int i,
const int k)
const
860 ASSERTL2(i >= 0 && i <= 5,
"face id is out of range");
861 ASSERTL2(k >= 0 && k <= 1,
"basis key id is out of range");
883 m_base[dir]->GetNumModes());
888 ASSERTL2((i >= 0)&&(i <= 11),
"edge id is out of range");
890 if((i == 0)||(i == 2)||(i==8)||(i==10))
894 else if((i == 1)||(i == 3)||(i == 9)||(i == 11))
917 for(
int k = 0; k < Qz; ++k ) {
918 for(
int j = 0; j < Qy; ++j ) {
919 for(
int i = 0; i < Qx; ++i ) {
920 int s = i + Qx*(j + Qy*k);
943 int nummodesA, nummodesB;
947 "Method only implemented if BasisType is indentical in "
951 "Method only implemented for Modified_A or GLL_Lagrange BasisType");
953 const int nummodes0 =
m_base[0]->GetNumModes();
954 const int nummodes1 =
m_base[1]->GetNumModes();
955 const int nummodes2 =
m_base[2]->GetNumModes();
961 nummodesA = nummodes0;
962 nummodesB = nummodes1;
966 nummodesA = nummodes0;
967 nummodesB = nummodes2;
971 nummodesA = nummodes1;
972 nummodesB = nummodes2;
976 bool CheckForZeroedModes =
false;
984 if((P != nummodesA)||(Q != nummodesB))
986 CheckForZeroedModes =
true;
990 int nFaceCoeffs = P*Q;
992 if(maparray.num_elements() != nFaceCoeffs)
997 if(signarray.num_elements() != nFaceCoeffs)
1003 fill( signarray.get() , signarray.get()+nFaceCoeffs, 1 );
1008 for(i = 0; i < Q; i++)
1010 for(j = 0; j < P; j++)
1012 if( faceOrient < 9 )
1014 arrayindx[i*P+j] = i*P+j;
1018 arrayindx[i*P+j] = j*Q+i;
1033 offset = nummodes0*nummodes1;
1037 offset = (nummodes2-1)*nummodes0*nummodes1;
1054 offset = nummodes0*(nummodes1-1);
1055 jump1 = nummodes0*nummodes1;
1060 jump1 = nummodes0*nummodes1;
1071 offset = nummodes0-1;
1072 jump1 = nummodes0*nummodes1;
1079 jump1 = nummodes0*nummodes1;
1084 ASSERTL0(
false,
"fid must be between 0 and 5");
1087 for(i = 0; i < Q; i++)
1089 for(j = 0; j < P; j++)
1091 maparray[ arrayindx[i*P+j] ]
1092 = i*jump1 + j*jump2 + offset;
1097 if(CheckForZeroedModes)
1103 for(i = 0; i < nummodesB; i++)
1105 for(j = nummodesA; j < P; j++)
1107 signarray[arrayindx[i*P+j]] = 0.0;
1108 maparray[arrayindx[i*P+j]] = maparray[0];
1112 for(i = nummodesB; i < Q; i++)
1114 for(j = 0; j < P; j++)
1116 signarray[arrayindx[i*P+j]] = 0.0;
1117 maparray[arrayindx[i*P+j]] = maparray[0];
1123 ASSERTL0(
false,
"Different trace space face dimention and element face dimention not possible for GLL-Lagrange bases");
1127 if( (faceOrient==6) || (faceOrient==8) ||
1128 (faceOrient==11) || (faceOrient==12) )
1134 for(i = 3; i < Q; i+=2)
1136 for(j = 0; j < P; j++)
1138 signarray[ arrayindx[i*P+j] ] *= -1;
1142 for(i = 0; i < P; i++)
1144 swap( maparray[i] , maparray[i+P] );
1145 swap( signarray[i] , signarray[i+P] );
1151 for(i = 0; i < P; i++)
1153 for(j = 0; j < Q/2; j++)
1155 swap( maparray[i + j*P],
1158 swap( signarray[i + j*P],
1169 for(i = 0; i < Q; i++)
1171 for(j = 3; j < P; j+=2)
1173 signarray[ arrayindx[i*P+j] ] *= -1;
1177 for(i = 0; i < Q; i++)
1179 swap( maparray[i] , maparray[i+Q] );
1180 swap( signarray[i] , signarray[i+Q] );
1186 for(i = 0; i < P; i++)
1188 for(j = 0; j < Q/2; j++)
1190 swap( maparray[i*Q + j],
1191 maparray[i*Q + Q -1 -j]);
1192 swap( signarray[i*Q + j],
1193 signarray[i*Q + Q -1 -j]);
1200 if( (faceOrient==7) || (faceOrient==8) ||
1201 (faceOrient==10) || (faceOrient==12) )
1207 for(i = 0; i < Q; i++)
1209 for(j = 3; j < P; j+=2)
1211 signarray[ arrayindx[i*P+j] ] *= -1;
1215 for(i = 0; i < Q; i++)
1217 swap( maparray[i*P],
1219 swap( signarray[i*P],
1225 for(i = 0; i < Q; i++)
1227 for(j = 0; j < P/2; j++)
1229 swap( maparray[i*P + j],
1230 maparray[i*P + P -1 -j]);
1231 swap( signarray[i*P + j],
1232 signarray[i*P + P -1 -j]);
1244 for(i = 3; i < Q; i+=2)
1246 for(j = 0; j < P; j++)
1248 signarray[ arrayindx[i*P+j] ] *= -1;
1252 for(i = 0; i < P; i++)
1254 swap( maparray[i*Q],
1256 swap( signarray[i*Q],
1262 for(i = 0; i < Q; i++)
1264 for(j = 0; j < P/2; j++)
1266 swap( maparray[i + j*Q] ,
1267 maparray[i+P*Q - Q -j*Q] );
1268 swap( signarray[i + j*Q] ,
1269 signarray[i+P*Q - Q -j*Q] );
1290 "BasisType is not a boundary interior form");
1293 "BasisType is not a boundary interior form");
1296 "BasisType is not a boundary interior form");
1298 ASSERTL1((localVertexId>=0)&&(localVertexId<8),
1299 "local vertex id must be between 0 and 7");
1306 int nummodes [3] = {
m_base[0]->GetNumModes(),
1307 m_base[1]->GetNumModes(),
1308 m_base[2]->GetNumModes()};
1310 if(useCoeffPacking ==
true)
1312 if(localVertexId > 3)
1324 switch(localVertexId % 4)
1371 if( (localVertexId % 4) % 3 > 0 )
1384 if( localVertexId % 4 > 1 )
1397 if( localVertexId > 3)
1410 return r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1427 "BasisType is not a boundary interior form");
1430 "BasisType is not a boundary interior form");
1433 "BasisType is not a boundary interior form");
1436 "local edge id must be between 0 and 11");
1440 if(maparray.num_elements()!=nEdgeIntCoeffs)
1445 if(signarray.num_elements() != nEdgeIntCoeffs)
1451 fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1454 int nummodes [3] = {
m_base[0]->GetNumModes(),
1455 m_base[1]->GetNumModes(),
1456 m_base[2]->GetNumModes()};
1462 bool reverseOrdering =
false;
1463 bool signChange =
false;
1465 int IdxRange [3][2] = {{0,0},{0,0},{0,0}};
1485 IdxRange[2][0] = nummodes[2] - 1;
1486 IdxRange[2][1] = nummodes[2];
1503 IdxRange[2][1] = nummodes[2] - 1;
1507 reverseOrdering =
true;
1513 IdxRange[2][1] = nummodes[2];
1542 IdxRange[1][0] = nummodes[1] - 1;
1543 IdxRange[1][1] = nummodes[1];
1558 IdxRange[1][1] = nummodes[1] - 1;
1562 reverseOrdering =
true;
1568 IdxRange[1][1] = nummodes[1];
1583 IdxRange[1][1] = nummodes[1] - 1;
1587 reverseOrdering =
true;
1593 IdxRange[1][1] = nummodes[1];
1622 IdxRange[0][0] = nummodes[0] - 1;
1623 IdxRange[0][1] = nummodes[0];
1638 IdxRange[0][1] = nummodes[0] - 1;
1642 reverseOrdering =
true;
1648 IdxRange[0][1] = nummodes[0];
1663 IdxRange[0][1] = nummodes[0] - 1;
1667 reverseOrdering =
true;
1673 IdxRange[0][1] = nummodes[0];
1687 for(r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1689 for(q = IdxRange[1][0]; q < IdxRange[1][1]; q++)
1691 for(p = IdxRange[0][0]; p < IdxRange[0][1]; p++)
1694 = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
1699 if( reverseOrdering )
1701 reverse( maparray.get() , maparray.get()+nEdgeIntCoeffs );
1706 for(p = 1; p < nEdgeIntCoeffs; p+=2)
1725 "BasisType is not a boundary interior form");
1728 "BasisType is not a boundary interior form");
1731 "BasisType is not a boundary interior form");
1734 "local face id must be between 0 and 5");
1738 if(maparray.num_elements()!=nFaceIntCoeffs)
1743 if(signarray.num_elements() != nFaceIntCoeffs)
1749 fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1752 int nummodes [3] = {
m_base[0]->GetNumModes(),
1753 m_base[1]->GetNumModes(),
1754 m_base[2]->GetNumModes()};
1770 nummodesA = nummodes[0];
1771 nummodesB = nummodes[1];
1777 nummodesA = nummodes[0];
1778 nummodesB = nummodes[2];
1784 nummodesA = nummodes[1];
1785 nummodesB = nummodes[2];
1794 for(i = 0; i < (nummodesB-2); i++)
1796 for(j = 0; j < (nummodesA-2); j++)
1798 if( faceOrient < 9 )
1800 arrayindx[i*(nummodesA-2)+j] = i*(nummodesA-2)+j;
1804 arrayindx[i*(nummodesA-2)+j] = j*(nummodesB-2)+i;
1809 int IdxRange [3][2];
1831 IdxRange[2][0] = nummodes[2] - 1;
1832 IdxRange[2][1] = nummodes[2];
1848 if( (((
int) faceOrient)-5) % 2 )
1850 IdxRange[2][0] = nummodes[2] - 2;
1858 IdxRange[2][1] = nummodes[2] - 1;
1865 IdxRange[2][1] = nummodes[2];
1868 if( (((
int) faceOrient)-5) % 2 )
1870 for(i = 3; i < nummodes[2]; i+=2)
1894 IdxRange[1][0] = nummodes[1] - 1;
1895 IdxRange[1][1] = nummodes[1];
1911 if( (((
int) faceOrient)-5) % 2 )
1913 IdxRange[1][0] = nummodes[1] - 2;
1921 IdxRange[1][1] = nummodes[1] - 1;
1928 IdxRange[1][1] = nummodes[1];
1931 if( (((
int) faceOrient)-5) % 2 )
1933 for(i = 3; i < nummodes[1]; i+=2)
1945 if( (((
int) faceOrient)-5) % 4 > 1 )
1947 IdxRange[1][0] = nummodes[1] - 2;
1955 IdxRange[1][1] = nummodes[1] - 1;
1962 IdxRange[1][1] = nummodes[1];
1965 if( (((
int) faceOrient)-5) % 4 > 1 )
1967 for(i = 3; i < nummodes[1]; i+=2)
1989 IdxRange[0][0] = nummodes[0] - 1;
1990 IdxRange[0][1] = nummodes[0];
2005 if( (((
int) faceOrient)-5) % 4 > 1 )
2007 IdxRange[0][0] = nummodes[0] - 2;
2015 IdxRange[0][1] = nummodes[0] - 1;
2022 IdxRange[0][1] = nummodes[0];
2025 if( (((
int) faceOrient)-5) % 4 > 1 )
2027 for(i = 3; i < nummodes[0]; i+=2)
2039 for(r = IdxRange[2][0]; r != IdxRange[2][1]; r+=Incr[2])
2041 for(q = IdxRange[1][0]; q != IdxRange[1][1]; q+=Incr[1])
2043 for(p = IdxRange[0][0]; p != IdxRange[0][1]; p+=Incr[0])
2045 maparray [ arrayindx[cnt ] ]
2046 = r*nummodes[0]*nummodes[1] + q*nummodes[0] + p;
2047 signarray[ arrayindx[cnt++] ]
2048 = sign0[p] * sign1[q] * sign2[r];
2062 "BasisType is not a boundary interior form");
2065 "BasisType is not a boundary interior form");
2068 "BasisType is not a boundary interior form");
2071 int nummodes [3] = {
m_base[0]->GetNumModes(),
2072 m_base[1]->GetNumModes(),
2073 m_base[2]->GetNumModes()};
2077 if(outarray.num_elements() != nIntCoeffs)
2091 for(i = 0; i < 3; i++)
2096 IntIdx[i][1] = nummodes[i];
2101 IntIdx[i][1] = nummodes[i]-1;
2105 for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2107 for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2109 for( p = IntIdx[0][0]; p < IntIdx[0][1]; p++)
2111 outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2126 "BasisType is not a boundary interior form");
2129 "BasisType is not a boundary interior form");
2132 "BasisType is not a boundary interior form");
2135 int nummodes [3] = {
m_base[0]->GetNumModes(),
2136 m_base[1]->GetNumModes(),
2137 m_base[2]->GetNumModes()};
2141 if(outarray.num_elements()!=nBndCoeffs)
2156 for(i = 0; i < 3; i++)
2164 IntIdx[i][1] = nummodes[i];
2168 BndIdx[i][1] = nummodes[i]-1;
2170 IntIdx[i][1] = nummodes[i]-1;
2175 for(i = 0; i < 2; i++)
2178 for( q = 0; q < nummodes[1]; q++)
2180 for( p = 0; p < nummodes[0]; p++)
2182 outarray[cnt++] = r*nummodes[0]*nummodes[1]+q*nummodes[0] + p;
2187 for(r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
2189 for( i = 0; i < 2; i++)
2192 for( p = 0; p < nummodes[0]; p++)
2194 outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2199 for( q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
2201 for( i = 0; i < 2; i++)
2204 outarray[cnt++] = r*nummodes[0]*nummodes[1] +
2210 sort(outarray.get(), outarray.get() + nBndCoeffs);
2278 if(inarray.get() == outarray.get())
2284 m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
2289 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
2298 int nquad0 =
m_base[0]->GetNumPoints();
2299 int nquad1 =
m_base[1]->GetNumPoints();
2300 int nquad2 =
m_base[2]->GetNumPoints();
2301 int nq01 = nquad0*nquad1;
2302 int nq12 = nquad1*nquad2;
2308 for(i = 0; i < nq12; ++i)
2311 w0.get(), 1, outarray.get()+i*nquad0,1);
2314 for(i = 0; i < nq12; ++i)
2316 Vmath::Smul(nquad0, w1[i%nquad1], outarray.get()+i*nquad0, 1,
2317 outarray.get()+i*nquad0, 1);
2320 for(i = 0; i < nquad2; ++i)
2322 Vmath::Smul(nq01, w2[i], outarray.get()+i*nq01, 1,
2323 outarray.get()+i*nq01, 1);
2331 int qa =
m_base[0]->GetNumPoints();
2332 int qb =
m_base[1]->GetNumPoints();
2333 int qc =
m_base[2]->GetNumPoints();
2334 int nmodes_a =
m_base[0]->GetNumModes();
2335 int nmodes_b =
m_base[1]->GetNumModes();
2336 int nmodes_c =
m_base[2]->GetNumModes();
2354 OrthoExp.
FwdTrans(array,orthocoeffs);
2358 int nmodes = max(nmodes_a,nmodes_b);
2359 nmodes = max(nmodes,nmodes_c);
2362 for(j = cutoff; j < nmodes; ++j)
2364 fac[j] = fabs((j-nmodes)/((
NekDouble) (j-cutoff+1.0)));
2368 for(i = 0; i < nmodes_a; ++i)
2370 for(j = 0; j < nmodes_b; ++j)
2372 for(k = 0; k < nmodes_c; ++k)
2374 if((i >= cutoff)||(j >= cutoff)||(k >= cutoff))
2376 orthocoeffs[i*nmodes_a*nmodes_b + j*nmodes_c + k] *= (SvvDiffCoeff*exp( -(fac[i]+fac[j]+fac[k]) ));
2380 orthocoeffs[i*nmodes_a*nmodes_b + j*nmodes_c + k] *= 0.0;
2387 OrthoExp.
BwdTrans(orthocoeffs,array);
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
virtual int v_GetEdgeNcoeffs(const int i) const
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
NekDouble GetConstFactor(const ConstFactorType &factor) const
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
#define ASSERTL0(condition, msg)
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual int v_NumBndryCoeffs() const
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Principle Modified Functions .
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
static Array< OneD, NekDouble > NullNekDouble1DArray
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetNedges() const
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Principle Modified Functions .
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetFaceNcoeffs(const int i) const
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
virtual bool v_IsBoundaryInteriorExpansion()
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
boost::shared_ptr< DNekMat > DNekMatSharedPtr
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
virtual int v_GetNfaces() const
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Principle Orthogonal Functions .
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
virtual int v_GetFaceNumPoints(const int i) const
int NumBndryCoeffs(void) const
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points...
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
The base class for all shapes.
virtual int v_GetTotalEdgeIntNcoeffs() const
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
int GetFaceIntNcoeffs(const int i) const
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Principle Modified Functions .
Class representing a hexehedral element in reference space.
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Principle Orthogonal Functions .
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Principle Orthogonal Functions .
Defines a specification for a set of points.
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int nummodesA=-1, int nummodesB=-1)
virtual int v_GetTotalFaceIntNcoeffs() const
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
virtual LibUtilities::ShapeType v_DetShapeType() const
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Differentiation Methods.
virtual int v_GetFaceIntNcoeffs(const int i) const
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Array< OneD, LibUtilities::BasisSharedPtr > m_base
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual int v_GetNverts() const
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multbyweights=true)
Describes the specification for a Basis.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...
virtual int v_NumDGBndryCoeffs() const