47 namespace LocalRegions
51 void Expansion3D::AddHDGHelmholtzFaceTerms(
60 int nquad_f = FaceExp->GetNumPoints(0)*FaceExp->GetNumPoints(1);
61 int order_f = FaceExp->GetNcoeffs();
62 int coordim = GetCoordim();
63 int ncoeffs = GetNcoeffs();
70 = GetFaceNormal(face);
79 GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
80 face, GetForient(face));
82 StdExpansion::GetIndexMap(ikey);
96 if ((x = varcoeffs.find(VarCoeff[0])) != varcoeffs.end())
98 GetPhysFaceVarCoeffsFromElement(face,FaceExp,x->second,varcoeff_work);
99 Vmath::Vmul(nquad_f,varcoeff_work,1,FaceExp->GetPhys(),1,FaceExp->UpdatePhys(),1);
106 FaceExp->IProductWRTBase(facePhys, outcoeff);
108 for(i = 0; i < order_f; ++i)
110 outarray[(*map)[i].index] += (*map)[i].sign*tau*outcoeff[i];
115 const NekDouble *data = invMass.GetRawPtr();
122 for(n = 0; n < coordim; ++n)
124 Vmath::Vmul(nquad_f, normals[n], 1, facePhys, 1, inval, 1);
126 if (m_negatedNormals[face])
142 FaceExp->IProductWRTBase(inval, outcoeff);
145 for(i = 0; i < ncoeffs; ++i)
148 for(j = 0; j < order_f; ++j)
150 tmpcoeff[i] += scale*data[i+(*map)[j].index*ncoeffs]*(*map)[j].sign*outcoeff[j];
155 Coeffs = Coeffs + Dmat*Tmpcoeff;
178 void Expansion3D::AddNormTraceInt(
187 int nfaces = GetNfaces();
190 for(f = 0; f < nfaces; ++f)
192 order_f = FaceExp[f]->GetNcoeffs();
193 nquad_f = FaceExp[f]->GetNumPoints(0)*FaceExp[f]->GetNumPoints(1);
199 for(i = 0; i < order_f; ++i)
201 faceCoeffs[i] = inarray[i+cnt];
205 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
221 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
223 if (m_negatedNormals[f])
228 AddFaceBoundaryInt(f, FaceExp[f], facePhys, outarray, varcoeffs);
233 void Expansion3D::AddNormTraceInt(
240 int order_f, nquad_f;
241 int nfaces = GetNfaces();
244 for(f = 0; f < nfaces; ++f)
246 order_f = FaceExp[f]->GetNcoeffs();
247 nquad_f = FaceExp[f]->GetNumPoints(0)*FaceExp[f]->GetNumPoints(1);
254 FaceExp[f]->BwdTrans(faceCoeffs[f], facePhys);
256 Vmath::Vmul(nquad_f, normals[dir], 1, facePhys, 1, facePhys, 1);
258 if (m_negatedNormals[f])
263 AddFaceBoundaryInt(f, FaceExp[f], facePhys, outarray);
270 void Expansion3D::AddFaceBoundaryInt(
278 int order_f = FaceExp->GetNcoeffs();
283 GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
284 face, GetForient(face));
286 StdExpansion::GetIndexMap(ikey);
301 FaceExp->IProductWRTBase(facePhys, coeff);
304 for(i = 0; i < order_f; ++i)
306 outarray[(*map)[i].index] += (*map)[i].sign*coeff[i];
313 void Expansion3D::SetFaceToGeomOrientation(
317 int nface = GetFaceNcoeffs(face);
325 GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
328 StdExpansion::GetIndexMap(ikey1);
331 GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
332 face, GetForient(face));
334 StdExpansion::GetIndexMap(ikey2);
336 ASSERTL1((*map1).num_elements() == (*map2).num_elements(),
337 "There is an error with the GetFaceToElementMap");
339 for(j = 0; j < (*map1).num_elements(); ++j)
342 for(k = 0; k < (*map2).num_elements(); ++k)
345 if((*map1)[j].index == (*map2)[k].index && k != j)
349 if((*map1)[j].sign != (*map2)[k].sign)
363 int nfaces = GetNfaces();
367 for(i = 0; i < nfaces; ++i)
369 SetFaceToGeomOrientation(i, f_tmp = inout + cnt);
370 cnt += GetFaceNcoeffs(i);
387 "Matrix construction is not implemented for variable "
388 "coefficients at the moment");
398 ASSERTL1(IsBoundaryInteriorExpansion(),
399 "HybridDGHelmholtz matrix not set up "
400 "for non boundary-interior expansions");
405 int ncoeffs = GetNcoeffs();
406 int nfaces = GetNfaces();
413 int order_f, coordim = GetCoordim();
427 for(i=0; i < coordim; ++i)
430 Mat = Mat + Dmat*invMass*
Transpose(Dmat);
454 Mat = Mat + lambdaval*Mass;
457 for(i = 0; i < nfaces; ++i)
459 FaceExp = GetFaceExp(i);
460 order_f = FaceExp->GetNcoeffs();
463 GetBasisNumModes(0), GetBasisNumModes(1),
464 GetBasisNumModes(2), i, GetForient(i));
466 StdExpansion::GetIndexMap(ikey);
486 for(j = 0; j < order_f; ++j)
488 for(k = 0; k < order_f; ++k)
490 Mat((*map)[j].index,(*map)[k].index) +=
491 tau*(*map)[j].sign*(*map)[k].sign*
eMass(j,k);
502 int nbndry = NumDGBndryCoeffs();
503 int ncoeffs = GetNcoeffs();
504 int nfaces = GetNfaces();
528 for(i = 0; i < nfaces; ++i)
530 FaceExp = GetFaceExp(i);
531 int nface = GetFaceNcoeffs(i);
537 for(j = 0; j < nface; ++j)
541 face_lambda[j] = 1.0;
543 SetFaceToGeomOrientation(i, face_lambda);
546 FaceExp->BwdTrans(face_lambda, tmp);
547 AddHDGHelmholtzFaceTerms(tau, i, tmp, mkey.
GetVarCoeffs(), f);
552 for(k = 0; k < ncoeffs; ++k)
554 Umat(k,bndry_cnt) = Ulam[k];
605 int nbndry = NumDGBndryCoeffs();
607 int ncoeffs = GetNcoeffs();
608 int nfaces = GetNfaces();
630 for(i = 0; i < nfaces; ++i)
632 FaceExp[i] = GetFaceExp(i);
652 ASSERTL0(
false,
"Direction not known");
659 for(j = 0; j < nbndry; ++j)
665 for(k = 0; k < ncoeffs; ++k)
667 Ulam[k] = lamToU(k,j);
674 SetTraceToGeomOrientation(lambda);
678 AddNormTraceInt(dir,lambda,FaceExp,f,mkey.
GetVarCoeffs());
684 Vmath::Vcopy(ncoeffs,&ulam[0],1,&(Qmat.GetPtr())[0]+j*ncoeffs,1);
692 int order_f, nquad_f;
693 int nbndry = NumDGBndryCoeffs();
694 int nfaces = GetNfaces();
717 LamToQ[0] = GetLocMatrix(LamToQ0key);
721 LamToQ[1] = GetLocMatrix(LamToQ1key);
725 LamToQ[2] = GetLocMatrix(LamToQ2key);
728 for(i = 0; i < nfaces; ++i)
730 FaceExp[i] = GetFaceExp(i);
734 for(i = 0; i < nbndry; ++i)
740 SetTraceToGeomOrientation(lam);
742 for(f = 0; f < nfaces; ++f)
744 order_f = FaceExp[f]->GetNcoeffs();
745 nquad_f = FaceExp[f]->GetNumPoints(0)*FaceExp[f]->GetNumPoints(1);
746 normals = GetFaceNormal(f);
753 GetBasisNumModes(0), GetBasisNumModes(1),
754 GetBasisNumModes(2), f, GetForient(f));
756 StdExpansion::GetIndexMap(ikey);
770 for(j = 0; j < order_f; ++j)
772 faceCoeffs[j] = (*map)[j].sign*(*LamToQ[0])((*map)[j].index,i);
775 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
787 Vmath::Vmul(nquad_f, normals[0], 1, facePhys, 1, work, 1);
790 for(j = 0; j < order_f; ++j)
792 faceCoeffs[j] = (*map)[j].sign*(*LamToQ[1])((*map)[j].index,i);
795 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
807 Vmath::Vvtvp(nquad_f, normals[1], 1, facePhys, 1, work, 1, work, 1);
810 for(j = 0; j < order_f; ++j)
812 faceCoeffs[j] = (*map)[j].sign*(*LamToQ[2])((*map)[j].index,i);
815 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
830 if (m_negatedNormals[f])
837 for(j = 0; j < order_f; ++j)
839 faceCoeffs[j] = (*map)[j].sign*LamToU((*map)[j].index,i) - lam[cnt+j];
842 FaceExp[f]->BwdTrans(faceCoeffs, facePhys);
858 FaceExp[f]->IProductWRTBase(work, faceCoeffs);
860 SetFaceToGeomOrientation(f, faceCoeffs);
862 for(j = 0; j < order_f; ++j)
864 BndMat(cnt+j,i) = faceCoeffs[j];
884 int nq = GetTotPoints();
888 IProductWRTBase(tmp, outarray);
891 &(lmat->GetPtr())[0],1);
898 ASSERTL0(
false,
"This matrix type cannot be generated from this class");
907 int nFaces = GetNfaces();
908 ASSERTL1(face < nFaces,
"Face is out of range.");
909 if (m_faceExp.size() < nFaces)
911 m_faceExp.resize(nFaces);
918 return m_faceExp[face].lock();
921 void Expansion3D::v_AddFaceNormBoundaryInt(
938 if (m_requireNeg.size() == 0)
940 m_requireNeg.resize(GetNfaces());
942 for (i = 0; i < GetNfaces(); ++i)
944 m_requireNeg[i] =
false;
945 if (m_negatedNormals[i])
947 m_requireNeg[i] =
true;
953 if (faceExp->GetRightAdjacentElementExp())
955 if (faceExp->GetRightAdjacentElementExp()->GetGeom3D()
956 ->GetGlobalID() == GetGeom3D()->GetGlobalID())
958 m_requireNeg[i] =
true;
966 GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
967 face, GetForient(face));
969 StdExpansion::GetIndexMap(ikey);
971 int order_e = (*map).num_elements();
972 int n_coeffs = FaceExp->GetNcoeffs();
976 if (n_coeffs != order_e)
981 FaceExp->FwdTrans(Fn, faceCoeffs);
983 int NumModesElementMax = FaceExp->GetBasis(0)->GetNumModes();
984 int NumModesElementMin = m_base[0]->GetNumModes();
986 FaceExp->ReduceOrderCoeffs(NumModesElementMin,
991 FaceExp->DetShapeType(),
993 FaceExp->MassMatrixOp(faceCoeffs, faceCoeffs, masskey);
996 int offset1 = 0, offset2 = 0;
1000 for (i = 0; i < NumModesElementMin; ++i)
1002 for (j = 0; j < NumModesElementMin; ++j)
1004 faceCoeffs[offset1+j] =
1005 faceCoeffs[offset2+j];
1007 offset1 += NumModesElementMin;
1008 offset2 += NumModesElementMax;
1012 for (i = NumModesElementMin; i < NumModesElementMax; ++i)
1014 for (j = NumModesElementMin; j < NumModesElementMax; ++j)
1016 faceCoeffs[i*NumModesElementMax+j] = 0.0;
1025 int offset1 = 0, offset2 = 0;
1027 for (i = 0; i < NumModesElementMin; ++i)
1029 for (j = 0; j < NumModesElementMin-i; ++j)
1031 faceCoeffs[offset1+j] =
1032 faceCoeffs[offset2+j];
1034 offset1 += NumModesElementMin-i;
1035 offset2 += NumModesElementMax-i;
1042 FaceExp->IProductWRTBase(Fn, faceCoeffs);
1045 if (m_requireNeg[face])
1047 for (i = 0; i < order_e; ++i)
1049 outarray[(*map)[i].index] -= (*map)[i].sign * faceCoeffs[i];
1054 for (i = 0; i < order_e; ++i)
1056 outarray[(*map)[i].index] += (*map)[i].sign * faceCoeffs[i];
1067 void Expansion3D::v_DGDeriv(
1074 int ncoeffs = GetNcoeffs();
1080 DNekScalMat &Dmat = *GetLocMatrix(DerivType[dir]);
1090 AddNormTraceInt(dir, FaceExp, faceCoeffs, coeffs);
1094 Out_d = InvMass*Coeffs;
1097 void Expansion3D::v_AddRobinMassMatrix(
1102 ASSERTL1(IsBoundaryInteriorExpansion(),
1103 "Not set up for non boundary-interior expansions");
1104 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
1105 "Assuming that input matrix was square");
1110 int order_f = faceExp->GetNcoeffs();
1119 faceExp->DetShapeType();
1127 DNekScalMat &facemat = *faceExp->GetLocMatrix(mkey);
1141 int rows = inoutmat->GetRows();
1143 if (rows == GetNcoeffs())
1145 GetFaceToElementMap(face,GetForient(face),map,sign);
1147 else if(rows == NumBndryCoeffs())
1149 int nbndry = NumBndryCoeffs();
1152 GetFaceToElementMap(face,GetForient(face),map,sign);
1153 GetBoundaryMap(bmap);
1155 for(i = 0; i < order_f; ++i)
1157 for(j = 0; j < nbndry; ++j)
1159 if(map[i] == bmap[j])
1165 ASSERTL1(j != nbndry,
"Did not find number in map");
1168 else if (rows == NumDGBndryCoeffs())
1177 GetBasisNumModes(0), GetBasisNumModes(1), GetBasisNumModes(2),
1178 face, GetForient(face));
1180 StdExpansion::GetIndexMap(ikey1);
1184 GetBasisNumModes(0),
1185 GetBasisNumModes(1),
1186 GetBasisNumModes(2),
1190 StdExpansion::GetIndexMap(ikey2);
1192 ASSERTL1((*map1).num_elements() == (*map2).num_elements(),
1193 "There is an error with the GetFaceToElementMap");
1195 for (i = 0; i < face; ++i)
1197 cnt += GetFaceNcoeffs(i);
1200 for(i = 0; i < (*map1).num_elements(); ++i)
1204 for(j = 0; j < (*map2).num_elements(); ++j)
1206 if((*map1)[i].index == (*map2)[j].index)
1213 ASSERTL2(idx >= 0,
"Index not found");
1214 map [i] = idx + cnt;
1215 sign[i] = (*map2)[idx].sign;
1220 ASSERTL0(
false,
"Could not identify matrix type from dimension");
1223 for(i = 0; i < order_f; ++i)
1226 for(j = 0; j < order_f; ++j)
1229 (*inoutmat)(id1,id2) += facemat(i,j)*sign[i]*sign[j];
1240 int nVerts, vid1, vid2, vMap1, vMap2;
1243 nVerts = GetNverts();
1247 nVerts, nVerts, 0.0, storage);
1248 DNekMat &VertexMat = (*m_vertexmatrix);
1250 for (vid1 = 0; vid1 < nVerts; ++vid1)
1252 vMap1 = GetVertexMap(vid1,
true);
1254 for (vid2 = 0; vid2 < nVerts; ++vid2)
1256 vMap2 = GetVertexMap(vid2,
true);
1257 VertexValue = (*r_bnd)(vMap1, vMap2);
1258 VertexMat.SetValue(vid1, vid2, VertexValue);
1262 return m_vertexmatrix;
1270 int eid, fid, vid, n, i;
1272 int nBndCoeffs = NumBndryCoeffs();
1277 nVerts = GetNverts();
1278 nEdges = GetNedges();
1303 int nConnectedEdges = 3;
1304 int nConnectedFaces = 3;
1308 MatEdgeLocation(nConnectedEdges);
1310 MatFaceLocation(nConnectedFaces);
1317 m_transformationmatrix =
1319 nBndCoeffs, nBndCoeffs, 0.0, storage);
1320 m_transposedtransformationmatrix =
1322 nBndCoeffs, nBndCoeffs, 0.0, storage);
1324 DNekMat &R = (*m_transformationmatrix);
1325 DNekMat &RT = (*m_transposedtransformationmatrix);
1330 for (vid = 0; vid < nVerts; ++vid)
1334 GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 0)) +
1335 GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 1)) +
1336 GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 2)) +
1337 GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
1338 GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
1339 GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 2)) - 6;
1341 int nedgemodesconnected =
1342 GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 0)) +
1343 GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 1)) +
1344 GetEdgeNcoeffs (geom->GetVertexEdgeMap(vid, 2)) - 6;
1347 int nfacemodesconnected =
1348 GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 0)) +
1349 GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 1)) +
1350 GetFaceIntNcoeffs(geom->GetVertexFaceMap(vid, 2));
1356 for (eid = 0; eid < nConnectedEdges; ++eid)
1358 MatEdgeLocation[eid] = GetEdgeInverseBoundaryMap(
1359 geom->GetVertexEdgeMap(vid, eid));
1360 nmodes = MatEdgeLocation[eid].num_elements();
1365 1, &edgemodearray[offset], 1);
1374 for (fid = 0; fid < nConnectedFaces; ++fid)
1376 MatFaceLocation[fid] = GetFaceInverseBoundaryMap(
1377 geom->GetVertexFaceMap(vid, fid));
1378 nmodes = MatFaceLocation[fid].num_elements();
1383 1, &facemodearray[offset], 1);
1390 1, efRow, 0.0, storage);
1391 DNekMat &Sveft = (*m_vertexedgefacetransformmatrix);
1395 1, efRow, 0.0, storage);
1396 DNekMat &Svef = (*m_vertexedgefacecoupling);
1399 for (n = 0; n < nedgemodesconnected; ++n)
1402 VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid),
1406 Svef.SetValue(0, n, VertexEdgeFaceValue);
1410 for (n = 0; n < nfacemodesconnected; ++n)
1413 VertexEdgeFaceValue = (*r_bnd)(GetVertexMap(vid),
1417 Svef.SetValue(0, n + nedgemodesconnected, VertexEdgeFaceValue);
1430 efRow, efRow, 0.0, storage);
1431 DNekMat &Sefef = (*m_edgefacecoupling);
1436 for (m = 0; m < nedgemodesconnected; ++m)
1438 for (n = 0; n < nedgemodesconnected; ++n)
1441 EdgeEdgeValue = (*r_bnd)(edgemodearray[n],
1445 Sefef.SetValue(n, m, EdgeEdgeValue);
1450 for (n = 0; n < nfacemodesconnected; ++n)
1452 for (m = 0; m < nfacemodesconnected; ++m)
1455 FaceFaceValue = (*r_bnd)(facemodearray[n],
1458 Sefef.SetValue(nedgemodesconnected + n,
1459 nedgemodesconnected + m, FaceFaceValue);
1464 for (n = 0; n < nedgemodesconnected; ++n)
1466 for (m = 0; m < nfacemodesconnected; ++m)
1469 FaceFaceValue = (*r_bnd)(edgemodearray[n],
1474 nedgemodesconnected + m,
1478 Sefef.SetValue(nedgemodesconnected + m,
1491 Sveft = -Svef * Sefef;
1495 for (n = 0; n < edgemodearray.num_elements(); ++n)
1497 RT.SetValue(edgemodearray[n], GetVertexMap(vid),
1499 R.SetValue(GetVertexMap(vid), edgemodearray[n],
1504 for (n = 0; n < facemodearray.num_elements(); ++n)
1506 RT.SetValue(facemodearray[n], GetVertexMap(vid),
1507 Sveft(0, n + nedgemodesconnected));
1508 R.SetValue(GetVertexMap(vid), facemodearray[n],
1509 Sveft(0, n + nedgemodesconnected));
1532 int efCol, efRow, nedgemodes;
1535 nConnectedFaces = 2;
1546 for (eid = 0; eid < nEdges; ++eid)
1549 efCol = GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
1550 GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
1551 efRow = GetEdgeNcoeffs(eid) - 2;
1556 efRow, efCol, 0.0, storage);
1557 DNekMat &Mef = (*m_efedgefacecoupling);
1562 efCol, efCol, 0.0, storage);
1563 DNekMat &Meff = (*m_effacefacecoupling);
1568 efRow, efCol, 0.0, storage);
1569 DNekMat &Meft = (*m_edgefacetransformmatrix);
1571 int nfacemodesconnected =
1572 GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 0)) +
1573 GetFaceIntNcoeffs(geom->GetEdgeFaceMap(eid, 1));
1575 facemodearray(nfacemodesconnected);
1579 = GetEdgeInverseBoundaryMap(eid);
1580 nedgemodes = GetEdgeNcoeffs(eid) - 2;
1586 1, &edgemodearray[0], 1);
1592 for (fid = 0; fid < nConnectedFaces; ++fid)
1594 MatFaceLocation[fid] = GetFaceInverseBoundaryMap(
1595 geom->GetEdgeFaceMap(eid, fid));
1596 nmodes = MatFaceLocation[fid].num_elements();
1601 1, &facemodearray[offset], 1);
1607 for (n = 0; n < nedgemodes; ++n)
1609 for (m = 0; m < nfacemodesconnected; ++m)
1612 EdgeFaceValue = (*r_bnd)(edgemodearray[n],
1616 Mef.SetValue(n, m, EdgeFaceValue);
1621 for (n = 0; n < nfacemodesconnected; ++n)
1623 for (m = 0; m < nfacemodesconnected; ++m)
1626 FaceFaceValue = (*r_bnd)(facemodearray[n],
1630 Meff.SetValue(n, m, FaceFaceValue);
1644 for (n = 0; n < Meft.GetRows(); ++n)
1646 for (m = 0; m < Meft.GetColumns(); ++m)
1648 R.SetValue(edgemodearray[n], facemodearray[m],
1650 RT.SetValue(facemodearray[m], edgemodearray[n],
1656 for (i = 0; i < R.GetRows(); ++i)
1658 R.SetValue(i, i, 1.0);
1659 RT.SetValue(i, i, 1.0);
1665 return m_transformationmatrix;
1670 return m_transposedtransformationmatrix;
1691 int i, j, n, eid = 0, fid = 0;
1692 int nCoeffs = NumBndryCoeffs();
1702 nCoeffs, nCoeffs, 0.0, storage);
1703 DNekMat &InvR = (*m_inversetransformationmatrix);
1705 int nVerts = GetNverts();
1706 int nEdges = GetNedges();
1707 int nFaces = GetNfaces();
1711 int nedgemodestotal = 0;
1712 int nfacemodestotal = 0;
1714 for (eid = 0; eid < nEdges; ++eid)
1716 nedgemodes = GetEdgeNcoeffs(eid) - 2;
1717 nedgemodestotal += nedgemodes;
1720 for (fid = 0; fid < nFaces; ++fid)
1722 nfacemodes = GetFaceIntNcoeffs(fid);
1723 nfacemodestotal += nfacemodes;
1727 edgemodearray(nedgemodestotal);
1729 facemodearray(nfacemodestotal);
1734 for (eid = 0; eid < nEdges; ++eid)
1737 = GetEdgeInverseBoundaryMap(eid);
1738 nedgemodes = GetEdgeNcoeffs(eid) - 2;
1744 &edgemodearray[offset], 1);
1747 offset += nedgemodes;
1753 for (fid = 0; fid < nFaces; ++fid)
1756 = GetFaceInverseBoundaryMap(fid);
1757 nfacemodes = GetFaceIntNcoeffs(fid);
1763 &facemodearray[offset], 1);
1766 offset += nfacemodes;
1770 for (i = 0; i < nVerts; ++i)
1772 for (j = 0; j < nedgemodestotal; ++j)
1775 GetVertexMap(i), edgemodearray[j],
1776 -R(GetVertexMap(i), edgemodearray[j]));
1778 for (j = 0; j < nfacemodestotal; ++j)
1781 GetVertexMap(i), facemodearray[j],
1782 -R(GetVertexMap(i), facemodearray[j]));
1783 for (n = 0; n < nedgemodestotal; ++n)
1785 MatrixValue = InvR.GetValue(GetVertexMap(i),
1787 + R(GetVertexMap(i), edgemodearray[n])
1788 * R(edgemodearray[n], facemodearray[j]);
1789 InvR.SetValue(GetVertexMap(i),
1797 for (i = 0; i < nedgemodestotal; ++i)
1799 for (j = 0; j < nfacemodestotal; ++j)
1802 edgemodearray[i], facemodearray[j],
1803 -R(edgemodearray[i], facemodearray[j]));
1807 for (i = 0; i < nCoeffs; ++i)
1809 InvR.SetValue(i, i, 1.0);
1812 return m_inversetransformationmatrix;
1820 int nBndCoeffs = NumBndryCoeffs();
1823 GetBoundaryMap(bmap);
1827 map<int, int> invmap;
1828 for (j = 0; j < nBndCoeffs; ++j)
1830 invmap[bmap[j]] = j;
1834 nEdgeCoeffs = GetEdgeNcoeffs(eid) - 2;
1840 geom->GetEorient(eid);
1847 GetEdgeInteriorMap(eid, eOrient, maparray, signarray);
1849 for (n = 0; n < nEdgeCoeffs; ++n)
1851 edgemaparray[n] = invmap[maparray[n]];
1854 return edgemaparray;
1858 Expansion3D::v_GetFaceInverseBoundaryMap(
1865 int nBndCoeffs = NumBndryCoeffs();
1868 GetBoundaryMap(bmap);
1872 map<int, int> reversemap;
1873 for (j = 0; j < bmap.num_elements(); ++j)
1875 reversemap[bmap[j]] = j;
1879 nFaceCoeffs = GetFaceIntNcoeffs(fid);
1890 fOrient = GetForient(fid);
1894 fOrient = faceOrient;
1898 GetFaceInteriorMap(fid, fOrient, maparray, signarray);
1900 for (n = 0; n < nFaceCoeffs; ++n)
1902 facemaparray[n] = reversemap[maparray[n]];
1905 return facemaparray;
1910 return m_geom->GetForient(face);
1917 void Expansion3D::v_GetTracePhysVals(
1924 v_GetFacePhysVals(face,FaceExp,inarray,outarray,orient);
1927 void Expansion3D::v_GetFacePhysVals(
1937 orient = GetForient(face);
1940 int nq0 = FaceExp->GetNumPoints(0);
1941 int nq1 = FaceExp->GetNumPoints(1);
1943 int nfacepts = GetFaceNumPoints(face);
1944 int dir0 = GetGeom3D()->GetDir(face,0);
1945 int dir1 = GetGeom3D()->GetDir(face,1);
1952 GetFacePhysMap(face,faceids);
1953 Vmath::Gathr(faceids.num_elements(),inarray,faceids,o_tmp);
1971 m_base[dir1]->GetPointsKey(),
1973 FaceExp->GetBasis(to_id0)->GetPointsKey(),
1974 FaceExp->GetBasis(to_id1)->GetPointsKey(),
1978 ReOrientFacePhysMap(FaceExp->GetNverts(),orient,nq0,nq1,faceids);
1982 void Expansion3D::ReOrientFacePhysMap(
const int nvert,
1984 const int nq0,
const int nq1,
1989 ReOrientTriFacePhysMap(orient,nq0,nq1,idmap);
1993 ReOrientQuadFacePhysMap(orient,nq0,nq1,idmap);
2002 if(idmap.num_elements() != nq0*nq1)
2011 for(
int i = 0; i < nq0*nq1; ++i)
2018 for (
int j = 0; j < nq1; ++j)
2020 for(
int i = 0; i < nq0; ++i)
2022 idmap[j*nq0+i] = nq0-1-i +j*nq0;
2027 ASSERTL0(
false,
"Case not supposed to be used in this function");
2037 if(idmap.num_elements() != nq0*nq1)
2047 for(
int i = 0; i < nq0*nq1; ++i)
2055 for (
int j = 0; j < nq1; j++)
2057 for (
int i =0; i < nq0; ++i)
2059 idmap[j*nq0+i] = nq0-1-i + j*nq0;
2067 for (
int j = 0; j < nq1; j++)
2069 for (
int i =0; i < nq0; ++i)
2071 idmap[j*nq0+i] = nq0*(nq1-1-j)+i;
2078 for (
int j = 0; j < nq1; j++)
2080 for (
int i =0; i < nq0; ++i)
2082 idmap[j*nq0+i] = nq0*nq1-1-j*nq0 -i;
2089 for (
int i =0; i < nq0; ++i)
2091 for (
int j = 0; j < nq1; ++j)
2093 idmap[i*nq1+j] = i + j*nq0;
2101 for (
int i =0; i < nq0; ++i)
2103 for (
int j = 0; j < nq1; ++j)
2105 idmap[i*nq1+j] = nq0-1-i + j*nq0;
2113 for (
int i =0; i < nq0; ++i)
2115 for (
int j = 0; j < nq1; ++j)
2117 idmap[i*nq1+j] = i+nq0*(nq1-1)-j*nq0;
2125 for (
int i =0; i < nq0; ++i)
2127 for (
int j = 0; j < nq1; ++j)
2129 idmap[i*nq1+j] = nq0*nq1-1-i-j*nq0;
2135 ASSERTL0(
false,
"Unknow orientation");
2140 void Expansion3D::v_NormVectorIProductWRTBase(
2144 NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
NekDouble GetConstFactor(const ConstFactorType &factor) const
#define ASSERTL0(condition, msg)
const ConstFactorMap & GetConstFactors() const
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
const VarCoeffMap & GetVarCoeffs() const
boost::shared_ptr< IndexMapValues > IndexMapValuesSharedPtr
MatrixType GetMatrixType() const
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
bool HasVarCoeff(const StdRegions::VarCoeffType &coeff) const
#define sign(a, b)
return the sign(b)*a
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
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 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
boost::shared_ptr< DNekMat > DNekMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
static DNekMatSharedPtr NullDNekMatSharedPtr
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::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< Expansion > ExpansionSharedPtr
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
static ConstFactorMap NullConstFactorMap
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr