56 int nquad_f = FaceExp->GetNumPoints(0) * FaceExp->GetNumPoints(1);
57 int order_f = FaceExp->GetNcoeffs();
91 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
93 GetPhysFaceVarCoeffsFromElement(face,FaceExp,x->second,varcoeff_work);
94 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp->GetPhys(),1,FaceExp->UpdatePhys(),1);
101 FaceExp->IProductWRTBase(facePhys, outcoeff);
103 for (i = 0; i < order_f; ++i)
105 outarray[(*map)[i].index] += (*map)[i].sign * tau * outcoeff[i];
114 for (n = 0; n < coordim; ++n)
129 Vmath::Vmul(nquad_f, ncdotMF_f, 1, facePhys, 1, inval, 1);
133 Vmath::Vmul(nquad_f, normals[n], 1, facePhys, 1, inval, 1);
137 const NekDouble *data = invMass.GetRawPtr();
150 FaceExp->IProductWRTBase(inval, outcoeff);
153 for (i = 0; i < ncoeffs; ++i)
156 for (j = 0; j < order_f; ++j)
158 tmpcoeff[i] += scale * data[i + (*map)[j].index * ncoeffs] *
159 (*map)[j].sign * outcoeff[j];
167 GetMF(n, coordim, varcoeffs);
177 Coeffs = Coeffs + Dmat * Tmpcoeff;
182 Coeffs = Coeffs + Dmat * Tmpcoeff;
219 for (
unsigned int i = 0; i < FaceExp->GetNcoeffs(); ++i)
221 facetmp[i] = tmp[emap[i]];
225 FaceExp->BwdTrans(facetmp, outarray);
239 int order_f, nquad_f;
243 for (f = 0; f < nfaces; ++f)
245 order_f = FaceExp[f]->GetNcoeffs();
246 nquad_f = FaceExp[f]->GetNumPoints(0) * FaceExp[f]->GetNumPoints(1);
253 for (i = 0; i < order_f; ++i)
255 faceCoeffs[i] = inarray[i + cnt];
259 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
276 StdRegions::VarCoeffMap::const_iterator x;
282 Vmath::Vmul(nquad_f, ncdotMF_f, 1, facePhys, 1, facePhys, 1);
286 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
303 for (f = 0; f < nfaces; ++f)
305 nquad_f = FaceExp[f]->GetNumPoints(0) * FaceExp[f]->GetNumPoints(1);
311 FaceExp[f]->BwdTrans(faceCoeffs[f], facePhys);
313 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
328 int order_f = FaceExp->GetNcoeffs();
336 FaceExp->IProductWRTBase(facePhys, coeff);
339 for (i = 0; i < order_f; ++i)
341 outarray[(*map)[i].index] += (*map)[i].sign * coeff[i];
367 ASSERTL1((*map1).size() == (*map2).size(),
368 "There is an error with the GetTraceToElementMap");
370 for (j = 0; j < (*map1).size(); ++j)
373 for (k = 0; k < (*map2).size(); ++k)
376 if ((*map1)[j].index == (*map2)[k].index && k != j)
380 if ((*map1)[j].
sign != (*map2)[k].
sign)
400 for (i = 0; i < nfaces; ++i)
413 "Geometric information is not set up");
500 int rows = deriv0.GetRows();
501 int cols = deriv1.GetColumns();
505 (*WeakDeriv) = df[3 * dir][0] * deriv0 +
506 df[3 * dir + 1][0] * deriv1 +
507 df[3 * dir + 2][0] * deriv2;
561 int rows = lap00.GetRows();
562 int cols = lap00.GetColumns();
567 (*lap) = gmat[0][0] * lap00 + gmat[4][0] * lap11 +
569 gmat[3][0] * (lap01 +
Transpose(lap01)) +
570 gmat[6][0] * (lap02 +
Transpose(lap02)) +
587 int rows = LapMat.GetRows();
588 int cols = LapMat.GetColumns();
594 (*helm) = LapMat + factor * MassMat;
610 "Need to specify eFactorGJP to construct "
611 "a HelmholtzGJP matrix");
615 factor /= HelmMat.Scale();
617 int ntot = HelmMat.GetRows() * HelmMat.GetColumns();
620 HelmMat.GetRawPtr(), 1, &NDTraceMat->GetPtr()[0], 1);
623 HelmMat.Scale(), NDTraceMat);
671 int rows = LapMat.GetRows();
672 int cols = LapMat.GetColumns();
678 (*adr) = LapMat - lambda * MassMat + AdvMat;
684 if (!massVarcoeffs.empty())
688 if (!lapVarcoeffs.empty())
743 "Need to specify eFactorGJP to construct "
744 "a LinearAdvectionDiffusionReactionGJP matrix");
746 int rows = LapMat.GetRows();
747 int cols = LapMat.GetColumns();
754 LapMat - lambda * MassMat + AdvMat + gjpfactor * NDTraceMat;
898 "Matrix construction is not implemented for variable "
899 "coefficients at the moment");
910 "HybridDGHelmholtz matrix not set up "
911 "for non boundary-interior expansions");
939 StdRegions::VarCoeffMap::const_iterator x;
942 for (i = 0; i < coordim; ++i)
949 GetMF(i, coordim, varcoeffs);
970 Mat = Mat + Dmat * invMass *
Transpose(Dmat);
975 Mat = Mat + Dmat * invMass *
Transpose(Dmat);
1000 Mat = Mat + lambdaval * Mass;
1003 for (i = 0; i < nfaces; ++i)
1006 order_f = FaceExp->GetNcoeffs();
1013 "HDGHelmholtz needs setting up for variable P");
1033 for (j = 0; j < order_f; ++j)
1035 for (k = 0; k < order_f; ++k)
1037 Mat((*map)[j].index, (*map)[k].index) +=
1038 tau * (*map)[j].sign * (*map)[k].sign *
eMass(j, k);
1078 for (i = 0; i < nfaces; ++i)
1088 for (j = 0; j < nface; ++j)
1092 face_lambda[j] = 1.0;
1097 FaceExp->BwdTrans(face_lambda, tmp);
1104 for (k = 0; k < ncoeffs; ++k)
1106 Umat(k, bndry_cnt) = Ulam[k];
1189 for (i = 0; i < nfaces; ++i)
1211 ASSERTL0(
false,
"Direction not known");
1217 StdRegions::VarCoeffMap::const_iterator x;
1224 GetMF(dir, coordim, varcoeffs);
1253 for (j = 0; j < nbndry; ++j)
1259 for (k = 0; k < ncoeffs; ++k)
1261 Ulam[k] = lamToU(k, j);
1279 &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
1287 int order_f, nquad_f;
1333 for (i = 0; i < nfaces; ++i)
1339 for (i = 0; i < nbndry; ++i)
1347 for (f = 0; f < nfaces; ++f)
1349 order_f = FaceExp[f]->GetNcoeffs();
1350 nquad_f = FaceExp[f]->GetNumPoints(0) *
1351 FaceExp[f]->GetNumPoints(1);
1375 for (j = 0; j < order_f; ++j)
1378 (*map)[j].sign * (*LamToQ[0])((*map)[j].index, i);
1381 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1397 0, f, FaceExp[f], normals, varcoeffs);
1399 Vmath::Vmul(nquad_f, ncdotMF, 1, facePhys, 1, work, 1);
1403 Vmath::Vmul(nquad_f, normals[0], 1, facePhys, 1, work,
1408 for (j = 0; j < order_f; ++j)
1411 (*map)[j].sign * (*LamToQ[1])((*map)[j].index, i);
1414 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1430 1, f, FaceExp[f], normals, varcoeffs);
1432 Vmath::Vvtvp(nquad_f, ncdotMF, 1, facePhys, 1, work, 1,
1437 Vmath::Vvtvp(nquad_f, normals[1], 1, facePhys, 1, work,
1442 for (j = 0; j < order_f; ++j)
1445 (*map)[j].sign * (*LamToQ[2])((*map)[j].index, i);
1448 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1464 2, f, FaceExp[f], normals, varcoeffs);
1466 Vmath::Vvtvp(nquad_f, ncdotMF, 1, facePhys, 1, work, 1,
1471 Vmath::Vvtvp(nquad_f, normals[2], 1, facePhys, 1, work,
1477 for (j = 0; j < order_f; ++j)
1480 (*map)[j].sign * LamToU((*map)[j].index, i) -
1484 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
1496 Vmath::Svtvp(nquad_f, -tau, facePhys, 1, work, 1, work, 1);
1499 FaceExp[f]->IProductWRTBase(work, faceCoeffs);
1503 for (j = 0; j < order_f; ++j)
1505 BndMat(cnt + j, i) = faceCoeffs[j];
1521 LapMat.GetRows(), LapMat.GetColumns());
1552 for (
int d = 0;
d < ncoords; ++
d)
1560 for (
int t = 0; t < ntraces; ++t)
1563 tracepts[t] = traceExp[t]->GetTotPoints();
1564 maxtpts = (maxtpts > tracepts[t]) ? maxtpts : tracepts[t];
1570 for (
int t = 0; t < ntraces; ++t)
1579 PhysDeriv(phys, Deriv[0], Deriv[1], Deriv[2]);
1581 for (
int t = 0; t < ntraces; ++t)
1588 traceExp[t]->GetBasis(0)->GetBasisKey();
1590 traceExp[t]->GetBasis(1)->GetBasisKey();
1592 (fromkey0 != tokey0) || (fromkey1 != tokey1);
1599 for (
int d = 0;
d < ncoords; ++
d)
1616 tmp = dphidn[t] + i * tracepts[t], 1,
1617 tmp1 = dphidn[t] + i * tracepts[t], 1);
1622 for (
int t = 0; t < ntraces; ++t)
1624 int nt = tracepts[t];
1630 (
p == 1) ? 0.02 * h * h : 0.8 * pow(
p + 1, -4.0) * h * h;
1637 dphidn[t] + j * nt, 1, val, 1);
1639 Mat(i, j) + scale * traceExp[t]->Integral(val);
1647 for (
int j = 0; j < i; ++j)
1649 Mat(i, j) = Mat(j, i);
1656 "This matrix type cannot be generated from this class");
1688 if (faceExp->GetRightAdjacentElementExp())
1690 if (faceExp->GetRightAdjacentElementExp()
1692 ->GetGlobalID() ==
GetGeom()->GetGlobalID())
1705 int order_e = (*map).size();
1706 int n_coeffs = FaceExp->GetNcoeffs();
1710 if (n_coeffs != order_e)
1712 FaceExp->FwdTrans(Fn, faceCoeffs);
1714 int NumModesElementMax = FaceExp->GetBasis(0)->GetNumModes();
1715 int NumModesElementMin =
m_base[0]->GetNumModes();
1717 FaceExp->ReduceOrderCoeffs(NumModesElementMin, faceCoeffs, faceCoeffs);
1720 FaceExp->DetShapeType(), *FaceExp);
1721 FaceExp->MassMatrixOp(faceCoeffs, faceCoeffs, masskey);
1724 int offset1 = 0, offset2 = 0;
1728 for (i = 0; i < NumModesElementMin; ++i)
1730 for (j = 0; j < NumModesElementMin; ++j)
1732 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1734 offset1 += NumModesElementMin;
1735 offset2 += NumModesElementMax;
1739 for (i = NumModesElementMin; i < NumModesElementMax; ++i)
1741 for (j = NumModesElementMin; j < NumModesElementMax; ++j)
1743 faceCoeffs[i * NumModesElementMax + j] = 0.0;
1752 int offset1 = 0, offset2 = 0;
1754 for (i = 0; i < NumModesElementMin; ++i)
1756 for (j = 0; j < NumModesElementMin - i; ++j)
1758 faceCoeffs[offset1 + j] = faceCoeffs[offset2 + j];
1760 offset1 += NumModesElementMin - i;
1761 offset2 += NumModesElementMax - i;
1767 FaceExp->IProductWRTBase(Fn, faceCoeffs);
1772 for (i = 0; i < order_e; ++i)
1774 outarray[(*map)[i].index] -= (*map)[i].sign * faceCoeffs[i];
1779 for (i = 0; i < order_e; ++i)
1781 outarray[(*map)[i].index] += (*map)[i].sign * faceCoeffs[i];
1817 Out_d = InvMass * Coeffs;
1825 "Not set up for non boundary-interior expansions");
1826 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1827 "Assuming that input matrix was square");
1832 int order_f = faceExp->GetNcoeffs();
1845 DNekScalMat &facemat = *faceExp->GetLocMatrix(mkey);
1862 int rows = inoutmat->GetRows();
1870 int nfvert = faceExp->GetNverts();
1877 GetLinStdExp()->GetTraceToElementMap(face, linmap, linsign,
1886 for (i = 0; i < nfvert; ++i)
1888 fmap = faceExp->GetVertexMap(i,
true);
1892 map[fmap] = linmap[i];
1903 for (i = 0; i < order_f; ++i)
1905 for (j = 0; j < nbndry; ++j)
1907 if (map[i] == bmap[j])
1913 ASSERTL1(j != nbndry,
"Did not find number in map");
1932 ASSERTL1((*map1).size() == (*map2).size(),
1933 "There is an error with the GetTraceToElementMap");
1935 for (i = 0; i < face; ++i)
1940 for (i = 0; i < (*map1).size(); ++i)
1944 for (j = 0; j < (*map2).size(); ++j)
1946 if ((*map1)[i].index == (*map2)[j].index)
1953 ASSERTL2(idx >= 0,
"Index not found");
1955 sign[i] = (*map2)[idx].sign;
1960 ASSERTL0(
false,
"Could not identify matrix type from dimension");
1963 for (i = 0; i < order_f; ++i)
1966 for (j = 0; j < order_f; ++j)
1969 (*inoutmat)(id1, id2) += facemat(i, j) *
sign[i] *
sign[j];
1980 int nVerts, vid1, vid2, vMap1, vMap2;
1987 DNekMat &VertexMat = (*vertexmatrix);
1989 for (vid1 = 0; vid1 < nVerts; ++vid1)
1993 for (vid2 = 0; vid2 < nVerts; ++vid2)
1996 VertexValue = (*r_bnd)(vMap1, vMap2);
1997 VertexMat.SetValue(vid1, vid2, VertexValue);
2001 return vertexmatrix;
2008 int eid, fid, vid, n, i;
2041 int nConnectedEdges = 3;
2042 int nConnectedFaces = 3;
2053 nBndCoeffs, nBndCoeffs, 0.0, storage);
2055 DNekMat &R = (*transformationmatrix);
2060 for (vid = 0; vid < nVerts; ++vid)
2070 int nedgemodesconnected =
2076 int nfacemodesconnected =
2085 for (eid = 0; eid < nConnectedEdges; ++eid)
2087 MatEdgeLocation[eid] =
2089 nmodes = MatEdgeLocation[eid].size();
2094 &edgemodearray[offset], 1);
2103 for (fid = 0; fid < nConnectedFaces; ++fid)
2105 MatFaceLocation[fid] =
2107 nmodes = MatFaceLocation[fid].size();
2112 &facemodearray[offset], 1);
2119 DNekMat &Sveft = (*vertexedgefacetransformmatrix);
2123 DNekMat &Svef = (*vertexedgefacecoupling);
2126 for (n = 0; n < nedgemodesconnected; ++n)
2129 VertexEdgeFaceValue = (*r_bnd)(
GetVertexMap(vid), edgemodearray[n]);
2132 Svef.SetValue(0, n, VertexEdgeFaceValue);
2136 for (n = 0; n < nfacemodesconnected; ++n)
2139 VertexEdgeFaceValue = (*r_bnd)(
GetVertexMap(vid), facemodearray[n]);
2142 Svef.SetValue(0, n + nedgemodesconnected, VertexEdgeFaceValue);
2156 DNekMat &Sefef = (*edgefacecoupling);
2161 for (m = 0; m < nedgemodesconnected; ++m)
2163 for (n = 0; n < nedgemodesconnected; ++n)
2166 EdgeEdgeValue = (*r_bnd)(edgemodearray[n], edgemodearray[m]);
2169 Sefef.SetValue(n, m, EdgeEdgeValue);
2174 for (n = 0; n < nfacemodesconnected; ++n)
2176 for (m = 0; m < nfacemodesconnected; ++m)
2179 FaceFaceValue = (*r_bnd)(facemodearray[n], facemodearray[m]);
2181 Sefef.SetValue(nedgemodesconnected + n, nedgemodesconnected + m,
2187 for (n = 0; n < nedgemodesconnected; ++n)
2189 for (m = 0; m < nfacemodesconnected; ++m)
2192 FaceFaceValue = (*r_bnd)(edgemodearray[n], facemodearray[m]);
2195 Sefef.SetValue(n, nedgemodesconnected + m, FaceFaceValue);
2197 FaceFaceValue = (*r_bnd)(facemodearray[m], edgemodearray[n]);
2200 Sefef.SetValue(nedgemodesconnected + m, n, FaceFaceValue);
2210 Sveft = -Svef * Sefef;
2214 for (n = 0; n < edgemodearray.size(); ++n)
2216 R.SetValue(
GetVertexMap(vid), edgemodearray[n], Sveft(0, n));
2220 for (n = 0; n < facemodearray.size(); ++n)
2223 Sveft(0, n + nedgemodesconnected));
2246 int efCol, efRow, nedgemodes;
2249 nConnectedFaces = 2;
2258 for (eid = 0; eid < nEdges; ++eid)
2269 DNekMat &Mef = (*efedgefacecoupling);
2275 DNekMat &Meff = (*effacefacecoupling);
2281 DNekMat &Meft = (*edgefacetransformmatrix);
2283 int nfacemodesconnected =
2295 Vmath::Vcopy(nedgemodes, &inedgearray[0], 1, &edgemodearray[0], 1);
2301 for (fid = 0; fid < nConnectedFaces; ++fid)
2303 MatFaceLocation[fid] =
2305 nmodes = MatFaceLocation[fid].size();
2310 &facemodearray[offset], 1);
2316 for (n = 0; n < nedgemodes; ++n)
2318 for (m = 0; m < nfacemodesconnected; ++m)
2321 EdgeFaceValue = (*r_bnd)(edgemodearray[n], facemodearray[m]);
2324 Mef.SetValue(n, m, EdgeFaceValue);
2329 for (n = 0; n < nfacemodesconnected; ++n)
2331 for (m = 0; m < nfacemodesconnected; ++m)
2334 FaceFaceValue = (*r_bnd)(facemodearray[n], facemodearray[m]);
2337 Meff.SetValue(n, m, FaceFaceValue);
2351 for (n = 0; n < Meft.GetRows(); ++n)
2353 for (m = 0; m < Meft.GetColumns(); ++m)
2355 R.SetValue(edgemodearray[n], facemodearray[m], Meft(n, m));
2360 for (i = 0; i < R.GetRows(); ++i)
2362 R.SetValue(i, i, 1.0);
2368 return transformationmatrix;
2389 int i, j, n, eid = 0, fid = 0;
2401 DNekMat &InvR = (*inversetransformationmatrix);
2409 int nedgemodestotal = 0;
2410 int nfacemodestotal = 0;
2412 for (eid = 0; eid < nEdges; ++eid)
2415 nedgemodestotal += nedgemodes;
2418 for (fid = 0; fid < nFaces; ++fid)
2421 nfacemodestotal += nfacemodes;
2430 for (eid = 0; eid < nEdges; ++eid)
2438 Vmath::Vcopy(nedgemodes, &edgearray[0], 1, &edgemodearray[offset],
2442 offset += nedgemodes;
2448 for (fid = 0; fid < nFaces; ++fid)
2456 Vmath::Vcopy(nfacemodes, &facearray[0], 1, &facemodearray[offset],
2460 offset += nfacemodes;
2464 for (i = 0; i < nVerts; ++i)
2466 for (j = 0; j < nedgemodestotal; ++j)
2471 for (j = 0; j < nfacemodestotal; ++j)
2475 for (n = 0; n < nedgemodestotal; ++n)
2477 MatrixValue = InvR.GetValue(
GetVertexMap(i), facemodearray[j]) +
2479 R(edgemodearray[n], facemodearray[j]);
2480 InvR.SetValue(
GetVertexMap(i), facemodearray[j], MatrixValue);
2486 for (i = 0; i < nedgemodestotal; ++i)
2488 for (j = 0; j < nfacemodestotal; ++j)
2490 InvR.SetValue(edgemodearray[i], facemodearray[j],
2491 -R(edgemodearray[i], facemodearray[j]));
2495 for (i = 0; i < nCoeffs; ++i)
2497 InvR.SetValue(i, i, 1.0);
2500 return inversetransformationmatrix;
2514 map<int, int> invmap;
2515 for (j = 0; j < nBndCoeffs; ++j)
2517 invmap[bmap[j]] = j;
2533 for (n = 0; n < nEdgeCoeffs; ++n)
2535 edgemaparray[n] = invmap[maparray[n]];
2538 return edgemaparray;
2554 map<int, int> reversemap;
2555 for (j = 0; j < bmap.size(); ++j)
2557 reversemap[bmap[j]] = j;
2573 fOrient = faceOrient;
2589 ASSERTL1(P1 <= locP1,
"Expect value of passed P1 to "
2590 "be lower or equal to face num modes");
2599 ASSERTL1(P2 <= locP2,
"Expect value of passed P2 to "
2600 "be lower or equal to face num modes");
2603 switch (
GetGeom3D()->GetFace(fid)->GetShapeType())
2607 if (((P1 - 3) > 0) && ((P2 - 3) > 0))
2614 for (n = 0; n < P1 - 3; ++n)
2616 for (
int m = 0; m < P2 - 3 - n; ++m, ++cnt)
2618 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2620 cnt1 += locP2 - 3 - n;
2627 if (((P1 - 2) > 0) && ((P2 - 2) > 0))
2634 for (n = 0; n < P2 - 2; ++n)
2636 for (
int m = 0; m < P1 - 2; ++m, ++cnt)
2638 facemaparray[cnt] = reversemap[maparray[cnt1 + m]];
2647 ASSERTL0(
false,
"Invalid shape type.");
2652 return facemaparray;
2671 map<int, int> reversemap;
2672 for (j = 0; j < bmap.size(); ++j)
2674 reversemap[bmap[j]] = j;
2679 for (n = 0; n < nverts; ++n)
2682 vmap[n] = reversemap[id];
2688 for (
int eid = 0; eid < nedges; ++eid)
2702 for (n = 0; n < nEdgeCoeffs; ++n)
2704 edgemaparray[n] = reversemap[maparray[n]];
2706 emap[eid] = edgemaparray;
2712 for (
int fid = 0; fid < nfaces; ++fid)
2726 for (n = 0; n < nFaceCoeffs; ++n)
2728 facemaparray[n] = reversemap[maparray[n]];
2731 fmap[fid] = facemaparray;
2737 return m_geom->GetForient(face);
2755 int nq0 = FaceExp->GetNumPoints(0);
2756 int nq1 = FaceExp->GetNumPoints(1);
2759 int dir0 =
GetGeom3D()->GetDir(face, 0);
2760 int dir1 =
GetGeom3D()->GetDir(face, 1);
2771 Vmath::Gathr(
static_cast<int>(faceids.size()), inarray, faceids, o_tmp);
2788 m_base[dir0]->GetPointsKey(),
m_base[dir1]->GetPointsKey(),
2789 o_tmp.data(), FaceExp->GetBasis(to_id0)->GetPointsKey(),
2790 FaceExp->GetBasis(to_id1)->GetPointsKey(), o_tmp2.data());
2800 if (faceGeom->GetNumVerts() == 3)
2804 m_geom->GetFace(traceid));
2810 m_geom->GetFace(traceid));
2819 if (idmap.size() != nq0 * nq1)
2830 for (
int i = 0; i < nq0 * nq1; ++i)
2837 for (
int j = 0; j < nq1; ++j)
2839 for (
int i = 0; i < nq0; ++i)
2841 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2847 "Case not supposed to be used in this function");
2856 for (
int i = 0; i < nq0 * nq1; ++i)
2864 for (
int j = 0; j < nq1; j++)
2866 for (
int i = 0; i < nq0; ++i)
2868 idmap[j * nq0 + i] = nq0 - 1 - i + j * nq0;
2876 for (
int j = 0; j < nq1; j++)
2878 for (
int i = 0; i < nq0; ++i)
2880 idmap[j * nq0 + i] = nq0 * (nq1 - 1 - j) + i;
2888 for (
int j = 0; j < nq1; j++)
2890 for (
int i = 0; i < nq0; ++i)
2892 idmap[j * nq0 + i] = nq0 * nq1 - 1 - j * nq0 - i;
2900 for (
int i = 0; i < nq0; ++i)
2902 for (
int j = 0; j < nq1; ++j)
2904 idmap[i * nq1 + j] = i + j * nq0;
2912 for (
int i = 0; i < nq0; ++i)
2914 for (
int j = 0; j < nq1; ++j)
2916 idmap[i * nq1 + j] = nq0 - 1 - i + j * nq0;
2924 for (
int i = 0; i < nq0; ++i)
2926 for (
int j = 0; j < nq1; ++j)
2928 idmap[i * nq1 + j] = i + nq0 * (nq1 - 1) - j * nq0;
2936 for (
int i = 0; i < nq0; ++i)
2938 for (
int j = 0; j < nq1; ++j)
2940 idmap[i * nq1 + j] = nq0 * nq1 - 1 - i - j * nq0;
2946 ASSERTL0(
false,
"Unknow orientation");
2965 int nquad_f = FaceExp_f->GetNumPoints(0) * FaceExp_f->GetNumPoints(1);
2968 int nquad0 =
m_base[0]->GetNumPoints();
2969 int nquad1 =
m_base[1]->GetNumPoints();
2970 int nquad2 =
m_base[2]->GetNumPoints();
2971 int nqtot = nquad0 * nquad1 * nquad2;
2983 StdRegions::VarCoeffMap::const_iterator MFdir;
2988 for (
int k = 0; k < coordim; k++)
2990 MFdir = varcoeffs.find(MMFCoeffs[dir * 5 + k]);
2991 tmp = MFdir->second.GetValue();
2995 Vmath::Vvtvp(nquad_f, &tmp_f[0], 1, &normals[k][0], 1, &nFacecdotMF[0],
2996 1, &nFacecdotMF[0], 1);
3006 int nverts = geom->GetFace(traceid)->GetNumVerts();
3009 tn1.
Sub(*(geom->GetFace(traceid)->GetVertex(1)),
3010 *(geom->GetFace(traceid)->GetVertex(0)));
3012 tn2.
Sub(*(geom->GetFace(traceid)->GetVertex(nverts - 1)),
3013 *(geom->GetFace(traceid)->GetVertex(0)));
3015 normal.
Mult(tn1, tn2);
3019 mag = 1.0 /
sqrt(mag);
3025 for (
int i = 0; i < nverts; ++i)
3028 int edgid = geom->GetEdgeNormalToFaceVert(traceid, i);
3031 Dx.
Sub(*(geom->GetEdge(edgid)->GetVertex(0)),
3032 *(geom->GetEdge(edgid)->GetVertex(1)));
3036 h += fabs(normal.
dot(Dx));
3042 int dir0 = geom->GetDir(traceid, 0);
3043 int dir1 = geom->GetDir(traceid, 1);
3045 for (dirn = 0; dirn < 3; ++dirn)
3047 if ((dirn != dir0) && (dirn != dir1))
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
#define sign(a, b)
return the sign(b)*a
Describes the specification for a Basis.
DNekMatSharedPtr v_BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType) override
void v_AddRobinMassMatrix(const int face, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat) override
void SetTraceToGeomOrientation(Array< OneD, NekDouble > &inout)
Align trace orientation with the geometry orientation.
void v_DGDeriv(const int dir, const Array< OneD, const NekDouble > &incoeffs, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &out_d) override
Evaluate coefficients of weak deriviative in the direction dir given the input coefficicents incoeffs...
StdRegions::Orientation v_GetTraceOrient(int face) override
void AddNormTraceInt(const int dir, Array< OneD, ExpansionSharedPtr > &FaceExp, Array< OneD, Array< OneD, NekDouble > > &faceCoeffs, Array< OneD, NekDouble > &outarray)
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Array< OneD, unsigned int > GetTraceInverseBoundaryMap(int fid, StdRegions::Orientation faceOrient=StdRegions::eNoOrientation, int P1=-1, int P2=-1)
std::vector< bool > m_requireNeg
void GetInverseBoundaryMaps(Array< OneD, unsigned int > &vmap, Array< OneD, Array< OneD, unsigned int > > &emap, Array< OneD, Array< OneD, unsigned int > > &fmap)
Array< OneD, NekDouble > GetnFacecdotMF(const int dir, const int face, ExpansionSharedPtr &FaceExp_f, const Array< OneD, const Array< OneD, NekDouble > > &normals, const StdRegions::VarCoeffMap &varcoeffs)
DNekMatSharedPtr v_BuildInverseTransformationMatrix(const DNekScalMatSharedPtr &transformationmatrix) override
Build inverse and inverse transposed transformation matrix: and .
void AddFaceBoundaryInt(const int face, ExpansionSharedPtr &FaceExp, Array< OneD, NekDouble > &facePhys, Array< OneD, NekDouble > &outarray, const StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
void GetPhysFaceVarCoeffsFromElement(const int face, ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &varcoeff, Array< OneD, NekDouble > &outarray)
void v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p) override
void v_ReOrientTracePhysMap(const StdRegions::Orientation orient, Array< OneD, int > &idmap, const int nq0, const int nq1) override
void v_NormVectorIProductWRTBase(const Array< OneD, const Array< OneD, NekDouble > > &Fvec, Array< OneD, NekDouble > &outarray) override
Array< OneD, unsigned int > GetEdgeInverseBoundaryMap(int eid)
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
void v_GenTraceExp(const int traceid, ExpansionSharedPtr &exp) override
void v_AddFaceNormBoundaryInt(const int face, const ExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
void AddHDGHelmholtzFaceTerms(const NekDouble tau, const int edge, Array< OneD, NekDouble > &facePhys, const StdRegions::VarCoeffMap &dirForcing, Array< OneD, NekDouble > &outarray)
DNekMatSharedPtr v_BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd) override
void v_GetTracePhysVals(const int face, const StdRegions::StdExpansionSharedPtr &FaceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient) override
Extract the physical values along face face from inarray into outarray following the local face orien...
void SetFaceToGeomOrientation(const int face, Array< OneD, NekDouble > &inout)
Align face orientation with the geometry orientation.
SpatialDomains::GeometrySharedPtr GetGeom() const
SpatialDomains::GeometrySharedPtr m_geom
void GetTracePhysMap(const int edge, Array< OneD, int > &outarray)
void DropLocMatrix(const LocalRegions::MatrixKey &mkey)
DNekMatSharedPtr BuildVertexMatrix(const DNekScalMatSharedPtr &r_bnd)
Array< OneD, NekDouble > GetMFMag(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
void GetTracePhysVals(const int trace, const StdRegions::StdExpansionSharedPtr &TraceExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient=StdRegions::eNoOrientation)
std::map< int, ExpansionWeakPtr > m_traceExp
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
ExpansionSharedPtr GetTraceExp(const int traceid)
StdRegions::Orientation GetTraceOrient(int trace)
IndexMapValuesSharedPtr GetIndexMap(const IndexMapKey &ikey)
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Array< OneD, NekDouble > GetMFDiv(const int dir, const StdRegions::VarCoeffMap &varcoeffs)
void TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
const NormalVector & GetTraceNormal(const int id)
Array< OneD, NekDouble > GetMF(const int dir, const int shapedim, const StdRegions::VarCoeffMap &varcoeffs)
DNekMatSharedPtr BuildTransformationMatrix(const DNekScalMatSharedPtr &r_bnd, const StdRegions::MatrixType matrixType)
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
boost::call_traits< DataType >::const_reference x() const
boost::call_traits< DataType >::const_reference z() const
boost::call_traits< DataType >::const_reference y() const
void Sub(PointGeom &a, PointGeom &b)
void Mult(PointGeom &a, PointGeom &b)
_this = a x b
NekDouble dot(PointGeom &a)
retun the dot product between this and input a
void UpdatePosition(NekDouble x, NekDouble y, NekDouble z)
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
int GetNedges() const
return the number of edges in 3D expansion
void GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards)
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
void FillMode(const int mode, Array< OneD, NekDouble > &outarray)
This function fills the array outarray with the mode-th mode of the expansion.
void DropLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
int GetTraceNumPoints(const int i) const
This function returns the number of quadrature points belonging to the i-th trace.
int NumBndryCoeffs(void) const
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
const LibUtilities::PointsKeyVector GetPointsKeys() const
int GetTraceIntNcoeffs(const int i) const
int NumDGBndryCoeffs(void) const
DNekScalBlkMatSharedPtr GetLocStaticCondMatrix(const LocalRegions::MatrixKey &mkey)
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
std::shared_ptr< StdExpansion > GetLinStdExp(void) 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...
int GetNtraces() const
Returns the number of trace elements connected to this element.
int GetNverts() const
This function returns the number of vertices of the expansion domain.
void GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
void GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eForwards)
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
void GetTraceNumModes(const int tid, int &numModes0, int &numModes1, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1, bool UseGLL=false) const
This function returns the basis key belonging to the i-th trace.
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
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)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
bool IsBoundaryInteriorExpansion() const
LibUtilities::ShapeType GetShapeType() const
const VarCoeffMap & GetVarCoeffs() const
MatrixType GetMatrixType() const
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
const ConstFactorMap & GetConstFactors() const
const Array< OneD, const NekDouble > & GetVarCoeff(const StdRegions::VarCoeffType &coeff) const
NekDouble GetConstFactor(const ConstFactorType &factor) const
bool ConstFactorExists(const ConstFactorType &factor) const
int getNumberOfCoefficients(int Na, int Nb)
int getNumberOfCoefficients(int Na, int Nb)
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
std::vector< PointsKey > PointsKeyVector
std::shared_ptr< Expansion > ExpansionSharedPtr
std::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
@ eNoGeomType
No type defined.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< Geometry > GeometrySharedPtr
std::shared_ptr< Geometry3D > Geometry3DSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
@ eLinearAdvectionDiffusionReaction
@ eLinearAdvectionDiffusionReactionGJP
@ eInvLaplacianWithUnityMean
static ConstFactorMap NullConstFactorMap
@ eDir1BwdDir2_Dir2BwdDir1
@ eDir1FwdDir1_Dir2FwdDir2
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir1_Dir2BwdDir2
@ eDir1BwdDir1_Dir2FwdDir2
@ eDir1FwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
static VarCoeffMap NullVarCoeffMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
std::vector< double > d(NPUPPER *NPUPPER)
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
static DNekMatSharedPtr NullDNekMatSharedPtr
static Array< OneD, NekDouble > NullNekDouble1DArray
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
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 Gathr(I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)