85 WARNINGL0(
false,
"LFR is deprecated, use LDG instead");
93 int nConvectiveFields = pFields.size();
94 int nDim = pFields[0]->GetCoordim(0);
95 int nSolutionPts = pFields[0]->GetTotPoints();
96 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
116 for (i = 0; i < nConvectiveFields; ++i)
132 for (j = 0; j < nDim; ++j)
171 int nLocalSolutionPts;
172 int nElements = pFields[0]->GetExpSize();
173 int nDimensions = pFields[0]->GetCoordim(0);
174 int nSolutionPts = pFields[0]->GetTotPoints();
175 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
178 for (i = 0; i < nDimensions; ++i)
197 for (n = 0; n < nElements; ++n)
199 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
200 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
201 phys_offset = pFields[0]->GetPhys_Offset(n);
208 for (i = 0; i < nLocalSolutionPts; ++i)
210 m_jac[i + phys_offset] = jac[0];
228 for (n = 0; n < nElements; ++n)
230 base = pFields[0]->GetExp(n)->GetBase();
231 nquad0 = base[0]->GetNumPoints();
232 nquad1 = base[1]->GetNumPoints();
240 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
242 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
244 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
246 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
249 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
250 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
251 phys_offset = pFields[0]->GetPhys_Offset(n);
264 ->GetDerivFactors(ptsKeys);
273 for (i = 0; i < nLocalSolutionPts; ++i)
275 m_jac[i + phys_offset] = jac[i];
276 m_gmat[0][i + phys_offset] = gmat[0][i];
277 m_gmat[1][i + phys_offset] = gmat[1][i];
278 m_gmat[2][i + phys_offset] = gmat[2][i];
279 m_gmat[3][i + phys_offset] = gmat[3][i];
284 for (i = 0; i < nLocalSolutionPts; ++i)
286 m_jac[i + phys_offset] = jac[0];
287 m_gmat[0][i + phys_offset] = gmat[0][0];
288 m_gmat[1][i + phys_offset] = gmat[1][0];
289 m_gmat[2][i + phys_offset] = gmat[2][0];
290 m_gmat[3][i + phys_offset] = gmat[3][0];
298 ASSERTL0(
false,
"3DFR Metric terms not implemented yet");
303 ASSERTL0(
false,
"Expansion dimension not recognised");
334 int nquad0, nquad1, nquad2;
335 int nmodes0, nmodes1, nmodes2;
338 int nElements = pFields[0]->GetExpSize();
339 int nDim = pFields[0]->GetCoordim(0);
348 for (n = 0; n < nElements; ++n)
350 base = pFields[0]->GetExp(n)->GetBase();
351 nquad0 = base[0]->GetNumPoints();
352 nmodes0 = base[0]->GetNumModes();
356 base[0]->GetZW(z0, w0);
368 int p0 = nmodes0 - 1;
375 std::tgamma(2 * p0 + 1) /
376 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
386 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
387 (ap0 * std::tgamma(p0 + 1)) *
388 (ap0 * std::tgamma(p0 + 1)));
392 c0 = 2.0 * (p0 + 1.0) /
393 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
394 (ap0 * std::tgamma(p0 + 1)));
399 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
400 (ap0 * std::tgamma(p0 + 1)));
404 c0 = 10000000000000000.0;
407 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
408 (ap0 * std::tgamma(p0 + 1)) *
409 (ap0 * std::tgamma(p0 + 1));
411 NekDouble overeta0 = 1.0 / (1.0 + etap0);
424 for (i = 0; i < nquad0; ++i)
434 for (i = 0; i < nquad0; ++i)
457 for (n = 0; n < nElements; ++n)
459 base = pFields[0]->GetExp(n)->GetBase();
460 nquad0 = base[0]->GetNumPoints();
461 nquad1 = base[1]->GetNumPoints();
462 nmodes0 = base[0]->GetNumModes();
463 nmodes1 = base[1]->GetNumModes();
470 base[0]->GetZW(z0, w0);
471 base[1]->GetZW(z1, w1);
488 int p0 = nmodes0 - 1;
489 int p1 = nmodes1 - 1;
497 std::tgamma(2 * p0 + 1) /
498 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
501 std::tgamma(2 * p1 + 1) /
502 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
513 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
514 (ap0 * std::tgamma(p0 + 1)) *
515 (ap0 * std::tgamma(p0 + 1)));
518 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
519 (ap1 * std::tgamma(p1 + 1)) *
520 (ap1 * std::tgamma(p1 + 1)));
524 c0 = 2.0 * (p0 + 1.0) /
525 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
526 (ap0 * std::tgamma(p0 + 1)));
528 c1 = 2.0 * (p1 + 1.0) /
529 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
530 (ap1 * std::tgamma(p1 + 1)));
535 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
536 (ap0 * std::tgamma(p0 + 1)));
539 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
540 (ap1 * std::tgamma(p1 + 1)));
544 c0 = 10000000000000000.0;
545 c1 = 10000000000000000.0;
548 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
549 (ap0 * std::tgamma(p0 + 1)) *
550 (ap0 * std::tgamma(p0 + 1));
552 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
553 (ap1 * std::tgamma(p1 + 1)) *
554 (ap1 * std::tgamma(p1 + 1));
556 NekDouble overeta0 = 1.0 / (1.0 + etap0);
557 NekDouble overeta1 = 1.0 / (1.0 + etap1);
576 for (i = 0; i < nquad0; ++i)
586 for (i = 0; i < nquad1; ++i)
596 for (i = 0; i < nquad0; ++i)
606 for (i = 0; i < nquad1; ++i)
626 for (n = 0; n < nElements; ++n)
628 base = pFields[0]->GetExp(n)->GetBase();
629 nquad0 = base[0]->GetNumPoints();
630 nquad1 = base[1]->GetNumPoints();
631 nquad2 = base[2]->GetNumPoints();
632 nmodes0 = base[0]->GetNumModes();
633 nmodes1 = base[1]->GetNumModes();
634 nmodes2 = base[2]->GetNumModes();
643 base[0]->GetZW(z0, w0);
644 base[1]->GetZW(z1, w1);
645 base[1]->GetZW(z2, w2);
667 int p0 = nmodes0 - 1;
668 int p1 = nmodes1 - 1;
669 int p2 = nmodes2 - 1;
678 std::tgamma(2 * p0 + 1) /
679 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
683 std::tgamma(2 * p1 + 1) /
684 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
688 std::tgamma(2 * p2 + 1) /
689 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
701 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
702 (ap0 * std::tgamma(p0 + 1)) *
703 (ap0 * std::tgamma(p0 + 1)));
706 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
707 (ap1 * std::tgamma(p1 + 1)) *
708 (ap1 * std::tgamma(p1 + 1)));
711 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
712 (ap2 * std::tgamma(p2 + 1)) *
713 (ap2 * std::tgamma(p2 + 1)));
717 c0 = 2.0 * (p0 + 1.0) /
718 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
719 (ap0 * std::tgamma(p0 + 1)));
721 c1 = 2.0 * (p1 + 1.0) /
722 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
723 (ap1 * std::tgamma(p1 + 1)));
725 c2 = 2.0 * (p2 + 1.0) /
726 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
727 (ap2 * std::tgamma(p2 + 1)));
732 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
733 (ap0 * std::tgamma(p0 + 1)));
736 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
737 (ap1 * std::tgamma(p1 + 1)));
740 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
741 (ap2 * std::tgamma(p2 + 1)));
745 c0 = 10000000000000000.0;
746 c1 = 10000000000000000.0;
747 c2 = 10000000000000000.0;
750 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
751 (ap0 * std::tgamma(p0 + 1)) *
752 (ap0 * std::tgamma(p0 + 1));
754 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
755 (ap1 * std::tgamma(p1 + 1)) *
756 (ap1 * std::tgamma(p1 + 1));
758 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
759 (ap2 * std::tgamma(p2 + 1)) *
760 (ap2 * std::tgamma(p2 + 1));
762 NekDouble overeta0 = 1.0 / (1.0 + etap0);
763 NekDouble overeta1 = 1.0 / (1.0 + etap1);
764 NekDouble overeta2 = 1.0 / (1.0 + etap2);
789 for (i = 0; i < nquad0; ++i)
799 for (i = 0; i < nquad1; ++i)
809 for (i = 0; i < nquad2; ++i)
819 for (i = 0; i < nquad0; ++i)
829 for (i = 0; i < nquad1; ++i)
839 for (i = 0; i < nquad2; ++i)
852 ASSERTL0(
false,
"Expansion dimension not recognised");
864 const std::size_t nConvectiveFields,
876 Basis = fields[0]->GetExp(0)->GetBase();
878 int nElements = fields[0]->GetExpSize();
879 int nDim = fields[0]->GetCoordim(0);
880 int nSolutionPts = fields[0]->GetTotPoints();
881 int nCoeffs = fields[0]->GetNcoeffs();
884 for (i = 0; i < nConvectiveFields; ++i)
897 for (i = 0; i < nConvectiveFields; ++i)
901 for (n = 0; n < nElements; n++)
903 phys_offset = fields[0]->GetPhys_Offset(n);
905 fields[i]->GetExp(n)->PhysDeriv(
906 0, auxArray1 = inarray[i] + phys_offset,
907 auxArray2 =
m_DU1[i][0] + phys_offset);
923 1, &
m_D1[i][0][0], 1);
933 for (i = 0; i < nConvectiveFields; ++i)
937 for (n = 0; n < nElements; n++)
939 phys_offset = fields[0]->GetPhys_Offset(n);
941 fields[i]->GetExp(n)->PhysDeriv(
942 0, auxArray1 =
m_D1[i][0] + phys_offset,
943 auxArray2 =
m_DD1[i][0] + phys_offset);
959 1, &outarray[i][0], 1);
962 if (!(
Basis[0]->Collocation()))
964 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
965 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
973 for (i = 0; i < nConvectiveFields; ++i)
975 for (j = 0; j < nDim; ++j)
984 &
m_gmat[0][0], 1, &u1_hat[0], 1);
990 &
m_gmat[1][0], 1, &u2_hat[0], 1);
998 &
m_gmat[2][0], 1, &u1_hat[0], 1);
1004 &
m_gmat[3][0], 1, &u2_hat[0], 1);
1010 for (n = 0; n < nElements; n++)
1012 phys_offset = fields[0]->GetPhys_Offset(n);
1014 fields[i]->GetExp(n)->StdPhysDeriv(
1015 0, auxArray1 = u1_hat + phys_offset,
1016 auxArray2 =
m_tmp1[i][j] + phys_offset);
1018 fields[i]->GetExp(n)->StdPhysDeriv(
1019 1, auxArray1 = u2_hat + phys_offset,
1020 auxArray2 =
m_tmp2[i][j] + phys_offset);
1028 &
m_DU1[i][j][0], 1);
1032 DerCFlux_2D(nConvectiveFields, j, fields, inarray[i],
1038 for (j = 0; j < nSolutionPts; ++j)
1059 for (j = 0; j < nSolutionPts; j++)
1071 for (j = 0; j < nDim; ++j)
1084 for (i = 0; i < nConvectiveFields; ++i)
1089 for (j = 0; j < nSolutionPts; j++)
1100 for (n = 0; n < nElements; n++)
1102 phys_offset = fields[0]->GetPhys_Offset(n);
1104 fields[0]->GetExp(n)->StdPhysDeriv(
1105 0, auxArray1 = f_hat + phys_offset,
1106 auxArray2 =
m_DD1[i][0] + phys_offset);
1108 fields[0]->GetExp(n)->StdPhysDeriv(
1109 1, auxArray1 = g_hat + phys_offset,
1110 auxArray2 =
m_DD1[i][1] + phys_offset);
1118 if (
Basis[0]->GetPointsType() ==
1120 Basis[1]->GetPointsType() ==
1134 &outarray[i][0], 1);
1139 &outarray[i][0], 1);
1142 if (!(
Basis[0]->Collocation()))
1144 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1145 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1153 ASSERTL0(
false,
"3D FRDG case not implemented yet");
1169 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1170 int nvariables = fields.size();
1171 int nDim = fields[0]->GetCoordim(0);
1179 for (i = 0; i < nDim; ++i)
1187 for (j = 0; j < nDim; ++j)
1189 for (i = 0; i < nvariables; ++i)
1192 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1202 fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1213 if (fields[0]->GetBndCondExpansions().size())
1244 int nBndEdgePts, nBndEdges;
1246 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1247 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1251 fields[var]->ExtractTracePhys(ufield, uplus);
1254 for (i = 0; i < nBndRegions; ++i)
1257 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1260 for (e = 0; e < nBndEdges; ++e)
1263 nBndEdgePts = fields[var]
1264 ->GetBndCondExpansions()[i]
1269 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1272 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1273 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1277 ->GetBndConditions()[i]
1282 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1283 1, &penaltyflux[id2], 1);
1286 else if ((fields[var]->GetBndConditions()[i])
1287 ->GetBoundaryConditionType() ==
1290 Vmath::Vcopy(nBndEdgePts, &uplus[id2], 1, &penaltyflux[id2], 1);
1307 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1308 int nvariables = fields.size();
1309 int nDim = fields[0]->GetCoordim(0);
1322 for (i = 0; i < nDim; ++i)
1329 for (i = 0; i < nvariables; ++i)
1332 for (j = 0; j < nDim; ++j)
1335 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1349 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1355 if (fields[0]->GetBndCondExpansions().size())
1364 Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1, qflux[i], 1);
1379 int nBndEdges, nBndEdgePts;
1380 int nBndRegions = fields[var]->GetBndCondExpansions().size();
1381 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1387 fields[var]->ExtractTracePhys(qfield, qtemp);
1389 for (i = 0; i < nBndRegions; ++i)
1391 nBndEdges = fields[var]->GetBndCondExpansions()[i]->GetExpSize();
1394 for (e = 0; e < nBndEdges; ++e)
1396 nBndEdgePts = fields[var]
1397 ->GetBndCondExpansions()[i]
1401 id1 = fields[var]->GetBndCondExpansions()[i]->GetPhys_Offset(e);
1403 id2 = fields[0]->GetTrace()->GetPhys_Offset(
1404 fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt++));
1409 ->GetBndConditions()[i]
1413 &qtemp[id2], 1, &penaltyflux[id2], 1);
1416 else if ((fields[var]->GetBndConditions()[i])
1417 ->GetBoundaryConditionType() ==
1422 &(fields[var]->GetBndCondExpansions()[i]->GetPhys())[id1],
1423 1, &penaltyflux[id2], 1);
1440 [[maybe_unused]]
const int nConvectiveFields,
1446 int nLocalSolutionPts, phys_offset, t_offset;
1454 Basis = fields[0]->GetExp(0)->GetBasis(0);
1456 int nElements = fields[0]->GetExpSize();
1457 int nSolutionPts = fields[0]->GetTotPoints();
1459 vector<bool> leftAdjacentTraces =
1460 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1461 ->GetLeftAdjacentTraces();
1472 fields[0]->GetTraceMap()->GetElmtToTrace();
1474 for (n = 0; n < nElements; ++n)
1476 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1477 phys_offset = fields[0]->GetPhys_Offset(n);
1481 Vmath::Vcopy(nLocalSolutionPts, &flux[phys_offset], 1, &tmparrayX1[0],
1484 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1487 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1488 elmtToTrace[n][0]->GetElmtId());
1490 if (leftAdjacentTraces[2 * n])
1492 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1496 JumpL[n] = iFlux[t_offset] - tmpFluxVertex[0];
1499 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1500 elmtToTrace[n][1]->GetElmtId());
1502 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1504 if (leftAdjacentTraces[2 * n + 1])
1506 JumpR[n] = iFlux[t_offset] - tmpFluxVertex[0];
1510 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex[0];
1514 for (n = 0; n < nElements; ++n)
1516 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1517 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1518 phys_offset = fields[0]->GetPhys_Offset(n);
1526 JumpL[n] = JumpL[n] * jac[0];
1527 JumpR[n] = JumpR[n] * jac[0];
1538 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1539 &derCFlux[phys_offset], 1);
1559 [[maybe_unused]]
const int nConvectiveFields,
const int direction,
1564 int n, e, i, j, cnt;
1568 int nElements = fields[0]->GetExpSize();
1569 int trace_offset, phys_offset;
1570 int nLocalSolutionPts;
1578 fields[0]->GetTraceMap()->GetElmtToTrace();
1581 for (n = 0; n < nElements; ++n)
1584 phys_offset = fields[0]->GetPhys_Offset(n);
1585 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1586 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1595 base = fields[0]->GetExp(n)->GetBase();
1596 nquad0 = base[0]->GetNumPoints();
1597 nquad1 = base[1]->GetNumPoints();
1605 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1608 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1614 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1615 elmtToTrace[n][e]->GetElmtId());
1620 fields[0]->GetExp(n)->GetTracePhysVals(
1621 e, elmtToTrace[n][e], flux + phys_offset, auxArray1 = tmparray);
1626 Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1, &tmparray[0], 1,
1631 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1640 ->as<LocalRegions::Expansion2D>()
1648 fields[0]->GetExp(n)->GetTracePhysVals(
1649 e, elmtToTrace[n][e], jac, auxArray1 = jacEdge);
1653 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1661 for (j = 0; j < nEdgePts; j++)
1663 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1671 Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1, fluxJumps, 1);
1680 for (i = 0; i < nquad0; ++i)
1682 for (j = 0; j < nquad1; ++j)
1684 cnt = i + j * nquad0;
1685 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
1690 for (i = 0; i < nquad1; ++i)
1692 for (j = 0; j < nquad0; ++j)
1694 cnt = (nquad0)*i + j;
1695 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1700 for (i = 0; i < nquad0; ++i)
1702 for (j = 0; j < nquad1; ++j)
1704 cnt = j * nquad0 + i;
1705 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1710 for (i = 0; i < nquad1; ++i)
1712 for (j = 0; j < nquad0; ++j)
1714 cnt = j + i * nquad0;
1715 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
1721 ASSERTL0(
false,
"edge value (< 3) is out of range");
1730 Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1, &divCFluxE3[0], 1,
1731 &derCFlux[phys_offset], 1);
1733 else if (direction == 1)
1735 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE2[0], 1,
1736 &derCFlux[phys_offset], 1);
1755 [[maybe_unused]]
const int nConvectiveFields,
1762 int n, e, i, j, cnt;
1764 int nElements = fields[0]->GetExpSize();
1766 int nLocalSolutionPts;
1777 fields[0]->GetTraceMap()->GetElmtToTrace();
1780 for (n = 0; n < nElements; ++n)
1783 phys_offset = fields[0]->GetPhys_Offset(n);
1784 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1786 base = fields[0]->GetExp(n)->GetBase();
1787 nquad0 = base[0]->GetNumPoints();
1788 nquad1 = base[1]->GetNumPoints();
1796 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1799 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1808 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1809 elmtToTrace[n][e]->GetElmtId());
1813 fields[0]->GetExp(n)->GetTraceNormal(e);
1817 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1818 fluxX1 + phys_offset,
1819 auxArray1 = tmparrayX1);
1823 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1824 fluxX2 + phys_offset,
1825 auxArray1 = tmparrayX2);
1828 for (i = 0; i < nEdgePts; ++i)
1835 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1839 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1843 auxArray2 = fluxJumps, 1);
1846 for (i = 0; i < nEdgePts; ++i)
1851 fluxJumps[i] = -fluxJumps[i];
1859 for (i = 0; i < nquad0; ++i)
1862 fluxJumps[i] = -(
m_Q2D_e0[n][i]) * fluxJumps[i];
1864 for (j = 0; j < nquad1; ++j)
1866 cnt = i + j * nquad0;
1867 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
1872 for (i = 0; i < nquad1; ++i)
1875 fluxJumps[i] = (
m_Q2D_e1[n][i]) * fluxJumps[i];
1877 for (j = 0; j < nquad0; ++j)
1879 cnt = (nquad0)*i + j;
1880 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1885 for (i = 0; i < nquad0; ++i)
1888 fluxJumps[i] = (
m_Q2D_e2[n][i]) * fluxJumps[i];
1890 for (j = 0; j < nquad1; ++j)
1892 cnt = j * nquad0 + i;
1893 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1898 for (i = 0; i < nquad1; ++i)
1901 fluxJumps[i] = -(
m_Q2D_e3[n][i]) * fluxJumps[i];
1902 for (j = 0; j < nquad0; ++j)
1904 cnt = j + i * nquad0;
1905 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
1911 ASSERTL0(
false,
"edge value (< 3) is out of range");
1917 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1918 &divCFlux[phys_offset], 1);
1920 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1921 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1923 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1924 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1943 [[maybe_unused]]
const int nConvectiveFields,
1950 int n, e, i, j, cnt;
1952 int nElements = fields[0]->GetExpSize();
1953 int nLocalSolutionPts;
1964 fields[0]->GetTraceMap()->GetElmtToTrace();
1967 for (n = 0; n < nElements; ++n)
1970 phys_offset = fields[0]->GetPhys_Offset(n);
1971 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1973 base = fields[0]->GetExp(n)->GetBase();
1974 nquad0 = base[0]->GetNumPoints();
1975 nquad1 = base[1]->GetNumPoints();
1983 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1986 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1995 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1996 elmtToTrace[n][e]->GetElmtId());
2000 fields[0]->GetExp(n)->GetTraceNormal(e);
2009 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2010 fluxX2 + phys_offset,
2011 auxArray1 = fluxN_D);
2016 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2020 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2024 auxArray2 = fluxN, 1);
2027 auxArray2 = fluxN_D, 1);
2031 for (i = 0; i < nquad0; ++i)
2034 fluxN_R[i] = (
m_Q2D_e0[n][i]) * fluxN[i];
2037 for (i = 0; i < nEdgePts; ++i)
2044 fluxN_R[i] = -fluxN_R[i];
2050 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2055 for (i = 0; i < nquad0; ++i)
2057 for (j = 0; j < nquad1; ++j)
2059 cnt = i + j * nquad0;
2060 divCFluxE0[cnt] = -fluxJumps[i] *
m_dGL_xi2[n][j];
2068 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2069 fluxX1 + phys_offset,
2070 auxArray1 = fluxN_D);
2073 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2077 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2081 auxArray2 = fluxN, 1);
2084 auxArray2 = fluxN_D, 1);
2088 for (i = 0; i < nquad1; ++i)
2091 fluxN_R[i] = (
m_Q2D_e1[n][i]) * fluxN[i];
2094 for (i = 0; i < nEdgePts; ++i)
2101 fluxN_R[i] = -fluxN_R[i];
2107 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2112 for (i = 0; i < nquad1; ++i)
2114 for (j = 0; j < nquad0; ++j)
2116 cnt = (nquad0)*i + j;
2117 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
2127 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2128 fluxX2 + phys_offset,
2129 auxArray1 = fluxN_D);
2132 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2136 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2140 auxArray2 = fluxN, 1);
2143 auxArray2 = fluxN_D, 1);
2147 for (i = 0; i < nquad0; ++i)
2150 fluxN_R[i] = (
m_Q2D_e2[n][i]) * fluxN[i];
2153 for (i = 0; i < nEdgePts; ++i)
2160 fluxN_R[i] = -fluxN_R[i];
2167 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2172 for (i = 0; i < nquad0; ++i)
2174 for (j = 0; j < nquad1; ++j)
2176 cnt = j * nquad0 + i;
2177 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
2186 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
2187 fluxX1 + phys_offset,
2188 auxArray1 = fluxN_D);
2192 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
2196 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
2200 auxArray2 = fluxN, 1);
2203 auxArray2 = fluxN_D, 1);
2207 for (i = 0; i < nquad1; ++i)
2210 fluxN_R[i] = (
m_Q2D_e3[n][i]) * fluxN[i];
2213 for (i = 0; i < nEdgePts; ++i)
2220 fluxN_R[i] = -fluxN_R[i];
2227 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
2232 for (i = 0; i < nquad1; ++i)
2234 for (j = 0; j < nquad0; ++j)
2236 cnt = j + i * nquad0;
2237 divCFluxE3[cnt] = -fluxJumps[i] *
m_dGL_xi1[n][j];
2242 ASSERTL0(
false,
"edge value (< 3) is out of range");
2248 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
2249 &divCFlux[phys_offset], 1);
2251 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2252 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2254 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2255 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
#define ASSERTL0(condition, msg)
#define WARNINGL0(condition, msg)
Represents a basis of a given type.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
const SpatialDomains::GeomFactorsSharedPtr & GetMetricInfo() const
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
void v_Diffuse(const size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override
Calculate FR Diffusion for the linear problems using an LDG interface flux.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Array< OneD, NekDouble > m_jac
void DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term....
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_divFC
void DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
void SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
LibUtilities::SessionReaderSharedPtr m_session
void WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
static std::string type[]
void WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
void SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72...
void DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
void NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
static DiffusionSharedPtr create(std::string diffType)
Array< OneD, Array< OneD, NekDouble > > m_gmat
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
void NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase DiffusionLFR objects and store them before starting the time-stepping.
Array< OneD, Array< OneD, NekDouble > > m_IF2
Array< OneD, Array< OneD, NekDouble > > m_divFD
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
void DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
DiffusionFactory & GetDiffusionFactory()
@ eDeformed
Geometry is curved or has non-constant factors.
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
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 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 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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vdiv(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 Reverse(int n, const T *x, const int incx, T *y, const int incy)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.