46 StdPyrExp::StdPyrExp()
66 ASSERTL0(
false,
"order in 'a' direction is higher "
67 "than order in 'c' direction");
71 ASSERTL0(
false,
"order in 'b' direction is higher "
72 "than order in 'c' direction");
90 for (
int i = 2; i <= P; ++i)
96 for (
int i = 2; i <= Q; ++i)
102 for (
int i = 2; i <= P; ++i)
108 for (
int i = 2; i <= Q; ++i)
114 for (
int i = 2; i <= R; ++i)
120 for (
int i = 2; i <= R; ++i)
126 for (
int i = 2; i <= R; ++i)
132 for (
int i = 2; i <= R; ++i)
138 for (
int j = 2; j <= Q; ++j)
140 for (
int i = 2; i <= P; ++i)
147 for (
int i = 2; i <= P; ++i)
149 for (
int j = 1; j <= R-i; ++j)
156 for (
int i = 2; i <= Q; ++i)
158 for (
int j = 1; j <= R-i; ++j)
165 for (
int i = 2; i <= P; ++i)
167 for (
int j = 1; j <= R-i; ++j)
174 for (
int i = 2; i <= Q; ++i)
176 for (
int j = 1; j <= R-i; ++j)
183 for (
int i = 2; i <= P+1; ++i)
185 for (
int j = 1; j <= Q-i+1; ++j)
187 for (
int k = 1; k <= R-i-j+1; ++k)
197 "Duplicate coefficient entries in map");
200 for (it =
m_map.begin(); it !=
m_map.end(); ++it)
202 const int p = it->first.get<0>();
203 const int q = it->first.get<1>();
204 const int r = it->first.get<2>();
205 const int rp = it->first.get<3>();
208 m_idxMap[p] = map<int, map<int, pair<int, int> > >();
213 m_idxMap[p][q] = map<int, pair<int, int> >();
218 m_idxMap[p][q][r] = pair<int, int>(it->second, rp);
260 int Qx =
m_base[0]->GetNumPoints();
261 int Qy =
m_base[1]->GetNumPoints();
262 int Qz =
m_base[2]->GetNumPoints();
270 eta_x =
m_base[0]->GetZ();
271 eta_y =
m_base[1]->GetZ();
272 eta_z =
m_base[2]->GetZ();
276 for (k = 0, n = 0; k < Qz; ++k)
278 for (j = 0; j < Qy; ++j)
280 for (i = 0; i < Qx; ++i, ++n)
282 if (out_dxi1.num_elements() > 0)
283 out_dxi1[n] = 2.0/(1.0 - eta_z[k]) * dEta_bar1[n];
284 if (out_dxi2.num_elements() > 0)
285 out_dxi2[n] = 2.0/(1.0 - eta_z[k]) * dXi2[n];
286 if (out_dxi3.num_elements() > 0)
287 out_dxi3[n] = (1.0+eta_x[i])/(1.0-eta_z[k])*dEta_bar1[n] +
288 (1.0+eta_y[j])/(1.0-eta_z[k])*dXi2[n] + dEta3[n];
323 ASSERTL1(
false,
"input dir is out of range");
372 if (
m_base[0]->Collocation() &&
373 m_base[1]->Collocation() &&
379 inarray, 1, outarray, 1);
398 inarray,outarray,wsp,
409 bool doCheckCollDir0,
410 bool doCheckCollDir1,
411 bool doCheckCollDir2)
413 const int Qx =
m_base[0]->GetNumPoints();
414 const int Qy =
m_base[1]->GetNumPoints();
415 const int Qz =
m_base[2]->GetNumPoints();
422 map<int, map<int, map<int, pair<int, int> > > >
::iterator it_p;
423 map<int, map<int, pair<int, int> > >
::iterator it_q;
429 for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
437 int i ,j, k, s = 0, cnt = 0, cnt2 = 0;
439 for (k = 0; k < Qz; ++k)
446 for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
449 for (it_r = it_q->second.begin(); it_r != it_q->second.end(); ++it_r)
451 sum += inarray[it_r->second.first] * bz[k + Qz*it_r->second.second];
457 for (j = 0; j < Qy; ++j)
466 for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
468 sum += by[j + Qy*it_q->first] * fpq[cnt++];
473 for (i = 0; i < Qx; ++i, ++s)
479 sum += bx[i + Qx*it_p->first] * fp[cnt2++];
481 sum += inarray[4]*(by1*(bx[i] + bx[i+Qx]) + by0*bx[i+Qx]);
540 if (
m_base[0]->Collocation() &&
541 m_base[1]->Collocation() &&
555 bool multiplybyweights)
559 if(multiplybyweights)
576 inarray,outarray,wsp,
588 bool doCheckCollDir0,
589 bool doCheckCollDir1,
590 bool doCheckCollDir2)
593 int Qx =
m_base[0]->GetNumPoints();
594 int Qy =
m_base[1]->GetNumPoints();
595 int Qz =
m_base[2]->GetNumPoints();
601 map<int, map<int, map<int, pair<int, int> > > >
::iterator it_p;
602 map<int, map<int, pair<int, int> > >
::iterator it_q;
610 const int p = it_p->first;
612 for (k = 0; k < Qz; ++k)
614 for (j = 0; j < Qy; ++j)
617 for (i = 0; i < Qx; ++i, ++s)
619 sum += bx[i + Qx*p]*inarray[s];
625 for (it_q = it_p->second.begin(); it_q != it_p->second.end(); ++it_q)
627 const int q = it_q->first;
629 for (k = 0; k < Qz; ++k)
632 for (j = 0; j < Qy; ++j)
634 sum += by[j + Qy*q]*f[j+Qy*k];
639 for (it_r = it_q->second.begin(); it_r != it_q->second.end(); ++it_r)
641 const int rpqr = it_r->second.second;
643 for (k = 0; k < Qz; ++k)
645 sum += bz[k + Qz*rpqr]*fb[k];
648 outarray[it_r->second.first] = sum;
655 for (k = 0; k < Qz; ++k)
657 for (j = 0; j < Qy; ++j)
659 for (i = 0; i < Qx; ++i, ++s)
661 outarray[4] += inarray[s] * bz[k+Qz]*(
690 int nquad0 =
m_base[0]->GetNumPoints();
691 int nquad1 =
m_base[1]->GetNumPoints();
692 int nquad2 =
m_base[2]->GetNumPoints();
693 int nqtot = nquad0*nquad1*nquad2;
706 for(i = 0; i < nquad0; ++i)
708 gfac0[i] = 0.5*(1+z0[i]);
712 for(i = 0; i < nquad1; ++i)
714 gfac1[i] = 0.5*(1+z1[i]);
718 for(i = 0; i < nquad2; ++i)
720 gfac2[i] = 2.0/(1-z2[i]);
724 const int nq01 = nquad0*nquad1;
725 for(i = 0; i < nquad2; ++i)
728 &inarray[0] + i*nq01, 1,
729 &tmp0 [0] + i*nq01, 1);
760 for(i = 0; i < nquad1*nquad2; ++i)
764 tmp0 .get() + i*nquad0, 1);
773 for(i = 0; i < nquad2; ++i)
776 &inarray[0] + i*nq01, 1,
777 &tmp0 [0] + i*nq01, 1);
779 for(i = 0; i < nquad1*nquad2; ++i)
782 &tmp0[0] + i*nquad0, 1,
783 &tmp0[0] + i*nquad0, 1);
806 ASSERTL1(
false,
"input dir is out of range");
831 eta[1] = 2.0*(1.0 + xi[1])/(1.0 - xi[2]) - 1.0;
832 eta[0] = 2.0*(1.0 + xi[0])/(1.0 - xi[2]) - 1.0;
848 for (
int k = 0; k < Qz; ++k )
850 for (
int j = 0; j < Qy; ++j)
852 for (
int i = 0; i < Qx; ++i)
854 int s = i + Qx*(j + Qy*k);
857 xi_y[s] = (1.0 + eta_y[j]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
858 xi_x[s] = (1.0 + etaBar_x[i]) * (1.0 - eta_z[k]) / 2.0 - 1.0;
900 "BasisType is not a boundary interior form");
903 "BasisType is not a boundary interior form");
906 "BasisType is not a boundary interior form");
908 int P =
m_base[0]->GetNumModes();
909 int Q =
m_base[1]->GetNumModes();
910 int R =
m_base[2]->GetNumModes();
918 ASSERTL2(i >= 0 && i <= 7,
"edge id is out of range");
920 if (i == 0 || i == 2)
924 else if (i == 1 || i == 3)
936 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
942 else if (i == 1 || i == 3)
945 return Q+1 + (P*(1 + 2*Q - P))/2;
950 return Q+1 + (P*(1 + 2*Q - P))/2;
956 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
958 int P =
m_base[0]->GetNumModes()-1;
959 int Q =
m_base[1]->GetNumModes()-1;
960 int R =
m_base[2]->GetNumModes()-1;
966 else if (i == 1 || i == 3)
968 return (P-1) * (2*(R-1) - (P-1) - 1) / 2;
972 return (Q-1) * (2*(R-1) - (Q-1) - 1) / 2;
978 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
982 return m_base[0]->GetNumPoints()*
983 m_base[1]->GetNumPoints();
985 else if (i == 1 || i == 3)
987 return m_base[0]->GetNumPoints()*
988 m_base[2]->GetNumPoints();
992 return m_base[1]->GetNumPoints()*
993 m_base[2]->GetNumPoints();
999 const int i,
const int k)
const
1001 ASSERTL2(i >= 0 && i <= 4,
"face id is out of range");
1002 ASSERTL2(k >= 0 && k <= 1,
"basis key id is out of range");
1011 m_base[k]->GetNumModes());
1020 m_base[2*k]->GetNumModes());
1028 m_base[k+1]->GetNumModes());
1037 const std::vector<unsigned int> &nummodes,
1041 nummodes[modes_offset],
1042 nummodes[modes_offset+1],
1043 nummodes[modes_offset+2]);
1051 ASSERTL2(i >= 0 && i <= 7,
"edge id is out of range");
1052 if (i == 0 || i == 2)
1056 else if (i == 1 || i == 3)
1080 "Method only implemented if BasisType is identical"
1081 "in x and y directions");
1084 "Method only implemented for Modified_A BasisType"
1085 "(x and y direction) and Modified_C BasisType (z "
1088 int i, j, k, p, q, r, nFaceCoeffs;
1089 int nummodesA, nummodesB;
1091 int order0 =
m_base[0]->GetNumModes();
1092 int order1 =
m_base[1]->GetNumModes();
1093 int order2 =
m_base[2]->GetNumModes();
1113 bool CheckForZeroedModes =
false;
1123 nFaceCoeffs = P*(2*Q-P+1)/2;
1124 CheckForZeroedModes =
true;
1129 CheckForZeroedModes =
true;
1133 if (maparray.num_elements() != nFaceCoeffs)
1138 if (signarray.num_elements() != nFaceCoeffs)
1144 fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1153 for (i = 0; i < Q; i++)
1155 for (j = 0; j < P; j++)
1159 arrayindx[i*P+j] = i*P+j;
1163 arrayindx[i*P+j] = j*Q+i;
1177 maparray[arrayindx[0]] = 0;
1178 maparray[arrayindx[1]] = 1;
1179 maparray[arrayindx[P+1]] = 2;
1180 maparray[arrayindx[P]] = 3;
1184 for (p = 2; p < P; ++p)
1186 maparray[arrayindx[p]] = p-2 + cnt;
1191 for (q = 2; q < Q; ++q)
1193 maparray[arrayindx[q*P+1]] = q-2 + cnt;
1198 for (p = 2; p < P; ++p)
1200 maparray[arrayindx[P+p]] = p-2 + cnt;
1205 for (q = 2; q < Q; ++q)
1207 maparray[arrayindx[q*P]] = q-2 + cnt;
1211 cnt += Q-2 + 4*(P-2);
1212 for (q = 2; q < Q; ++q)
1214 for (p = 2; p < P; ++p)
1216 maparray[arrayindx[q*P+p]] = cnt + (q-2)*P+(p-2);
1232 for (p = 2; p < P; q += Q-p, ++p)
1234 maparray[q] = cnt++;
1235 if ((
int)faceOrient == 7)
1237 signarray[q] = p % 2 ? -1 : 1;
1242 cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1243 for (q = 2; q < Q; ++q)
1245 maparray[q] = cnt++;
1246 if ((
int)faceOrient == 7)
1248 signarray[q] = q % 2 ? -1 : 1;
1253 cnt = 5 + 2*(order0-2) + 2*(order1-2);
1254 for (q = 2; q < Q; ++q)
1256 maparray[Q+q-1] = cnt++;
1257 if ((
int)faceOrient == 7)
1259 signarray[Q+q-1] = q % 2 ? -1 : 1;
1264 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1267 for (p = 2; p < P; ++p)
1269 for (r = 2; r < Q-p; ++r)
1271 maparray[cnt2] = cnt++;
1272 if ((
int)faceOrient == 7 && p > 1)
1274 signarray[cnt2++] = p % 2 ? -1 : 1;
1288 cnt = 5 + (order0-2);
1290 for (p = 2; p < P; q += Q-p, ++p)
1292 maparray[q] = cnt++;
1293 if ((
int)faceOrient == 7)
1295 signarray[q] = p % 2 ? -1 : 1;
1300 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1301 for (q = 2; q < Q; ++q)
1303 maparray[q] = cnt++;
1304 if ((
int)faceOrient == 7)
1306 signarray[q] = q % 2 ? -1 : 1;
1311 cnt = 5 + 2*(order0-2) + 2*(order1-2) + (order2-2);
1312 for (q = 2; q < Q; ++q)
1314 maparray[Q+q-1] = cnt++;
1315 if ((
int)faceOrient == 7)
1317 signarray[Q+q-1] = q % 2 ? -1 : 1;
1322 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1325 for (p = 2; p < P; ++p)
1327 for (r = 2; r < Q-p; ++r)
1329 maparray[cnt2] = cnt++;
1330 if ((
int)faceOrient == 7 && p > 1)
1332 signarray[cnt2++] = p % 2 ? -1 : 1;
1346 cnt = 5 + (order0-2) + (order1-2);
1348 for (p = 2; p < P; q += Q-p, ++p)
1350 maparray[q] = cnt++;
1351 if ((
int)faceOrient == 7)
1353 signarray[q] = p % 2 ? -1 : 1;
1358 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 2*(order2-2);
1359 for (q = 2; q < Q; ++q)
1361 maparray[q] = cnt++;
1362 if ((
int)faceOrient == 7)
1364 signarray[q] = q % 2 ? -1 : 1;
1369 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1370 for (q = 2; q < Q; ++q)
1372 maparray[Q+q-1] = cnt++;
1373 if ((
int)faceOrient == 7)
1375 signarray[Q+q-1] = q % 2 ? -1 : 1;
1380 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1384 for (p = 2; p < P; ++p)
1386 for (r = 2; r < Q-p; ++r)
1388 maparray[cnt2] = cnt++;
1389 if ((
int)faceOrient == 7 && p > 1)
1391 signarray[cnt2++] = p % 2 ? -1 : 1;
1405 cnt = 5 + 2*(order0-2) + (order1-2);
1407 for (p = 2; p < P; q += Q-p, ++p)
1409 maparray[q] = cnt++;
1410 if ((
int)faceOrient == 7)
1412 signarray[q] = p % 2 ? -1 : 1;
1417 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 3*(order2-2);
1418 for (q = 2; q < Q; ++q)
1420 maparray[q] = cnt++;
1421 if ((
int)faceOrient == 7)
1423 signarray[q] = q % 2 ? -1 : 1;
1428 cnt = 5 + 2*(order0-2) + 2*(order1-2);
1429 for (q = 2; q < Q; ++q)
1431 maparray[Q+q-1] = cnt++;
1432 if ((
int)faceOrient == 7)
1434 signarray[Q+q-1] = q % 2 ? -1 : 1;
1439 cnt = 5 + 2*(order0-2) + 2*(order1-2) + 4*(order2-2)
1443 for (p = 2; p < P; ++p)
1445 for (r = 2; r < Q-p; ++r)
1447 maparray[cnt2] = cnt++;
1448 if ((
int)faceOrient == 7 && p > 1)
1450 signarray[cnt2++] = p % 2 ? -1 : 1;
1458 ASSERTL0(
false,
"Face to element map unavailable.");
1464 if(CheckForZeroedModes)
1469 for (j = 0; j < P; ++j)
1472 for (k = Q-j; k < Q-j; ++k)
1474 signarray[idx] = 0.0;
1475 maparray[idx++] = maparray[0];
1479 for (j = P; j < P; ++j)
1481 for (k = 0; k < Q-j; ++k)
1483 signarray[idx] = 0.0;
1484 maparray[idx++] = maparray[0];
1491 if ((
int)faceOrient == 7)
1493 swap(maparray[0], maparray[P]);
1494 for (i = 1; i < P-1; ++i)
1496 swap(maparray[i+1], maparray[P+i]);
1502 if(CheckForZeroedModes)
1506 for (j = 0; j < P; ++j)
1508 for (k = Q; k < Q; ++k)
1510 signarray[arrayindx[j+k*P]] = 0.0;
1511 maparray[arrayindx[j+k*P]] = maparray[0];
1515 for (j = P; j < P; ++j)
1517 for (k = 0; k < Q; ++k)
1519 signarray[arrayindx[j+k*P]] = 0.0;
1520 maparray[arrayindx[j+k*P]] = maparray[0];
1529 if (faceOrient == 6 || faceOrient == 8 ||
1530 faceOrient == 11 || faceOrient == 12)
1534 for (i = 3; i < Q; i += 2)
1536 for (j = 0; j < P; j++)
1538 signarray[arrayindx[i*P+j]] *= -1;
1542 for (i = 0; i < P; i++)
1544 swap(maparray [i], maparray [i+P]);
1545 swap(signarray[i], signarray[i+P]);
1550 for (i = 0; i < Q; i++)
1552 for (j = 3; j < P; j += 2)
1554 signarray[arrayindx[i*P+j]] *= -1;
1558 for (i = 0; i < Q; i++)
1560 swap (maparray [i], maparray [i+Q]);
1561 swap (signarray[i], signarray[i+Q]);
1566 if (faceOrient == 7 || faceOrient == 8 ||
1567 faceOrient == 10 || faceOrient == 12)
1571 for (i = 0; i < Q; i++)
1573 for (j = 3; j < P; j += 2)
1575 signarray[arrayindx[i*P+j]] *= -1;
1579 for(i = 0; i < Q; i++)
1581 swap(maparray [i*P], maparray [i*P+1]);
1582 swap(signarray[i*P], signarray[i*P+1]);
1587 for (i = 3; i < Q; i += 2)
1589 for (j = 0; j < P; j++)
1591 signarray[arrayindx[i*P+j]] *= -1;
1595 for (i = 0; i < P; i++)
1597 swap(maparray [i*Q], maparray [i*Q+1]);
1598 swap(signarray[i*Q], signarray[i*Q+1]);
1610 "Mapping not defined for this type of basis");
1622 const int P =
m_base[0]->GetNumModes() - 2;
1623 const int Q =
m_base[1]->GetNumModes() - 2;
1624 const int R =
m_base[2]->GetNumModes() - 2;
1627 if (maparray.num_elements() != nEdgeIntCoeffs)
1632 if(signarray.num_elements() != nEdgeIntCoeffs)
1638 fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
1664 offset += 2*(P+Q)+R;
1667 offset += 2*(P+Q+R);
1670 offset += 2*(P+Q)+3*R;
1673 ASSERTL0(
false,
"Edge not defined.");
1677 for (i = 0; i < nEdgeIntCoeffs; ++i)
1679 maparray[i] = offset + i;
1684 for (i = 1; i < nEdgeIntCoeffs; i += 2)
1697 const int P =
m_base[0]->GetNumModes() - 1;
1698 const int Q =
m_base[1]->GetNumModes() - 1;
1699 const int R =
m_base[2]->GetNumModes() - 1;
1701 int p, q, r, idx = 0;
1706 if (maparray.num_elements() != nFaceIntCoeffs)
1711 if (signarray.num_elements() != nFaceIntCoeffs)
1717 fill(signarray.get(), signarray.get()+nFaceIntCoeffs, 1);
1728 for (i = 0; i < nummodesB; i++)
1730 for (j = 0; j < nummodesA; j++)
1734 arrayindx[i*nummodesA+j] = i*nummodesA+j;
1738 arrayindx[i*nummodesA+j] = j*nummodesB+i;
1744 int offset = 5 + 2*(P-1) + 2*(Q-1) + 4*(R-1);
1746 for (i = 0; i < fid; ++i)
1754 for (q = 2; q <= Q; ++q)
1756 for (p = 2; p <= P; ++p)
1758 maparray[arrayindx[(q-2)*nummodesA+(p-2)]]
1759 = offset + (q-2)*nummodesA+(p-2);
1766 for (p = 2; p <= P; ++p)
1768 for (r = 1; r <= R-p; ++r, ++idx)
1770 if ((
int)faceOrient == 7)
1772 signarray[idx] = p % 2 ? -1 : 1;
1774 maparray[idx] = offset + idx;
1781 for (q = 2; q <= Q; ++q)
1783 for (r = 1; r <= R-q; ++r, ++idx)
1785 if ((
int)faceOrient == 7)
1787 signarray[idx] = q % 2 ? -1 : 1;
1789 maparray[idx] = offset + idx;
1795 ASSERTL0(
false,
"Face interior map not available.");
1805 if (faceOrient == 6 || faceOrient == 8 ||
1806 faceOrient == 11 || faceOrient == 12)
1810 for (i = 1; i < nummodesB; i += 2)
1812 for (j = 0; j < nummodesA; j++)
1814 signarray[arrayindx[i*nummodesA+j]] *= -1;
1820 for (i = 0; i < nummodesB; i++)
1822 for (j = 1; j < nummodesA; j += 2)
1824 signarray[arrayindx[i*nummodesA+j]] *= -1;
1830 if (faceOrient == 7 || faceOrient == 8 ||
1831 faceOrient == 10 || faceOrient == 12)
1835 for (i = 0; i < nummodesB; i++)
1837 for (j = 1; j < nummodesA; j += 2)
1839 signarray[arrayindx[i*nummodesA+j]] *= -1;
1845 for (i = 1; i < nummodesB; i += 2)
1847 for (j = 0; j < nummodesA; j++)
1849 signarray[arrayindx[i*nummodesA+j]] *= -1;
1860 "BasisType is not a boundary interior form");
1863 "BasisType is not a boundary interior form");
1866 "BasisType is not a boundary interior form");
1871 if (outarray.num_elements() != nIntCoeffs)
1878 for (p = nBndCoeffs; p <
m_ncoeffs; ++p)
1880 outarray[idx++] = p;
1888 "BasisType is not a boundary interior form");
1891 "BasisType is not a boundary interior form");
1894 "BasisType is not a boundary interior form");
1898 for (idx = 0; idx < nBndry; ++idx)
1900 maparray[idx] = idx;
1925 const int R =
m_base[2]->GetNumModes();
1927 for (i = 0; i < I; ++i)
1929 cnt += (R-i)*(R-i+1)/2;
1933 for (j = 0; j < J; ++j)
1948 int nquad0 =
m_base[0]->GetNumPoints();
1949 int nquad1 =
m_base[1]->GetNumPoints();
1950 int nquad2 =
m_base[2]->GetNumPoints();
1959 for(i = 0; i < nquad1*nquad2; ++i)
1962 w0.get(), 1, outarray.get()+i*nquad0,1);
1966 for(j = 0; j < nquad2; ++j)
1968 for(i = 0; i < nquad1; ++i)
1970 Blas::Dscal(nquad0,w1[i], &outarray[0]+i*nquad0 +
1982 for(i = 0; i < nquad2; ++i)
1984 Blas::Dscal(nquad0*nquad1, 0.25*w2[i],
1985 &outarray[0]+i*nquad0*nquad1, 1);
1990 for(i = 0; i < nquad2; ++i)
1992 Blas::Dscal(nquad0*nquad1,0.125*(1-z2[i])*(1-z2[i])*w2[i],
1993 &outarray[0]+i*nquad0*nquad1,1);
virtual LibUtilities::ShapeType v_DetShapeType() const
virtual void v_GetFaceInteriorMap(const int fid, const Orientation faceOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
#define ASSERTL0(condition, msg)
Principle Modified Functions .
virtual int v_GetEdgeNcoeffs(const int i) const
virtual void v_GetCoords(Array< OneD, NekDouble > &xi_x, Array< OneD, NekDouble > &xi_y, Array< OneD, NekDouble > &xi_z)
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
static Array< OneD, NekDouble > NullNekDouble1DArray
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
virtual const LibUtilities::BasisKey v_DetFaceBasisKey(const int i, const int k) const
void MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
virtual int v_GetNedges() 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)
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Principle Modified Functions .
virtual void v_GetEdgeInteriorMap(const int eid, const Orientation edgeOrient, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray)
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
boost::shared_ptr< DNekMat > DNekMatSharedPtr
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
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)
Calculate the derivative of the physical points.
virtual int v_GetFaceNcoeffs(const int i) const
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
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...
virtual int v_GetNfaces() const
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)
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
static const NekDouble kNekZeroTol
virtual int v_GetFaceIntNcoeffs(const int i) 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)
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
boost::tuple< unsigned int, unsigned int, unsigned int, unsigned int > Mode
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
virtual int v_NumBndryCoeffs() const
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
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_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)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
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)
std::map< Mode, unsigned int, cmpop > m_map
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
int getNumberOfCoefficients(int Na, int Nb, int Nc)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
int GetFaceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th face.
int GetTetMode(int I, int J, int K)
Number tetrahedral modes in r-direction. Much the same as StdTetExp::GetTetMode but slightly simplifi...
virtual int v_GetFaceNumPoints(const int i) const
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Gauss Radau pinned at x=-1, .
virtual LibUtilities::BasisType v_GetEdgeBasisType(const int i) const
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Inner product of inarray over region with respect to the expansion basis m_base[0]->GetBdata(),m_base[1]->GetBdata(), m_base[2]->GetBdata() and return in outarray.
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
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)
int GetNumModes() const
Returns the order of the basis.
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual int v_GetNverts() const
#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 void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward transformation is evaluated at the quadrature points.
Describes the specification for a Basis.
std::map< int, std::map< int, std::map< int, std::pair< int, int > > > > m_idxMap
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.