88 WARNINGL0(
false,
"FR is deprecated, use DG instead");
118 int nLocalSolutionPts;
119 int nElements = pFields[0]->GetExpSize();
120 int nDimensions = pFields[0]->GetCoordim(0);
121 int nSolutionPts = pFields[0]->GetTotPoints();
122 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
137 for (n = 0; n < nElements; ++n)
139 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
140 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
141 phys_offset = pFields[0]->GetPhys_Offset(n);
148 for (i = 0; i < nLocalSolutionPts; ++i)
150 m_jac[i + phys_offset] = jac[0];
170 for (n = 0; n < nElements; ++n)
172 base = pFields[0]->GetExp(n)->GetBase();
173 nquad0 = base[0]->GetNumPoints();
174 nquad1 = base[1]->GetNumPoints();
182 pFields[0]->GetExp(n)->GetTraceQFactors(0, auxArray1 =
184 pFields[0]->GetExp(n)->GetTraceQFactors(1, auxArray1 =
186 pFields[0]->GetExp(n)->GetTraceQFactors(2, auxArray1 =
188 pFields[0]->GetExp(n)->GetTraceQFactors(3, auxArray1 =
191 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
192 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
193 phys_offset = pFields[0]->GetPhys_Offset(n);
206 ->GetDerivFactors(ptsKeys);
215 for (i = 0; i < nLocalSolutionPts; ++i)
217 m_jac[i + phys_offset] = jac[i];
218 m_gmat[0][i + phys_offset] = gmat[0][i];
219 m_gmat[1][i + phys_offset] = gmat[1][i];
220 m_gmat[2][i + phys_offset] = gmat[2][i];
221 m_gmat[3][i + phys_offset] = gmat[3][i];
226 for (i = 0; i < nLocalSolutionPts; ++i)
228 m_jac[i + phys_offset] = jac[0];
229 m_gmat[0][i + phys_offset] = gmat[0][0];
230 m_gmat[1][i + phys_offset] = gmat[1][0];
231 m_gmat[2][i + phys_offset] = gmat[2][0];
232 m_gmat[3][i + phys_offset] = gmat[3][0];
238 for (i = 0; i < nDimensions; ++i)
248 ASSERTL0(
false,
"3DFR Metric terms not implemented yet");
253 ASSERTL0(
false,
"Expansion dimension not recognised");
283 int nquad0, nquad1, nquad2;
284 int nmodes0, nmodes1, nmodes2;
287 int nElements = pFields[0]->GetExpSize();
288 int nDimensions = pFields[0]->GetCoordim(0);
297 for (n = 0; n < nElements; ++n)
299 base = pFields[0]->GetExp(n)->GetBase();
300 nquad0 = base[0]->GetNumPoints();
301 nmodes0 = base[0]->GetNumModes();
305 base[0]->GetZW(z0, w0);
317 int p0 = nmodes0 - 1;
324 std::tgamma(2 * p0 + 1) /
325 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
335 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
336 (ap0 * std::tgamma(p0 + 1)) *
337 (ap0 * std::tgamma(p0 + 1)));
341 c0 = 2.0 * (p0 + 1.0) /
342 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
343 (ap0 * std::tgamma(p0 + 1)));
348 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
349 (ap0 * std::tgamma(p0 + 1)));
353 c0 = 10000000000000000.0;
356 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
357 (ap0 * std::tgamma(p0 + 1)) *
358 (ap0 * std::tgamma(p0 + 1));
360 NekDouble overeta0 = 1.0 / (1.0 + etap0);
373 for (i = 0; i < nquad0; ++i)
383 for (i = 0; i < nquad0; ++i)
406 for (n = 0; n < nElements; ++n)
408 base = pFields[0]->GetExp(n)->GetBase();
409 nquad0 = base[0]->GetNumPoints();
410 nquad1 = base[1]->GetNumPoints();
411 nmodes0 = base[0]->GetNumModes();
412 nmodes1 = base[1]->GetNumModes();
419 base[0]->GetZW(z0, w0);
420 base[1]->GetZW(z1, w1);
437 int p0 = nmodes0 - 1;
438 int p1 = nmodes1 - 1;
446 std::tgamma(2 * p0 + 1) /
447 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
450 std::tgamma(2 * p1 + 1) /
451 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
462 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
463 (ap0 * std::tgamma(p0 + 1)) *
464 (ap0 * std::tgamma(p0 + 1)));
467 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
468 (ap1 * std::tgamma(p1 + 1)) *
469 (ap1 * std::tgamma(p1 + 1)));
473 c0 = 2.0 * (p0 + 1.0) /
474 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
475 (ap0 * std::tgamma(p0 + 1)));
477 c1 = 2.0 * (p1 + 1.0) /
478 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
479 (ap1 * std::tgamma(p1 + 1)));
484 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
485 (ap0 * std::tgamma(p0 + 1)));
488 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
489 (ap1 * std::tgamma(p1 + 1)));
493 c0 = 10000000000000000.0;
494 c1 = 10000000000000000.0;
497 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
498 (ap0 * std::tgamma(p0 + 1)) *
499 (ap0 * std::tgamma(p0 + 1));
501 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
502 (ap1 * std::tgamma(p1 + 1)) *
503 (ap1 * std::tgamma(p1 + 1));
505 NekDouble overeta0 = 1.0 / (1.0 + etap0);
506 NekDouble overeta1 = 1.0 / (1.0 + etap1);
525 for (i = 0; i < nquad0; ++i)
535 for (i = 0; i < nquad1; ++i)
545 for (i = 0; i < nquad0; ++i)
555 for (i = 0; i < nquad1; ++i)
575 for (n = 0; n < nElements; ++n)
577 base = pFields[0]->GetExp(n)->GetBase();
578 nquad0 = base[0]->GetNumPoints();
579 nquad1 = base[1]->GetNumPoints();
580 nquad2 = base[2]->GetNumPoints();
581 nmodes0 = base[0]->GetNumModes();
582 nmodes1 = base[1]->GetNumModes();
583 nmodes2 = base[2]->GetNumModes();
592 base[0]->GetZW(z0, w0);
593 base[1]->GetZW(z1, w1);
594 base[1]->GetZW(z2, w2);
616 int p0 = nmodes0 - 1;
617 int p1 = nmodes1 - 1;
618 int p2 = nmodes2 - 1;
626 std::tgamma(2 * p0 + 1) /
627 (pow(2.0, p0) * std::tgamma(p0 + 1) * std::tgamma(p0 + 1));
631 std::tgamma(2 * p1 + 1) /
632 (pow(2.0, p1) * std::tgamma(p1 + 1) * std::tgamma(p1 + 1));
636 std::tgamma(2 * p2 + 1) /
637 (pow(2.0, p2) * std::tgamma(p2 + 1) * std::tgamma(p2 + 1));
649 ((2.0 * p0 + 1.0) * (p0 + 1.0) *
650 (ap0 * std::tgamma(p0 + 1)) *
651 (ap0 * std::tgamma(p0 + 1)));
654 ((2.0 * p1 + 1.0) * (p1 + 1.0) *
655 (ap1 * std::tgamma(p1 + 1)) *
656 (ap1 * std::tgamma(p1 + 1)));
659 ((2.0 * p2 + 1.0) * (p2 + 1.0) *
660 (ap2 * std::tgamma(p2 + 1)) *
661 (ap2 * std::tgamma(p2 + 1)));
665 c0 = 2.0 * (p0 + 1.0) /
666 ((2.0 * p0 + 1.0) * p0 * (ap0 * std::tgamma(p0 + 1)) *
667 (ap0 * std::tgamma(p0 + 1)));
669 c1 = 2.0 * (p1 + 1.0) /
670 ((2.0 * p1 + 1.0) * p1 * (ap1 * std::tgamma(p1 + 1)) *
671 (ap1 * std::tgamma(p1 + 1)));
673 c2 = 2.0 * (p2 + 1.0) /
674 ((2.0 * p2 + 1.0) * p2 * (ap2 * std::tgamma(p2 + 1)) *
675 (ap2 * std::tgamma(p2 + 1)));
680 -2.0 / ((2.0 * p0 + 1.0) * (ap0 * std::tgamma(p0 + 1)) *
681 (ap0 * std::tgamma(p0 + 1)));
684 -2.0 / ((2.0 * p1 + 1.0) * (ap1 * std::tgamma(p1 + 1)) *
685 (ap1 * std::tgamma(p1 + 1)));
688 -2.0 / ((2.0 * p2 + 1.0) * (ap2 * std::tgamma(p2 + 1)) *
689 (ap2 * std::tgamma(p2 + 1)));
693 c0 = 10000000000000000.0;
694 c1 = 10000000000000000.0;
695 c2 = 10000000000000000.0;
698 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0) *
699 (ap0 * std::tgamma(p0 + 1)) *
700 (ap0 * std::tgamma(p0 + 1));
702 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0) *
703 (ap1 * std::tgamma(p1 + 1)) *
704 (ap1 * std::tgamma(p1 + 1));
706 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0) *
707 (ap2 * std::tgamma(p2 + 1)) *
708 (ap2 * std::tgamma(p2 + 1));
710 NekDouble overeta0 = 1.0 / (1.0 + etap0);
711 NekDouble overeta1 = 1.0 / (1.0 + etap1);
712 NekDouble overeta2 = 1.0 / (1.0 + etap2);
737 for (i = 0; i < nquad0; ++i)
747 for (i = 0; i < nquad1; ++i)
757 for (i = 0; i < nquad2; ++i)
767 for (i = 0; i < nquad0; ++i)
777 for (i = 0; i < nquad1; ++i)
787 for (i = 0; i < nquad2; ++i)
800 ASSERTL0(
false,
"Expansion dimension not recognised");
819 const int nConvectiveFields,
834 Basis = fields[0]->GetExp(0)->GetBase();
836 int nElements = fields[0]->GetExpSize();
837 int nDimensions = fields[0]->GetCoordim(0);
838 int nSolutionPts = fields[0]->GetTotPoints();
839 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
840 int nCoeffs = fields[0]->GetNcoeffs();
853 for (i = 0; i < nConvectiveFields; ++i)
863 for (i = 0; i < nConvectiveFields; ++i)
869 fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
888 for (i = 0; i < nConvectiveFields; ++i)
890 for (n = 0; n < nElements; n++)
892 phys_offset = fields[0]->GetPhys_Offset(n);
894 fields[i]->GetExp(n)->PhysDeriv(
895 0, fluxvector[i][0] + phys_offset,
896 auxArray2 = DfluxvectorX1 + phys_offset);
900 DivCFlux_1D(nConvectiveFields, fields, fluxvector[i][0],
908 Vmath::Vadd(nSolutionPts, &outarray[i][0], 1, &DfluxvectorX1[0],
909 1, &outarray[i][0], 1);
912 if (!(
Basis[0]->Collocation()))
914 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
915 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
931 for (i = 0; i < nConvectiveFields; ++i)
938 &fluxvector[i][0][0], 1, &
m_gmat[2][0], 1,
939 &fluxvector[i][1][0], 1, &f_hat[0], 1);
945 &fluxvector[i][0][0], 1, &
m_gmat[3][0], 1,
946 &fluxvector[i][1][0], 1, &g_hat[0], 1);
952 for (n = 0; n < nElements; n++)
954 phys_offset = fields[0]->GetPhys_Offset(n);
955 fields[0]->GetExp(n)->StdPhysDeriv(
956 0, auxArray1 = f_hat + phys_offset,
957 auxArray2 = DfluxvectorX1 + phys_offset);
958 fields[0]->GetExp(n)->StdPhysDeriv(
959 1, auxArray1 = g_hat + phys_offset,
960 auxArray2 = DfluxvectorX2 + phys_offset);
964 Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1, DfluxvectorX2, 1,
968 if (
Basis[0]->GetPointsType() ==
970 Basis[1]->GetPointsType() ==
978 DivCFlux_2D(nConvectiveFields, fields, fluxvector[i][0],
979 fluxvector[i][1], numflux[i], divFC);
983 Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1, outarray[i], 1);
990 if (!(
Basis[0]->Collocation()))
992 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
993 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1001 ASSERTL0(
false,
"3D FRDG case not implemented yet");
1018 [[maybe_unused]]
const int nConvectiveFields,
1025 int nLocalSolutionPts, phys_offset, t_offset;
1028 Basis = fields[0]->GetExp(0)->GetBasis(0);
1030 int nElements = fields[0]->GetExpSize();
1031 int nSolutionPts = fields[0]->GetTotPoints();
1033 vector<bool> leftAdjacentTraces =
1034 std::static_pointer_cast<MultiRegions::DisContField>(fields[0])
1035 ->GetLeftAdjacentTraces();
1046 fields[0]->GetTraceMap()->GetElmtToTrace();
1048 for (n = 0; n < nElements; ++n)
1050 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1051 phys_offset = fields[0]->GetPhys_Offset(n);
1055 Vmath::Vcopy(nLocalSolutionPts, &fluxX1[phys_offset], 1, &tmparrayX1[0],
1058 fields[0]->GetExp(n)->GetTracePhysVals(0, elmtToTrace[n][0], tmparrayX1,
1061 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1062 elmtToTrace[n][0]->GetElmtId());
1064 if (leftAdjacentTraces[2 * n])
1066 JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1070 JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1073 fields[0]->GetExp(n)->GetTracePhysVals(1, elmtToTrace[n][1], tmparrayX1,
1076 t_offset = fields[0]->GetTrace()->GetPhys_Offset(
1077 elmtToTrace[n][1]->GetElmtId());
1079 if (leftAdjacentTraces[2 * n + 1])
1081 JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex[0];
1085 JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex[0];
1089 for (n = 0; n < nElements; ++n)
1091 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1092 phys_offset = fields[0]->GetPhys_Offset(n);
1103 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1104 &divCFlux[phys_offset], 1);
1121 [[maybe_unused]]
const int nConvectiveFields,
1128 int n, e, i, j, cnt;
1130 int nElements = fields[0]->GetExpSize();
1132 int nLocalSolutionPts, nEdgePts;
1133 int trace_offset, phys_offset;
1140 fields[0]->GetTraceMap()->GetElmtToTrace();
1143 for (n = 0; n < nElements; ++n)
1146 phys_offset = fields[0]->GetPhys_Offset(n);
1147 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1149 base = fields[0]->GetExp(n)->GetBase();
1150 nquad0 = base[0]->GetNumPoints();
1151 nquad1 = base[1]->GetNumPoints();
1159 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1162 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1171 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1172 elmtToTrace[n][e]->GetElmtId());
1176 fields[0]->GetExp(n)->GetTraceNormal(e);
1180 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1181 fluxX1 + phys_offset,
1182 auxArray1 = tmparrayX1);
1186 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1187 fluxX2 + phys_offset,
1188 auxArray1 = tmparrayX2);
1197 Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1, &fluxN[0], 1,
1201 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1207 for (i = 0; i < nEdgePts; ++i)
1212 fluxJumps[i] = -fluxJumps[i];
1220 for (i = 0; i < nquad0; ++i)
1223 fluxJumps[i] = -(
m_Q2D_e0[n][i]) * fluxJumps[i];
1225 for (j = 0; j < nquad1; ++j)
1227 cnt = i + j * nquad0;
1228 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
1233 for (i = 0; i < nquad1; ++i)
1236 fluxJumps[i] = (
m_Q2D_e1[n][i]) * fluxJumps[i];
1238 for (j = 0; j < nquad0; ++j)
1240 cnt = (nquad0)*i + j;
1241 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1246 for (i = 0; i < nquad0; ++i)
1249 fluxJumps[i] = (
m_Q2D_e2[n][i]) * fluxJumps[i];
1251 for (j = 0; j < nquad1; ++j)
1253 cnt = j * nquad0 + i;
1254 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1259 for (i = 0; i < nquad1; ++i)
1262 fluxJumps[i] = -(
m_Q2D_e3[n][i]) * fluxJumps[i];
1263 for (j = 0; j < nquad0; ++j)
1265 cnt = j + i * nquad0;
1266 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
1271 ASSERTL0(
false,
"edge value (< 3) is out of range");
1277 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1278 &divCFlux[phys_offset], 1);
1280 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1281 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1283 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1284 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1303 [[maybe_unused]]
const int nConvectiveFields,
1310 int n, e, i, j, cnt;
1312 int nElements = fields[0]->GetExpSize();
1313 int nLocalSolutionPts;
1324 fields[0]->GetTraceMap()->GetElmtToTrace();
1327 for (n = 0; n < nElements; ++n)
1330 phys_offset = fields[0]->GetPhys_Offset(n);
1331 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1333 base = fields[0]->GetExp(n)->GetBase();
1334 nquad0 = base[0]->GetNumPoints();
1335 nquad1 = base[1]->GetNumPoints();
1343 for (e = 0; e < fields[0]->GetExp(n)->GetNtraces(); ++e)
1346 nEdgePts = fields[0]->GetExp(n)->GetTraceNumPoints(e);
1355 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1356 elmtToTrace[n][e]->GetElmtId());
1360 fields[0]->GetExp(n)->GetTraceNormal(e);
1369 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1370 fluxX2 + phys_offset,
1371 auxArray1 = fluxN_D);
1376 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1380 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1384 auxArray2 = fluxN, 1);
1387 auxArray2 = fluxN_D, 1);
1391 for (i = 0; i < nquad0; ++i)
1394 fluxN_R[i] = (
m_Q2D_e0[n][i]) * fluxN[i];
1397 for (i = 0; i < nEdgePts; ++i)
1404 fluxN_R[i] = -fluxN_R[i];
1410 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1415 for (i = 0; i < nquad0; ++i)
1417 for (j = 0; j < nquad1; ++j)
1419 cnt = i + j * nquad0;
1420 divCFluxE0[cnt] = -fluxJumps[i] *
m_dGL_xi2[n][j];
1428 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1429 fluxX1 + phys_offset,
1430 auxArray1 = fluxN_D);
1433 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1437 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1441 auxArray2 = fluxN, 1);
1444 auxArray2 = fluxN_D, 1);
1448 for (i = 0; i < nquad1; ++i)
1451 fluxN_R[i] = (
m_Q2D_e1[n][i]) * fluxN[i];
1454 for (i = 0; i < nEdgePts; ++i)
1461 fluxN_R[i] = -fluxN_R[i];
1467 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1472 for (i = 0; i < nquad1; ++i)
1474 for (j = 0; j < nquad0; ++j)
1476 cnt = (nquad0)*i + j;
1477 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1487 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1488 fluxX2 + phys_offset,
1489 auxArray1 = fluxN_D);
1492 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1496 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1500 auxArray2 = fluxN, 1);
1503 auxArray2 = fluxN_D, 1);
1507 for (i = 0; i < nquad0; ++i)
1510 fluxN_R[i] = (
m_Q2D_e2[n][i]) * fluxN[i];
1513 for (i = 0; i < nEdgePts; ++i)
1520 fluxN_R[i] = -fluxN_R[i];
1527 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1532 for (i = 0; i < nquad0; ++i)
1534 for (j = 0; j < nquad1; ++j)
1536 cnt = j * nquad0 + i;
1537 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1546 fields[0]->GetExp(n)->GetTracePhysVals(e, elmtToTrace[n][e],
1547 fluxX1 + phys_offset,
1548 auxArray1 = fluxN_D);
1552 Vmath::Vcopy(nEdgePts, &numericalFlux[trace_offset], 1,
1556 if (fields[0]->GetExp(n)->GetTraceOrient(e) ==
1560 auxArray2 = fluxN, 1);
1563 auxArray2 = fluxN_D, 1);
1567 for (i = 0; i < nquad1; ++i)
1570 fluxN_R[i] = (
m_Q2D_e3[n][i]) * fluxN[i];
1573 for (i = 0; i < nEdgePts; ++i)
1580 fluxN_R[i] = -fluxN_R[i];
1587 Vmath::Vsub(nEdgePts, &fluxN_R[0], 1, &fluxN_D[0], 1,
1592 for (i = 0; i < nquad1; ++i)
1594 for (j = 0; j < nquad0; ++j)
1596 cnt = j + i * nquad0;
1597 divCFluxE3[cnt] = -fluxJumps[i] *
m_dGL_xi1[n][j];
1602 ASSERTL0(
false,
"edge value (< 3) is out of range");
1608 Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1, &divCFluxE1[0], 1,
1609 &divCFlux[phys_offset], 1);
1611 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1612 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1614 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1615 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1633 [[maybe_unused]]
const int nConvectiveFields,
#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
static AdvectionSharedPtr create(std::string advType)
void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Compute the advection term at each time-step using the Flux Reconstruction approach (FR).
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
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...
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Array< OneD, NekDouble > m_jac
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
AdvectionFR(std::string advType)
AdvectionFR uses the Flux Reconstruction (FR) approach to compute the advection term....
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initiliase AdvectionFR objects and store them before starting the time-stepping.
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
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".
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.
void DivCFlux_3D(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 > &fluxX3, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 3D problems.
static std::string type[]
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
void DivCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 1D problems.
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_gmat
int m_spaceDim
Storage for space dimension. Used for homogeneous extension.
AdvectionFluxVecCB m_fluxVector
Callback function to the flux vector (set when advection is in conservative form).
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
RiemannSolverSharedPtr m_riemann
Riemann solver for DG-type schemes.
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
@ 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 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 Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
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.