47 StdTetExp::StdTetExp()
67 "order in 'a' direction is higher than order "
70 "order in 'a' direction is higher than order "
73 "order in 'b' direction is higher than order "
114 int Q0 =
m_base[0]->GetNumPoints();
115 int Q1 =
m_base[1]->GetNumPoints();
116 int Q2 =
m_base[2]->GetNumPoints();
124 bool Do_2 = (out_dxi2.num_elements() > 0)?
true:
false;
125 bool Do_1 = (out_dxi1.num_elements() > 0)?
true:
false;
142 eta_0 =
m_base[0]->GetZ();
143 eta_1 =
m_base[1]->GetZ();
144 eta_2 =
m_base[2]->GetZ();
150 for(
int k=0; k< Q2; ++k)
152 for(
int j=0; j<Q1; ++j,dEta0+=Q0)
154 Vmath::Smul(Q0,2.0/(1.0-eta_1[j]),dEta0,1,dEta0,1);
156 fac = 1.0/(1.0-eta_2[k]);
157 Vmath::Smul(Q0*Q1,fac,&out_dEta0[0]+k*Q0*Q1,1,&out_dEta0[0]+k*Q0*Q1,1);
160 if (out_dxi0.num_elements() > 0)
173 for(
int k = 0; k < Q1*Q2; ++k)
175 Vmath::Vmul(Q0,&Fac0[0],1,&out_dEta0[0]+k*Q0,1,&out_dEta0[0]+k*Q0,1);
178 for(
int k = 0; k < Q2; ++k)
180 Vmath::Smul(Q0*Q1,2.0/(1.0-eta_2[k]),&out_dEta1[0]+k*Q0*Q1,1,
181 &out_dEta1[0]+k*Q0*Q1,1);
188 Vmath::Vadd(Qtot,out_dEta0,1,out_dEta1,1,out_dxi1,1);
196 for(
int k=0; k< Q2; ++k)
198 for(
int j=0; j<Q1; ++j,dEta1+=Q0)
200 Vmath::Smul(Q0,(1.0+eta_1[j])/2.0,dEta1,1,dEta1,1);
207 Vmath::Vadd(Qtot,out_dEta0,1,out_dEta1,1,out_dxi2,1);
208 Vmath::Vadd(Qtot,out_dEta2,1,out_dxi2 ,1,out_dxi2,1);
247 ASSERTL1(
false,
"input dir is out of range");
301 "Basis[1] is not a general tensor type");
305 "Basis[2] is not a general tensor type");
308 &&
m_base[2]->Collocation())
313 inarray, 1, outarray, 1);
329 int nquad1 =
m_base[1]->GetNumPoints();
330 int nquad2 =
m_base[2]->GetNumPoints();
331 int order0 =
m_base[0]->GetNumModes();
332 int order1 =
m_base[1]->GetNumModes();
335 nquad2*nquad1*order0);
340 inarray,outarray,wsp,
365 bool doCheckCollDir0,
366 bool doCheckCollDir1,
367 bool doCheckCollDir2)
369 int nquad0 =
m_base[0]->GetNumPoints();
370 int nquad1 =
m_base[1]->GetNumPoints();
371 int nquad2 =
m_base[2]->GetNumPoints();
373 int order0 =
m_base[0]->GetNumModes();
374 int order1 =
m_base[1]->GetNumModes();
375 int order2 =
m_base[2]->GetNumModes();
380 int i, j, mode, mode1, cnt;
383 mode = mode1 = cnt = 0;
384 for(i = 0; i < order0; ++i)
386 for(j = 0; j < order1-i; ++j, ++cnt)
388 Blas::Dgemv(
'N', nquad2, order2-i-j,
389 1.0, base2.get()+mode*nquad2, nquad2,
390 inarray.get()+mode1, 1,
391 0.0, tmp.get()+cnt*nquad2, 1);
396 for(j = order1-i; j < order2-i; ++j)
408 Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
412 Blas::Daxpy(nquad2,inarray[1],base2.get()+nquad2,1,
413 &tmp[0]+order1*nquad2,1);
418 for(i = 0; i < order0; ++i)
420 Blas::Dgemm(
'N',
'T', nquad1, nquad2, order1-i,
421 1.0, base1.get()+mode*nquad1, nquad1,
422 tmp.get()+mode*nquad2, nquad2,
423 0.0, tmp1.get()+i*nquad1*nquad2, nquad1);
434 for(i = 0; i < nquad2; ++i)
436 Blas::Daxpy(nquad1,tmp[nquad2+i], base1.get()+nquad1,1,
437 &tmp1[nquad1*nquad2]+i*nquad1,1);
442 Blas::Dgemm(
'N',
'T', nquad0, nquad1*nquad2, order0,
443 1.0, base0.get(), nquad0,
444 tmp1.get(), nquad1*nquad2,
445 0.0, outarray.get(), nquad0);
511 "Basis[1] is not a general tensor type");
515 "Basis[2] is not a general tensor type");
517 if(
m_base[0]->Collocation() &&
m_base[1]->Collocation())
536 Blas::Dgemv(
'N',
m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
537 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
550 bool multiplybyweights)
552 int nquad0 =
m_base[0]->GetNumPoints();
553 int nquad1 =
m_base[1]->GetNumPoints();
554 int nquad2 =
m_base[2]->GetNumPoints();
555 int order0 =
m_base[0]->GetNumModes();
556 int order1 =
m_base[1]->GetNumModes();
559 nquad2*order0*(2*order1-order0+1)/2);
561 if(multiplybyweights)
570 tmp, outarray, wsp,
true,
true,
true);
578 inarray, outarray, wsp,
true,
true,
true);
590 bool doCheckCollDir0,
591 bool doCheckCollDir1,
592 bool doCheckCollDir2)
594 int nquad0 =
m_base[0]->GetNumPoints();
595 int nquad1 =
m_base[1]->GetNumPoints();
596 int nquad2 =
m_base[2]->GetNumPoints();
598 int order0 =
m_base[0]->GetNumModes();
599 int order1 =
m_base[1]->GetNumModes();
600 int order2 =
m_base[2]->GetNumModes();
605 int i,j, mode,mode1, cnt;
608 Blas::Dgemm(
'T',
'N', nquad1*nquad2, order0, nquad0,
609 1.0, inarray.get(), nquad0,
611 0.0, tmp1.get(), nquad1*nquad2);
614 for(mode=i=0; i < order0; ++i)
616 Blas::Dgemm(
'T',
'N', nquad2, order1-i, nquad1,
617 1.0, tmp1.get()+i*nquad1*nquad2, nquad1,
618 base1.get()+mode*nquad1, nquad1,
619 0.0, tmp2.get()+mode*nquad2, nquad2);
628 Blas::Dgemv(
'T', nquad1, nquad2,
629 1.0, tmp1.get()+nquad1*nquad2, nquad1,
630 base1.get()+nquad1, 1,
631 1.0, tmp2.get()+nquad2, 1);
635 mode = mode1 = cnt = 0;
636 for(i = 0; i < order0; ++i)
638 for(j = 0; j < order1-i; ++j, ++cnt)
640 Blas::Dgemv(
'T', nquad2, order2-i-j,
641 1.0, base2.get()+mode*nquad2, nquad2,
642 tmp2.get()+cnt*nquad2, 1,
643 0.0, outarray.get()+mode1, 1);
648 for(j = order1-i; j < order2-i; ++j)
659 outarray[1] +=
Blas::Ddot(nquad2,base2.get()+nquad2,1,
663 outarray[1] +=
Blas::Ddot(nquad2,base2.get()+nquad2,1,
664 &tmp2[nquad2*order1],1);
683 ASSERTL0((dir==0)||(dir==1)||(dir==2),
"input dir is out of range");
704 Blas::Dgemv(
'N',
m_ncoeffs,nq,1.0,iprodmat->GetPtr().get(),
705 m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
721 int nquad0 =
m_base[0]->GetNumPoints();
722 int nquad1 =
m_base[1]->GetNumPoints();
723 int nquad2 =
m_base[2]->GetNumPoints();
724 int nqtot = nquad0*nquad1*nquad2;
725 int nmodes0 =
m_base[0]->GetNumModes();
726 int nmodes1 =
m_base[1]->GetNumModes();
727 int wspsize = nquad0 + nquad1 + nquad2 + max(nqtot,
m_ncoeffs)
728 + nquad1*nquad2*nmodes0 + nquad2*nmodes0*(2*nmodes1-nmodes0+1)/2;
741 for(i = 0; i < nquad0; ++i)
743 gfac0[i] = 0.5*(1+z0[i]);
747 for(i = 0; i < nquad1; ++i)
749 gfac1[i] = 2.0/(1-z1[i]);
753 for(i = 0; i < nquad2; ++i)
755 gfac2[i] = 2.0/(1-z2[i]);
759 for(i = 0; i < nquad1*nquad2; ++i)
761 Vmath::Smul(nquad0,gfac1[i%nquad1],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
763 for(i = 0; i < nquad2; ++i)
765 Vmath::Smul(nquad0*nquad1,gfac2[i],&tmp0[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
785 for(i = 0; i < nquad1*nquad2; ++i)
787 Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
796 for(i = 0; i < nquad2; ++i)
798 Vmath::Smul(nquad0*nquad1,gfac2[i],&inarray[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
813 for(i = 0; i < nquad1; ++i)
815 gfac1[i] = (1+z1[i])/2;
818 for(i = 0; i < nquad1*nquad2; ++i)
820 Vmath::Vmul(nquad0,&gfac0[0],1,&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
828 for(i = 0; i < nquad2; ++i)
830 Vmath::Smul(nquad0*nquad1,gfac2[i],&inarray[0]+i*nquad0*nquad1,1,&tmp0[0]+i*nquad0*nquad1,1);
832 for(i = 0; i < nquad1*nquad2; ++i)
834 Vmath::Smul(nquad0,gfac1[i%nquad1],&tmp0[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
856 ASSERTL1(
false,
"input dir is out of range");
893 eta[0] = 2.0*(1.0+xi[0])/(-xi[1]-xi[2]) - 1.0;
895 eta[1] = 2.0*(1.0+xi[1])/(1.0-xi[2]) - 1.0;
938 "BasisType is not a boundary interior form");
941 "BasisType is not a boundary interior form");
944 "BasisType is not a boundary interior form");
946 int P =
m_base[0]->GetNumModes();
947 int Q =
m_base[1]->GetNumModes();
948 int R =
m_base[2]->GetNumModes();
958 "BasisType is not a boundary interior form");
961 "BasisType is not a boundary interior form");
964 "BasisType is not a boundary interior form");
966 int P =
m_base[0]->GetNumModes()-1;
967 int Q =
m_base[1]->GetNumModes()-1;
968 int R =
m_base[2]->GetNumModes()-1;
971 return (Q+1) + P*(1 + 2*Q - P)/2
972 + (R+1) + P*(1 + 2*R - P)/2
973 + 2*(R+1) + Q*(1 + 2*R - Q);
978 ASSERTL2((i >= 0) && (i <= 5),
"edge id is out of range");
979 int P =
m_base[0]->GetNumModes();
980 int Q =
m_base[1]->GetNumModes();
981 int R =
m_base[2]->GetNumModes();
987 else if (i == 1 || i == 2)
999 int P =
m_base[0]->GetNumModes()-2;
1000 int Q =
m_base[1]->GetNumModes()-2;
1001 int R =
m_base[2]->GetNumModes()-2;
1008 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1009 int nFaceCoeffs = 0;
1010 int nummodesA, nummodesB, P, Q;
1016 else if ((i == 1) || (i == 2))
1028 nFaceCoeffs = Q+1 + (P*(1 + 2*Q - P))/2;
1034 ASSERTL2((i >= 0) && (i <= 3),
"face id is out of range");
1035 int Pi =
m_base[0]->GetNumModes() - 2;
1036 int Qi =
m_base[1]->GetNumModes() - 2;
1037 int Ri =
m_base[2]->GetNumModes() - 2;
1041 return Pi * (2*Qi - Pi - 1) / 2;
1045 return Pi * (2*Ri - Pi - 1) / 2;
1049 return Qi * (2*Ri - Qi - 1) / 2;
1055 int Pi =
m_base[0]->GetNumModes() - 2;
1056 int Qi =
m_base[1]->GetNumModes() - 2;
1057 int Ri =
m_base[2]->GetNumModes() - 2;
1059 return Pi * (2*Qi - Pi - 1) / 2 +
1060 Pi * (2*Ri - Pi - 1) / 2 +
1061 Qi * (2*Ri - Qi - 1);
1066 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1070 return m_base[0]->GetNumPoints()*
1071 m_base[1]->GetNumPoints();
1075 return m_base[0]->GetNumPoints()*
1076 m_base[2]->GetNumPoints();
1080 return m_base[1]->GetNumPoints()*
1081 m_base[2]->GetNumPoints();
1086 const int i,
const int j)
const
1088 ASSERTL2(i >= 0 && i <= 3,
"face id is out of range");
1089 ASSERTL2(j == 0 || j == 1,
"face direction is out of range");
1093 return m_base[j]->GetPointsKey();
1097 return m_base[2*j]->GetPointsKey();
1101 return m_base[j+1]->GetPointsKey();
1106 const std::vector<unsigned int>& nummodes,
1110 nummodes[modes_offset],
1111 nummodes[modes_offset+1],
1112 nummodes[modes_offset+2]);
1119 const int i,
const int k)
const
1121 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
1122 ASSERTL2(k == 0 || k == 1,
"face direction out of range");
1142 m_base[dir]->GetNumModes());
1150 ASSERTL2(i >= 0 && i <= 5,
"edge id is out of range");
1156 else if (i == 1 || i == 2)
1180 for(
int k = 0; k < Qz; ++k ) {
1181 for(
int j = 0; j < Qy; ++j ) {
1182 for(
int i = 0; i < Qx; ++i ) {
1183 int s = i + Qx*(j + Qy*k);
1184 xi_x[s] = (eta_x[i] + 1.0) * (1.0 - eta_y[j]) * (1.0 - eta_z[k]) / 4 - 1.0;
1185 xi_y[s] = (eta_y[j] + 1.0) * (1.0 - eta_z[k]) / 2 - 1.0;
1216 int nummodesA,nummodesB, i, j, k, idx;
1219 "Method only implemented for Modified_A BasisType (x "
1220 "direction), Modified_B BasisType (y direction), and "
1221 "Modified_C BasisType(z direction)");
1223 int nFaceCoeffs = 0;
1228 nummodesA =
m_base[0]->GetNumModes();
1229 nummodesB =
m_base[1]->GetNumModes();
1232 nummodesA =
m_base[0]->GetNumModes();
1233 nummodesB =
m_base[2]->GetNumModes();
1237 nummodesA =
m_base[1]->GetNumModes();
1238 nummodesB =
m_base[2]->GetNumModes();
1242 bool CheckForZeroedModes =
false;
1250 CheckForZeroedModes =
true;
1253 nFaceCoeffs = P*(2*Q-P+1)/2;
1256 if(maparray.num_elements() != nFaceCoeffs)
1261 if(signarray.num_elements() != nFaceCoeffs)
1267 fill(signarray.get(),signarray.get()+nFaceCoeffs, 1 );
1274 for (i = 0; i < P; ++i)
1276 for (j = 0; j < Q-i; ++j)
1278 if ((
int)faceOrient == 7 && i > 1)
1280 signarray[idx] = (i%2 ? -1 : 1);
1282 maparray[idx++] =
GetMode(i,j,0);
1288 for (i = 0; i < P; ++i)
1290 for (k = 0; k < Q-i; ++k)
1292 if ((
int)faceOrient == 7 && i > 1)
1294 signarray[idx] = (i%2 ? -1: 1);
1296 maparray[idx++] =
GetMode(i,0,k);
1302 for (j = 0; j < P-1; ++j)
1304 for (k = 0; k < Q-1-j; ++k)
1306 if ((
int)faceOrient == 7 && j > 1)
1308 signarray[idx] = ((j+1)%2 ? -1: 1);
1310 maparray[idx++] =
GetMode(1,j,k);
1312 if (j == 0 && k == 0)
1314 maparray[idx++] =
GetMode(0,0,1);
1316 if (j == 0 && k == Q-2)
1318 for (
int r = 0; r < Q-1; ++r)
1320 maparray[idx++] =
GetMode(0,1,r);
1328 for (j = 0; j < P; ++j)
1330 for (k = 0; k < Q-j; ++k)
1332 if ((
int)faceOrient == 7 && j > 1)
1334 signarray[idx] = (j%2 ? -1: 1);
1336 maparray[idx++] =
GetMode(0,j,k);
1341 ASSERTL0(
false,
"Element map not available.");
1344 if ((
int)faceOrient == 7)
1346 swap(maparray[0], maparray[Q]);
1348 for (i = 1; i < Q-1; ++i)
1350 swap(maparray[i+1], maparray[Q+i]);
1354 if(CheckForZeroedModes)
1359 for (j = 0; j < nummodesA; ++j)
1362 for (k = nummodesB-j; k < Q-j; ++k)
1364 signarray[idx] = 0.0;
1365 maparray[idx++] = maparray[0];
1369 for (j = nummodesA; j < P; ++j)
1371 for (k = 0; k < Q-j; ++k)
1373 signarray[idx] = 0.0;
1374 maparray[idx++] = maparray[0];
1385 "Mapping not defined for this type of basis");
1388 if(useCoeffPacking ==
true)
1390 switch(localVertexId)
1414 ASSERTL0(
false,
"Vertex ID must be between 0 and 3");
1421 switch(localVertexId)
1445 ASSERTL0(
false,
"Vertex ID must be between 0 and 3");
1465 const int P =
m_base[0]->GetNumModes();
1466 const int Q =
m_base[1]->GetNumModes();
1467 const int R =
m_base[2]->GetNumModes();
1471 if(maparray.num_elements() != nEdgeIntCoeffs)
1477 fill( maparray.get(), maparray.get() + nEdgeIntCoeffs, 0);
1480 if(signarray.num_elements() != nEdgeIntCoeffs)
1486 fill( signarray.get() , signarray.get()+nEdgeIntCoeffs, 1 );
1492 for (i = 0; i < P-2; ++i)
1494 maparray[i] =
GetMode(i+2, 0, 0);
1498 for(i = 1; i < nEdgeIntCoeffs; i+=2)
1505 for (i = 0; i < Q-2; ++i)
1507 maparray[i] =
GetMode(1, i+1, 0);
1511 for(i = 1; i < nEdgeIntCoeffs; i+=2)
1518 for (i = 0; i < Q-2; ++i)
1520 maparray[i] =
GetMode(0, i+2, 0);
1524 for(i = 1; i < nEdgeIntCoeffs; i+=2)
1531 for (i = 0; i < R-2; ++i)
1533 maparray[i] =
GetMode(0, 0, i+2);
1537 for(i = 1; i < nEdgeIntCoeffs; i+=2)
1544 for (i = 0; i < R - 2; ++i)
1546 maparray[i] =
GetMode(1, 0, i+1);
1550 for(i = 1; i < nEdgeIntCoeffs; i+=2)
1557 for (i = 0; i < R-2; ++i)
1559 maparray[i] =
GetMode(0, 1, i+1);
1563 for(i = 1; i < nEdgeIntCoeffs; i+=2)
1570 ASSERTL0(
false,
"Edge not defined.");
1582 const int P =
m_base[0]->GetNumModes();
1583 const int Q =
m_base[1]->GetNumModes();
1584 const int R =
m_base[2]->GetNumModes();
1588 if(maparray.num_elements() != nFaceIntCoeffs)
1593 if(signarray.num_elements() != nFaceIntCoeffs)
1599 fill( signarray.get() , signarray.get()+nFaceIntCoeffs, 1 );
1606 for (i = 2; i < P-1; ++i)
1608 for (j = 1; j < Q-i; ++j)
1610 if ((
int)faceOrient == 7)
1612 signarray[idx] = (i%2 ? -1 : 1);
1614 maparray[idx++] =
GetMode(i,j,0);
1620 for (i = 2; i < P; ++i)
1622 for (k = 1; k < R-i; ++k)
1624 if ((
int)faceOrient == 7)
1626 signarray[idx] = (i%2 ? -1: 1);
1628 maparray[idx++] =
GetMode(i,0,k);
1634 for (j = 1; j < Q-2; ++j)
1636 for (k = 1; k < R-1-j; ++k)
1638 if ((
int)faceOrient == 7)
1640 signarray[idx] = ((j+1)%2 ? -1: 1);
1642 maparray[idx++] =
GetMode(1,j,k);
1648 for (j = 2; j < Q-1; ++j)
1650 for (k = 1; k < R-j; ++k)
1652 if ((
int)faceOrient == 7)
1654 signarray[idx] = (j%2 ? -1: 1);
1656 maparray[idx++] =
GetMode(0,j,k);
1661 ASSERTL0(
false,
"Face interior map not available.");
1673 "BasisType is not a boundary interior form");
1676 "BasisType is not a boundary interior form");
1679 "BasisType is not a boundary interior form");
1681 int P =
m_base[0]->GetNumModes();
1682 int Q =
m_base[1]->GetNumModes();
1683 int R =
m_base[2]->GetNumModes();
1687 if(outarray.num_elements() != nIntCoeffs)
1693 for (
int i = 2; i < P-2; ++i)
1695 for (
int j = 1; j < Q-i-1; ++j)
1697 for (
int k = 1; k < R-i-j; ++k)
1699 outarray[idx++] =
GetMode(i,j,k);
1712 "BasisType is not a boundary interior form");
1715 "BasisType is not a boundary interior form");
1718 "BasisType is not a boundary interior form");
1720 int P =
m_base[0]->GetNumModes();
1721 int Q =
m_base[1]->GetNumModes();
1722 int R =
m_base[2]->GetNumModes();
1729 if (outarray.num_elements() != nBnd)
1734 for (i = 0; i < P; ++i)
1739 for (j = 0; j < Q-i; j++)
1741 for (k = 0; k < R-i-j; ++k)
1743 outarray[idx++] =
GetMode(i,j,k);
1751 for (k = 0; k < R-i; ++k)
1753 outarray[idx++] =
GetMode(i,0,k);
1755 for (j = 1; j < Q-i; ++j)
1757 outarray[idx++] =
GetMode(i,j,0);
1778 int nq0 =
m_base[0]->GetNumPoints();
1779 int nq1 =
m_base[1]->GetNumPoints();
1780 int nq2 =
m_base[2]->GetNumPoints();
1790 nq = max(nq0,max(nq1,nq2));
1804 for(
int i = 0; i < nq; ++i)
1806 for(
int j = 0; j < nq-i; ++j)
1808 for(
int k = 0; k < nq-i-j; ++k,++cnt)
1811 coords[cnt][0] = -1.0 + 2*k/(
NekDouble)(nq-1);
1812 coords[cnt][1] = -1.0 + 2*j/(
NekDouble)(nq-1);
1813 coords[cnt][2] = -1.0 + 2*i/(
NekDouble)(nq-1);
1818 for(
int i = 0; i < neq; ++i)
1822 I[0] =
m_base[0]->GetI(coll);
1823 I[1] =
m_base[1]->GetI(coll+1);
1824 I[2] =
m_base[2]->GetI(coll+2);
1828 for(
int k = 0; k < nq2; ++k)
1830 for (
int j = 0; j < nq1; ++j)
1833 fac = (I[1]->GetPtr())[j]*(I[2]->GetPtr())[k];
1888 const int Q =
m_base[1]->GetNumModes();
1889 const int R =
m_base[2]->GetNumModes();
1891 int i,j,q_hat,k_hat;
1895 for (i = 0; i < I; ++i)
1900 k_hat = max(R-Q-i,0);
1901 cnt += q_hat*(q_hat+1)/2 + k_hat*Q;
1906 for (j = 0; j < J; ++j)
1924 int nquad0 =
m_base[0]->GetNumPoints();
1925 int nquad1 =
m_base[1]->GetNumPoints();
1926 int nquad2 =
m_base[2]->GetNumPoints();
1936 for(i = 0; i < nquad1*nquad2; ++i)
1939 w0.get(),1, &outarray[0]+i*nquad0,1);
1946 for(j = 0; j < nquad2; ++j)
1948 for(i = 0; i < nquad1; ++i)
1950 Blas::Dscal(nquad0,0.5*w1[i], &outarray[0]+i*nquad0+
1957 for(j = 0; j < nquad2; ++j)
1959 for(i = 0; i < nquad1; ++i)
1962 0.5*(1-z1[i])*w1[i],
1963 &outarray[0]+i*nquad0 + j*nquad0*nquad1,
1974 for(i = 0; i < nquad2; ++i)
1976 Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1977 &outarray[0]+i*nquad0*nquad1, 1);
1982 for(i = 0; i < nquad2; ++i)
1984 Blas::Dscal(nquad0*nquad1, 0.25*(1-z2[i])*w2[i],
1985 &outarray[0]+i*nquad0*nquad1, 1);
1989 for(i = 0; i < nquad2; ++i)
1991 Blas::Dscal(nquad0*nquad1,0.25*(1-z2[i])*(1-z2[i])*w2[i],
1992 &outarray[0]+i*nquad0*nquad1,1);
2009 int qa =
m_base[0]->GetNumPoints();
2010 int qb =
m_base[1]->GetNumPoints();
2011 int qc =
m_base[2]->GetNumPoints();
2012 int nmodes_a =
m_base[0]->GetNumModes();
2013 int nmodes_b =
m_base[1]->GetNumModes();
2014 int nmodes_c =
m_base[2]->GetNumModes();
2039 int cutoff_a = (int) (SVVCutOff*nmodes_a);
2040 int cutoff_b = (int) (SVVCutOff*nmodes_b);
2041 int cutoff_c = (int) (SVVCutOff*nmodes_c);
2042 int nmodes = min(min(nmodes_a,nmodes_b),nmodes_c);
2043 NekDouble cutoff = min(min(cutoff_a,cutoff_b),cutoff_c);
2047 OrthoExp.
FwdTrans(array,orthocoeffs);
2050 for(i = 0; i < nmodes_a; ++i)
2052 for(j = 0; j < nmodes_b-i; ++j)
2054 for(k = 0; k < nmodes_c-i-j; ++k)
2056 if(i + j + k >= cutoff)
2058 orthocoeffs[cnt] *= ((SvvDiffCoeff)*exp(-(i+j+k-nmodes)*(i+j+k-nmodes)/((
NekDouble)((i+j+k-cutoff+epsilon)*(i+j+k-cutoff+epsilon)))));
2062 orthocoeffs[cnt] *= 0.0;
2069 OrthoExp.
BwdTrans(orthocoeffs,array);
2078 int nquad0 =
m_base[0]->GetNumPoints();
2079 int nquad1 =
m_base[1]->GetNumPoints();
2080 int nquad2 =
m_base[2]->GetNumPoints();
2081 int nqtot = nquad0 * nquad1 * nquad2;
2082 int nmodes0 =
m_base[0]->GetNumModes();
2083 int nmodes1 =
m_base[1]->GetNumModes();
2084 int nmodes2 =
m_base[2]->GetNumModes();
2085 int numMax = nmodes0;
2113 OrthoTetExp->FwdTrans(phys_tmp, coeff);
2119 for (
int u = 0; u < numMin; ++u)
2121 for (
int i = 0; i < numMin-u; ++i)
2124 tmp2 = coeff_tmp1 + cnt, 1);
2125 cnt += numMax - u - i;
2127 for (
int i = numMin; i < numMax-u; ++i)
2129 cnt += numMax - u - i;
2133 OrthoTetExp->BwdTrans(coeff_tmp1,phys_tmp);
2142 int np0 =
m_base[0]->GetNumPoints();
2143 int np1 =
m_base[1]->GetNumPoints();
2144 int np2 =
m_base[2]->GetNumPoints();
2145 int np = max(np0,max(np1,np2));
2157 for(
int i = 0; i < np-1; ++i)
2159 planep1 += (np-i)*(np-i+1)/2;
2164 for(
int j = 0; j < np-i-1; ++j)
2168 for(
int k = 0; k < np-i-j-2; ++k)
2170 conn[cnt++] = plane + row +k+1;
2171 conn[cnt++] = plane + row +k;
2172 conn[cnt++] = plane + rowp1 +k;
2173 conn[cnt++] = planep1 + row1 +k;
2175 conn[cnt++] = plane + row +k+1;
2176 conn[cnt++] = plane + rowp1 +k+1;
2177 conn[cnt++] = planep1 + row1 +k+1;
2178 conn[cnt++] = planep1 + row1 +k;
2180 conn[cnt++] = plane + rowp1 +k+1;
2181 conn[cnt++] = plane + row +k+1;
2182 conn[cnt++] = plane + rowp1 +k;
2183 conn[cnt++] = planep1 + row1 +k;
2185 conn[cnt++] = planep1 + row1 +k;
2186 conn[cnt++] = planep1 + row1p1+k;
2187 conn[cnt++] = plane + rowp1 +k;
2188 conn[cnt++] = plane + rowp1 +k+1;
2190 conn[cnt++] = planep1 + row1 +k;
2191 conn[cnt++] = planep1 + row1p1+k;
2192 conn[cnt++] = planep1 + row1 +k+1;
2193 conn[cnt++] = plane + rowp1 +k+1;
2197 conn[cnt++] = plane + rowp1 +k+1;
2198 conn[cnt++] = planep1 + row1p1 +k+1;
2199 conn[cnt++] = planep1 + row1 +k+1;
2200 conn[cnt++] = planep1 + row1p1 +k;
2204 conn[cnt++] = plane + row + np-i-j-1;
2205 conn[cnt++] = plane + row + np-i-j-2;
2206 conn[cnt++] = plane + rowp1 + np-i-j-2;
2207 conn[cnt++] = planep1 + row1 + np-i-j-2;
2212 plane += (np-i)*(np-i+1)/2;
int GetMode(const int i, const int j, const int k)
Compute the mode number in the expansion for a particular tensorial combination.
NekDouble GetConstFactor(const ConstFactorType &factor) const
#define ASSERTL0(condition, msg)
virtual int v_NumBndryCoeffs() const
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
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)
Principle Modified Functions .
MatrixType GetMatrixType() const
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
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)
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
static Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetEdgeNcoeffs(const int i) const
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Principle Modified Functions .
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dx, Array< OneD, NekDouble > &out_dy, Array< OneD, NekDouble > &out_dz)
Calculate the derivative of the physical points.
virtual int v_GetNverts() const
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)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
boost::shared_ptr< DNekMat > DNekMatSharedPtr
virtual int v_NumDGBndryCoeffs() const
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
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)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
virtual int v_GetFaceIntNcoeffs(const int i) const
LibUtilities::ShapeType DetShapeType() const
Gauss Radau pinned at x=-1, .
virtual LibUtilities::ShapeType v_DetShapeType() const
virtual int v_GetTotalFaceIntNcoeffs() const
Principle Orthogonal Functions .
bool ConstFactorExists(const ConstFactorType &factor) const
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
int NumBndryCoeffs(void) const
LibUtilities::BasisType GetEdgeBasisType(const int i) const
This function returns the type of expansion basis on the i-th edge.
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...
The base class for all shapes.
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)
static const NekDouble kNekZeroTol
virtual LibUtilities::PointsKey v_GetFacePointsKey(const int i, const int j) const
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
int getNumberOfCoefficients(int Na)
virtual int v_GetTotalEdgeIntNcoeffs() const
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Principle Modified Functions .
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Principle Orthogonal Functions .
virtual void v_MultiplyByStdQuadratureMetric(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.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
virtual void v_GetFaceToElementMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, int P=-1, int Q=-1)
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
virtual bool v_IsBoundaryInteriorExpansion()
virtual int v_GetNfaces() const
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTBase_MatOp(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...
virtual int v_GetFaceNcoeffs(const int i) const
Gauss Radau pinned at x=-1, .
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
virtual int v_GetFaceNumPoints(const int i) const
boost::shared_ptr< StdTetExp > StdTetExpSharedPtr
int GetNumModes() const
Returns the order of the basis.
void Zero(int n, T *x, const int incx)
Zero vector.
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
virtual int v_GetNedges() const
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Describes the specification for a Basis.
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
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.
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space...