41 #include <boost/math/special_functions/gamma.hpp>
51 std::string DiffusionLFR::type[] = {
53 "LFRDG", DiffusionLFR::create),
55 "LFRSD", DiffusionLFR::create),
57 "LFRHU", DiffusionLFR::create),
59 "LFRcmin", DiffusionLFR::create),
61 "LFRcinf", DiffusionLFR::create)};
71 DiffusionLFR::DiffusionLFR(std::string diffType):m_diffType(diffType)
95 int nConvectiveFields = pFields.num_elements();
96 int nDim = pFields[0]->GetCoordim(0);
97 int nSolutionPts = pFields[0]->GetTotPoints();
98 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
125 for (i = 0; i < nConvectiveFields; ++i)
141 for (j = 0; j < nDim; ++j)
181 int nLocalSolutionPts;
182 int nElements = pFields[0]->GetExpSize();
183 int nDimensions = pFields[0]->GetCoordim(0);
184 int nSolutionPts = pFields[0]->GetTotPoints();
185 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
188 for (i = 0; i < nDimensions; ++i)
207 for (n = 0; n < nElements; ++n)
209 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
210 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
211 phys_offset = pFields[0]->GetPhys_Offset(n);
212 jac = pFields[0]->GetExp(n)
215 for (i = 0; i < nLocalSolutionPts; ++i)
217 m_jac[i+phys_offset] = jac[0];
235 for (n = 0; n < nElements; ++n)
237 base = pFields[0]->GetExp(n)->GetBase();
238 nquad0 = base[0]->GetNumPoints();
239 nquad1 = base[1]->GetNumPoints();
247 pFields[0]->GetExp(n)->GetEdgeQFactors(
248 0, auxArray1 = m_Q2D_e0[n]);
249 pFields[0]->GetExp(n)->GetEdgeQFactors(
250 1, auxArray1 = m_Q2D_e1[n]);
251 pFields[0]->GetExp(n)->GetEdgeQFactors(
252 2, auxArray1 = m_Q2D_e2[n]);
253 pFields[0]->GetExp(n)->GetEdgeQFactors(
254 3, auxArray1 = m_Q2D_e3[n]);
256 ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
257 nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
258 phys_offset = pFields[0]->GetPhys_Offset(n);
260 jac = pFields[0]->GetExp(n)
263 gmat = pFields[0]->GetExp(n)
267 if (pFields[0]->GetExp(n)
268 ->as<LocalRegions::Expansion2D>()->GetGeom2D()
269 ->GetMetricInfo()->GetGtype()
272 for (i = 0; i < nLocalSolutionPts; ++i)
274 m_jac[i+phys_offset] = jac[i];
275 m_gmat[0][i+phys_offset] = gmat[0][i];
276 m_gmat[1][i+phys_offset] = gmat[1][i];
277 m_gmat[2][i+phys_offset] = gmat[2][i];
278 m_gmat[3][i+phys_offset] = gmat[3][i];
283 for (i = 0; i < nLocalSolutionPts; ++i)
285 m_jac[i+phys_offset] = jac[0];
286 m_gmat[0][i+phys_offset] = gmat[0][0];
287 m_gmat[1][i+phys_offset] = gmat[1][0];
288 m_gmat[2][i+phys_offset] = gmat[2][0];
289 m_gmat[3][i+phys_offset] = gmat[3][0];
297 ASSERTL0(
false,
"3DFR Metric terms not implemented yet");
302 ASSERTL0(
false,
"Expansion dimension not recognised");
332 int nquad0, nquad1, nquad2;
333 int nmodes0, nmodes1, nmodes2;
336 int nElements = pFields[0]->GetExpSize();
337 int nDim = pFields[0]->GetCoordim(0);
346 for (n = 0; n < nElements; ++n)
348 base = pFields[0]->GetExp(n)->GetBase();
349 nquad0 = base[0]->GetNumPoints();
350 nmodes0 = base[0]->GetNumModes();
354 base[0]->GetZW(z0, w0);
366 int p0 = nmodes0 - 1;
372 NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
374 * boost::math::tgamma(p0 + 1)
375 * boost::math::tgamma(p0 + 1));
384 c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
385 * (ap0 * boost::math::tgamma(p0 + 1))
386 * (ap0 * boost::math::tgamma(p0 + 1)));
390 c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
391 * (ap0 * boost::math::tgamma(p0 + 1))
392 * (ap0 * boost::math::tgamma(p0 + 1)));
396 c0 = -2.0 / ((2.0 * p0 + 1.0)
397 * (ap0 * boost::math::tgamma(p0 + 1))
398 * (ap0 * boost::math::tgamma(p0 + 1)));
402 c0 = 10000000000000000.0;
405 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
406 * (ap0 * boost::math::tgamma(p0 + 1))
407 * (ap0 * boost::math::tgamma(p0 + 1));
409 NekDouble overeta0 = 1.0 / (1.0 + etap0);
420 for(i = 0; i < nquad0; ++i)
426 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
430 for(i = 0; i < nquad0; ++i)
432 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
433 m_dGR_xi1[n][i] += dLpp0[i];
434 m_dGR_xi1[n][i] *= overeta0;
435 m_dGR_xi1[n][i] += dLp0[i];
436 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
453 for (n = 0; n < nElements; ++n)
455 base = pFields[0]->GetExp(n)->GetBase();
456 nquad0 = base[0]->GetNumPoints();
457 nquad1 = base[1]->GetNumPoints();
458 nmodes0 = base[0]->GetNumModes();
459 nmodes1 = base[1]->GetNumModes();
466 base[0]->GetZW(z0, w0);
467 base[1]->GetZW(z1, w1);
484 int p0 = nmodes0 - 1;
485 int p1 = nmodes1 - 1;
492 NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
494 * boost::math::tgamma(p0 + 1)
495 * boost::math::tgamma(p0 + 1));
497 NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
499 * boost::math::tgamma(p1 + 1)
500 * boost::math::tgamma(p1 + 1));
510 c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
511 * (ap0 * boost::math::tgamma(p0 + 1))
512 * (ap0 * boost::math::tgamma(p0 + 1)));
514 c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
515 * (ap1 * boost::math::tgamma(p1 + 1))
516 * (ap1 * boost::math::tgamma(p1 + 1)));
520 c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
521 * (ap0 * boost::math::tgamma(p0 + 1))
522 * (ap0 * boost::math::tgamma(p0 + 1)));
524 c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
525 * (ap1 * boost::math::tgamma(p1 + 1))
526 * (ap1 * boost::math::tgamma(p1 + 1)));
530 c0 = -2.0 / ((2.0 * p0 + 1.0)
531 * (ap0 * boost::math::tgamma(p0 + 1))
532 * (ap0 * boost::math::tgamma(p0 + 1)));
534 c1 = -2.0 / ((2.0 * p1 + 1.0)
535 * (ap1 * boost::math::tgamma(p1 + 1))
536 * (ap1 * boost::math::tgamma(p1 + 1)));
540 c0 = 10000000000000000.0;
541 c1 = 10000000000000000.0;
544 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
545 * (ap0 * boost::math::tgamma(p0 + 1))
546 * (ap0 * boost::math::tgamma(p0 + 1));
548 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
549 * (ap1 * boost::math::tgamma(p1 + 1))
550 * (ap1 * boost::math::tgamma(p1 + 1));
552 NekDouble overeta0 = 1.0 / (1.0 + etap0);
553 NekDouble overeta1 = 1.0 / (1.0 + etap1);
568 for(i = 0; i < nquad0; ++i)
574 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
578 for(i = 0; i < nquad1; ++i)
580 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
581 m_dGL_xi2[n][i] += dLpp1[i];
582 m_dGL_xi2[n][i] *= overeta1;
583 m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
584 m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
588 for(i = 0; i < nquad0; ++i)
590 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
591 m_dGR_xi1[n][i] += dLpp0[i];
592 m_dGR_xi1[n][i] *= overeta0;
593 m_dGR_xi1[n][i] += dLp0[i];
594 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
598 for(i = 0; i < nquad1; ++i)
600 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
601 m_dGR_xi2[n][i] += dLpp1[i];
602 m_dGR_xi2[n][i] *= overeta1;
603 m_dGR_xi2[n][i] += dLp1[i];
604 m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
618 for (n = 0; n < nElements; ++n)
620 base = pFields[0]->GetExp(n)->GetBase();
621 nquad0 = base[0]->GetNumPoints();
622 nquad1 = base[1]->GetNumPoints();
623 nquad2 = base[2]->GetNumPoints();
624 nmodes0 = base[0]->GetNumModes();
625 nmodes1 = base[1]->GetNumModes();
626 nmodes2 = base[2]->GetNumModes();
635 base[0]->GetZW(z0, w0);
636 base[1]->GetZW(z1, w1);
637 base[1]->GetZW(z2, w2);
659 int p0 = nmodes0 - 1;
660 int p1 = nmodes1 - 1;
661 int p2 = nmodes2 - 1;
669 NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
671 * boost::math::tgamma(p0 + 1)
672 * boost::math::tgamma(p0 + 1));
675 NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
677 * boost::math::tgamma(p1 + 1)
678 * boost::math::tgamma(p1 + 1));
681 NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
683 * boost::math::tgamma(p2 + 1)
684 * boost::math::tgamma(p2 + 1));
695 c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
696 * (ap0 * boost::math::tgamma(p0 + 1))
697 * (ap0 * boost::math::tgamma(p0 + 1)));
699 c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
700 * (ap1 * boost::math::tgamma(p1 + 1))
701 * (ap1 * boost::math::tgamma(p1 + 1)));
703 c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
704 * (ap2 * boost::math::tgamma(p2 + 1))
705 * (ap2 * boost::math::tgamma(p2 + 1)));
709 c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
710 * (ap0 * boost::math::tgamma(p0 + 1))
711 * (ap0 * boost::math::tgamma(p0 + 1)));
713 c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
714 * (ap1 * boost::math::tgamma(p1 + 1))
715 * (ap1 * boost::math::tgamma(p1 + 1)));
717 c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
718 * (ap2 * boost::math::tgamma(p2 + 1))
719 * (ap2 * boost::math::tgamma(p2 + 1)));
723 c0 = -2.0 / ((2.0 * p0 + 1.0)
724 * (ap0 * boost::math::tgamma(p0 + 1))
725 * (ap0 * boost::math::tgamma(p0 + 1)));
727 c1 = -2.0 / ((2.0 * p1 + 1.0)
728 * (ap1 * boost::math::tgamma(p1 + 1))
729 * (ap1 * boost::math::tgamma(p1 + 1)));
731 c2 = -2.0 / ((2.0 * p2 + 1.0)
732 * (ap2 * boost::math::tgamma(p2 + 1))
733 * (ap2 * boost::math::tgamma(p2 + 1)));
737 c0 = 10000000000000000.0;
738 c1 = 10000000000000000.0;
739 c2 = 10000000000000000.0;
742 NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
743 * (ap0 * boost::math::tgamma(p0 + 1))
744 * (ap0 * boost::math::tgamma(p0 + 1));
746 NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
747 * (ap1 * boost::math::tgamma(p1 + 1))
748 * (ap1 * boost::math::tgamma(p1 + 1));
750 NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
751 * (ap2 * boost::math::tgamma(p2 + 1))
752 * (ap2 * boost::math::tgamma(p2 + 1));
754 NekDouble overeta0 = 1.0 / (1.0 + etap0);
755 NekDouble overeta1 = 1.0 / (1.0 + etap1);
756 NekDouble overeta2 = 1.0 / (1.0 + etap2);
775 for(i = 0; i < nquad0; ++i)
781 m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
785 for(i = 0; i < nquad1; ++i)
787 m_dGL_xi2[n][i] = etap1 * dLpm1[i];
788 m_dGL_xi2[n][i] += dLpp1[i];
789 m_dGL_xi2[n][i] *= overeta1;
790 m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
791 m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
795 for(i = 0; i < nquad2; ++i)
797 m_dGL_xi3[n][i] = etap2 * dLpm2[i];
798 m_dGL_xi3[n][i] += dLpp2[i];
799 m_dGL_xi3[n][i] *= overeta2;
800 m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
801 m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
805 for(i = 0; i < nquad0; ++i)
807 m_dGR_xi1[n][i] = etap0 * dLpm0[i];
808 m_dGR_xi1[n][i] += dLpp0[i];
809 m_dGR_xi1[n][i] *= overeta0;
810 m_dGR_xi1[n][i] += dLp0[i];
811 m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
815 for(i = 0; i < nquad1; ++i)
817 m_dGR_xi2[n][i] = etap1 * dLpm1[i];
818 m_dGR_xi2[n][i] += dLpp1[i];
819 m_dGR_xi2[n][i] *= overeta1;
820 m_dGR_xi2[n][i] += dLp1[i];
821 m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
825 for(i = 0; i < nquad2; ++i)
827 m_dGR_xi3[n][i] = etap2 * dLpm2[i];
828 m_dGR_xi3[n][i] += dLpp2[i];
829 m_dGR_xi3[n][i] *= overeta2;
830 m_dGR_xi3[n][i] += dLp2[i];
831 m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
838 ASSERTL0(
false,
"Expansion dimension not recognised");
850 const int nConvectiveFields,
864 Basis = fields[0]->GetExp(0)->GetBase();
866 int nElements = fields[0]->GetExpSize();
867 int nDim = fields[0]->GetCoordim(0);
868 int nSolutionPts = fields[0]->GetTotPoints();
869 int nCoeffs = fields[0]->GetNcoeffs();
872 for (i = 0; i < nConvectiveFields; ++i)
885 for (i = 0; i < nConvectiveFields; ++i)
889 for (n = 0; n < nElements; n++)
891 phys_offset = fields[0]->GetPhys_Offset(n);
893 fields[i]->GetExp(n)->PhysDeriv(0,
894 auxArray1 = inarray[i] + phys_offset,
895 auxArray2 =
m_DU1[i][0] + phys_offset);
921 for (i = 0; i < nConvectiveFields; ++i)
925 for (n = 0; n < nElements; n++)
927 phys_offset = fields[0]->GetPhys_Offset(n);
929 fields[i]->GetExp(n)->PhysDeriv(0,
930 auxArray1 =
m_D1[i][0] + phys_offset,
931 auxArray2 =
m_DD1[i][0] + phys_offset);
947 &
m_DD1[i][0][0], 1, &outarray[i][0], 1);
950 if (!(Basis[0]->Collocation()))
952 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
953 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
961 for(i = 0; i < nConvectiveFields; ++i)
963 for (j = 0; j < nDim; ++j)
972 &
m_gmat[0][0], 1, &u1_hat[0], 1);
975 &
m_jac[0], 1, &u1_hat[0], 1);
978 &
m_gmat[1][0], 1, &u2_hat[0], 1);
981 &
m_jac[0], 1, &u2_hat[0], 1);
986 &
m_gmat[2][0], 1, &u1_hat[0], 1);
989 &
m_jac[0], 1, &u1_hat[0], 1);
992 &
m_gmat[3][0], 1, &u2_hat[0], 1);
995 &
m_jac[0], 1, &u2_hat[0], 1);
998 for (n = 0; n < nElements; n++)
1000 phys_offset = fields[0]->GetPhys_Offset(n);
1002 fields[i]->GetExp(n)->StdPhysDeriv(0,
1003 auxArray1 = u1_hat + phys_offset,
1004 auxArray2 =
m_tmp1[i][j] + phys_offset);
1006 fields[i]->GetExp(n)->StdPhysDeriv(1,
1007 auxArray1 = u2_hat + phys_offset,
1008 auxArray2 =
m_tmp2[i][j] + phys_offset);
1013 &
m_DU1[i][j][0], 1);
1022 inarray[i],
m_IF1[i][j],
1028 for (j = 0; j < nSolutionPts; ++j)
1051 for (j = 0; j < nSolutionPts; j++)
1060 m_gmat[0][j] * m_BD1[i][1][j] -
1061 m_gmat[1][j] * m_BD1[i][0][j];
1065 for (j = 0; j < nDim; ++j)
1079 for (i = 0; i < nConvectiveFields; ++i)
1084 for (j = 0; j < nSolutionPts; j++)
1093 for (n = 0; n < nElements; n++)
1095 phys_offset = fields[0]->GetPhys_Offset(n);
1097 fields[0]->GetExp(n)->StdPhysDeriv(0,
1098 auxArray1 = f_hat + phys_offset,
1099 auxArray2 =
m_DD1[i][0] + phys_offset);
1101 fields[0]->GetExp(n)->StdPhysDeriv(1,
1102 auxArray1 = g_hat + phys_offset,
1103 auxArray2 =
m_DD1[i][1] + phys_offset);
1113 if (Basis[0]->GetPointsType() ==
1115 Basis[1]->GetPointsType() ==
1131 &
m_divFC[i][0], 1, &outarray[i][0], 1);
1136 &
m_jac[0], 1, &outarray[i][0], 1);
1139 if (!(Basis[0]->Collocation()))
1141 fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1142 fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1150 ASSERTL0(
false,
"3D FRDG case not implemented yet");
1166 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1167 int nvariables = fields.num_elements();
1168 int nDim = fields[0]->GetCoordim(0);
1176 for (i = 0; i < nDim; ++i)
1185 for (j = 0; j < nDim; ++j)
1187 for (i = 0; i < nvariables ; ++i)
1190 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1200 fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1211 if(fields[0]->GetBndCondExpansions().num_elements())
1243 int nBndEdgePts, nBndEdges;
1245 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1246 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1250 fields[var]->ExtractTracePhys(ufield, uplus);
1253 for (i = 0; i < nBndRegions; ++i)
1256 nBndEdges = fields[var]->
1257 GetBndCondExpansions()[i]->GetExpSize();
1260 for (e = 0; e < nBndEdges ; ++e)
1263 nBndEdgePts = fields[var]->
1264 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1268 GetBndCondExpansions()[i]->GetPhys_Offset(e);
1271 id2 = fields[0]->GetTrace()->
1272 GetPhys_Offset(fields[0]->GetTraceMap()->
1273 GetBndCondTraceToGlobalTraceMap(cnt++));
1276 if (fields[var]->GetBndConditions()[i]->
1281 GetBndCondExpansions()[i]->
1283 &penaltyflux[id2], 1);
1286 else if ((fields[var]->GetBndConditions()[i])->
1291 &penaltyflux[id2], 1);
1308 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1309 int nvariables = fields.num_elements();
1310 int nDim = fields[0]->GetCoordim(0);
1323 for(i = 0; i < nDim; ++i)
1331 for (i = 0; i < nvariables; ++i)
1334 for (j = 0; j < nDim; ++j)
1337 fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1351 fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1376 if (fields[0]->GetBndCondExpansions().num_elements())
1406 int nBndEdges, nBndEdgePts;
1407 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1408 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1414 fields[var]->ExtractTracePhys(qfield, qtemp);
1416 for (i = 0; i < nBndRegions; ++i)
1418 nBndEdges = fields[var]->
1419 GetBndCondExpansions()[i]->GetExpSize();
1422 for (e = 0; e < nBndEdges ; ++e)
1424 nBndEdgePts = fields[var]->
1425 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1428 GetBndCondExpansions()[i]->GetPhys_Offset(e);
1430 id2 = fields[0]->GetTrace()->
1431 GetPhys_Offset(fields[0]->GetTraceMap()->
1432 GetBndCondTraceToGlobalTraceMap(cnt++));
1436 if (fields[var]->GetBndConditions()[i]->
1440 &qtemp[id2], 1, &penaltyflux[id2], 1);
1443 else if ((fields[var]->GetBndConditions()[i])->
1447 &(fields[var]->GetBndCondExpansions()[i]->
1448 GetPhys())[id1], 1, &penaltyflux[id2], 1);
1465 const int nConvectiveFields,
1472 int nLocalSolutionPts, phys_offset, t_offset;
1480 Basis = fields[0]->GetExp(0)->GetBasis(0);
1482 int nElements = fields[0]->GetExpSize();
1483 int nSolutionPts = fields[0]->GetTotPoints();
1496 &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1498 for (n = 0; n < nElements; ++n)
1500 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1501 phys_offset = fields[0]->GetPhys_Offset(n);
1506 &flux[phys_offset], 1,
1509 fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1512 t_offset = fields[0]->GetTrace()
1513 ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1515 if(negatedFluxNormal[2*n])
1517 JumpL[n] = iFlux[t_offset] - tmpFluxVertex;
1521 JumpL[n] = -iFlux[t_offset] - tmpFluxVertex;
1524 t_offset = fields[0]->GetTrace()
1525 ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1527 fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1529 if(negatedFluxNormal[2*n+1])
1531 JumpR[n] = -iFlux[t_offset] - tmpFluxVertex;
1535 JumpR[n] = iFlux[t_offset] - tmpFluxVertex;
1539 for (n = 0; n < nElements; ++n)
1541 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1542 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1543 phys_offset = fields[0]->GetPhys_Offset(n);
1544 jac = fields[0]->GetExp(n)
1548 JumpL[n] = JumpL[n] * jac[0];
1549 JumpR[n] = JumpR[n] * jac[0];
1560 Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1561 &derCFlux[phys_offset], 1);
1581 const int nConvectiveFields,
1582 const int direction,
1588 int n, e, i, j, cnt;
1593 int nElements = fields[0]->GetExpSize();
1595 int trace_offset, phys_offset;
1596 int nLocalSolutionPts;
1604 &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1607 for (n = 0; n < nElements; ++n)
1610 phys_offset = fields[0]->GetPhys_Offset(n);
1611 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1612 ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1617 base = fields[0]->GetExp(n)->GetBase();
1618 nquad0 = base[0]->GetNumPoints();
1619 nquad1 = base[1]->GetNumPoints();
1627 for (e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1630 nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1636 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1637 elmtToTrace[n][e]->GetElmtId());
1646 fields[0]->GetExp(n)->GetEdgePhysVals(e, elmtToTrace[n][e],
1648 auxArray1 = tmparray);
1654 &tmparray[0], 1, &fluxJumps[0], 1);
1658 if (fields[0]->GetExp(n)->GetEorient(e) ==
1667 if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1668 ->GetGeom2D()->GetMetricInfo()->GetGtype()
1674 fields[0]->GetExp(n)->GetEdgePhysVals(
1675 e, elmtToTrace[n][e],
1676 jac, auxArray1 = jacEdge);
1680 if (fields[0]->GetExp(n)->GetEorient(e) ==
1690 for (j = 0; j < nEdgePts; j++)
1692 fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1710 for (i = 0; i < nquad0; ++i)
1715 for (j = 0; j < nquad1; ++j)
1718 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
1723 for (i = 0; i < nquad1; ++i)
1728 for (j = 0; j < nquad0; ++j)
1730 cnt = (nquad0)*i + j;
1731 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1736 for (i = 0; i < nquad0; ++i)
1741 for (j = 0; j < nquad1; ++j)
1743 cnt = (nquad0 - 1) + j*nquad0 - i;
1744 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1749 for (i = 0; i < nquad1; ++i)
1753 for (j = 0; j < nquad0; ++j)
1755 cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1756 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
1762 ASSERTL0(
false,
"edge value (< 3) is out of range");
1772 &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1774 else if (direction == 1)
1777 &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1796 const int nConvectiveFields,
1803 int n, e, i, j, cnt;
1805 int nElements = fields[0]->GetExpSize();
1807 int nLocalSolutionPts;
1818 &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1821 for(n = 0; n < nElements; ++n)
1824 phys_offset = fields[0]->GetPhys_Offset(n);
1825 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1827 base = fields[0]->GetExp(n)->GetBase();
1828 nquad0 = base[0]->GetNumPoints();
1829 nquad1 = base[1]->GetNumPoints();
1837 for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1840 nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1849 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1850 elmtToTrace[n][e]->GetElmtId());
1854 fields[0]->GetExp(n)->GetEdgeNormal(e);
1858 fields[0]->GetExp(n)->GetEdgePhysVals(
1859 e, elmtToTrace[n][e],
1860 fluxX1 + phys_offset,
1861 auxArray1 = tmparrayX1);
1865 fields[0]->GetExp(n)->GetEdgePhysVals(
1866 e, elmtToTrace[n][e],
1867 fluxX2 + phys_offset,
1868 auxArray1 = tmparrayX2);
1871 for (i = 0; i < nEdgePts; ++i)
1880 &numericalFlux[trace_offset], 1,
1881 &fluxN[0], 1, &fluxJumps[0], 1);
1884 if (fields[0]->GetExp(n)->GetEorient(e) ==
1888 auxArray1 = fluxJumps, 1,
1889 auxArray2 = fluxJumps, 1);
1892 NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1895 for (i = 0; i < nEdgePts; ++i)
1900 fluxJumps[i] = -fluxJumps[i];
1908 for (i = 0; i < nquad0; ++i)
1911 fluxJumps[i] = -(
m_Q2D_e0[n][i]) * fluxJumps[i];
1913 for (j = 0; j < nquad1; ++j)
1916 divCFluxE0[cnt] = fluxJumps[i] *
m_dGL_xi2[n][j];
1921 for (i = 0; i < nquad1; ++i)
1924 fluxJumps[i] = (
m_Q2D_e1[n][i]) * fluxJumps[i];
1926 for (j = 0; j < nquad0; ++j)
1928 cnt = (nquad0)*i + j;
1929 divCFluxE1[cnt] = fluxJumps[i] *
m_dGR_xi1[n][j];
1934 for (i = 0; i < nquad0; ++i)
1937 fluxJumps[i] = (
m_Q2D_e2[n][i]) * fluxJumps[i];
1939 for (j = 0; j < nquad1; ++j)
1941 cnt = (nquad0 - 1) + j*nquad0 - i;
1942 divCFluxE2[cnt] = fluxJumps[i] *
m_dGR_xi2[n][j];
1947 for (i = 0; i < nquad1; ++i)
1950 fluxJumps[i] = -(
m_Q2D_e3[n][i]) * fluxJumps[i];
1951 for (j = 0; j < nquad0; ++j)
1953 cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1954 divCFluxE3[cnt] = fluxJumps[i] *
m_dGL_xi1[n][j];
1960 ASSERTL0(
false,
"edge value (< 3) is out of range");
1967 &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1969 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1970 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1972 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1973 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1992 const int nConvectiveFields,
1999 int n, e, i, j, cnt;
2001 int nElements = fields[0]->GetExpSize();
2002 int nLocalSolutionPts;
2013 &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2016 for(n = 0; n < nElements; ++n)
2019 phys_offset = fields[0]->GetPhys_Offset(n);
2020 nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2022 base = fields[0]->GetExp(n)->GetBase();
2023 nquad0 = base[0]->GetNumPoints();
2024 nquad1 = base[1]->GetNumPoints();
2032 for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
2035 nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
2044 trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2045 elmtToTrace[n][e]->GetElmtId());
2049 fields[0]->GetExp(n)->GetEdgeNormal(e);
2058 fields[0]->GetExp(n)->GetEdgePhysVals(
2059 e, elmtToTrace[n][e],
2060 fluxX2 + phys_offset,
2061 auxArray1 = fluxN_D);
2067 &numericalFlux[trace_offset], 1,
2071 if (fields[0]->GetExp(n)->GetEorient(e) ==
2075 auxArray1 = fluxN, 1,
2076 auxArray2 = fluxN, 1);
2079 auxArray1 = fluxN_D, 1,
2080 auxArray2 = fluxN_D, 1);
2084 for (i = 0; i < nquad0; ++i)
2087 fluxN_R[i] = (
m_Q2D_e0[n][i]) * fluxN[i];
2090 for (i = 0; i < nEdgePts; ++i)
2097 fluxN_R[i] = -fluxN_R[i];
2105 &fluxN_D[0], 1, &fluxJumps[0], 1);
2109 for (i = 0; i < nquad0; ++i)
2111 for (j = 0; j < nquad1; ++j)
2123 fields[0]->GetExp(n)->GetEdgePhysVals(
2124 e, elmtToTrace[n][e],
2125 fluxX1 + phys_offset,
2126 auxArray1 = fluxN_D);
2130 &numericalFlux[trace_offset], 1,
2134 if (fields[0]->GetExp(n)->GetEorient(e) ==
2138 auxArray1 = fluxN, 1,
2139 auxArray2 = fluxN, 1);
2142 auxArray1 = fluxN_D, 1,
2143 auxArray2 = fluxN_D, 1);
2147 for (i = 0; i < nquad1; ++i)
2150 fluxN_R[i] = (
m_Q2D_e1[n][i]) * fluxN[i];
2153 for (i = 0; i < nEdgePts; ++i)
2160 fluxN_R[i] = -fluxN_R[i];
2168 &fluxN_D[0], 1, &fluxJumps[0], 1);
2172 for (i = 0; i < nquad1; ++i)
2174 for (j = 0; j < nquad0; ++j)
2176 cnt = (nquad0)*i + j;
2188 fields[0]->GetExp(n)->GetEdgePhysVals(
2189 e, elmtToTrace[n][e],
2190 fluxX2 + phys_offset,
2191 auxArray1 = fluxN_D);
2195 &numericalFlux[trace_offset], 1,
2199 if (fields[0]->GetExp(n)->GetEorient(e) ==
2203 auxArray1 = fluxN, 1,
2204 auxArray2 = fluxN, 1);
2207 auxArray1 = fluxN_D, 1,
2208 auxArray2 = fluxN_D, 1);
2212 for (i = 0; i < nquad0; ++i)
2215 fluxN_R[i] = (
m_Q2D_e2[n][i]) * fluxN[i];
2218 for (i = 0; i < nEdgePts; ++i)
2225 fluxN_R[i] = -fluxN_R[i];
2234 &fluxN_D[0], 1, &fluxJumps[0], 1);
2238 for (i = 0; i < nquad0; ++i)
2240 for (j = 0; j < nquad1; ++j)
2242 cnt = (nquad0 - 1) + j*nquad0 - i;
2253 fields[0]->GetExp(n)->GetEdgePhysVals(
2254 e, elmtToTrace[n][e],
2255 fluxX1 + phys_offset,
2256 auxArray1 = fluxN_D);
2261 &numericalFlux[trace_offset], 1,
2265 if (fields[0]->GetExp(n)->GetEorient(e) ==
2269 auxArray1 = fluxN, 1,
2270 auxArray2 = fluxN, 1);
2273 auxArray1 = fluxN_D, 1,
2274 auxArray2 = fluxN_D, 1);
2278 for (i = 0; i < nquad1; ++i)
2281 fluxN_R[i] = (
m_Q2D_e3[n][i]) * fluxN[i];
2284 for (i = 0; i < nEdgePts; ++i)
2291 fluxN_R[i] = -fluxN_R[i];
2300 &fluxN_D[0], 1, &fluxJumps[0], 1);
2304 for (i = 0; i < nquad1; ++i)
2306 for (j = 0; j < nquad0; ++j)
2308 cnt = (nquad0*nquad1 - nquad0) + j
2316 ASSERTL0(
false,
"edge value (< 3) is out of range");
2324 &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2326 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2327 &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2329 Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2330 &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
virtual void v_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_dGL_xi3
#define ASSERTL0(condition, msg)
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
std::vector< PointsKey > PointsKeyVector
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
virtual void v_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.
DiffusionFactory & GetDiffusionFactory()
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
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
1D Gauss-Gauss-Legendre quadrature points
Array< OneD, NekDouble > m_jac
Array< OneD, Array< OneD, NekDouble > > m_divFC
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initiliase DiffusionLFR objects and store them before starting the time-stepping. ...
virtual void v_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–7...
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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.
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
virtual void v_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.
virtual void v_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.
This class is the abstraction of a global discontinuous two- dimensional spectral/hp element expansio...
void Neg(int n, T *x, const int incx)
Negate x = -x.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
virtual void v_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...
virtual void v_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.
virtual void v_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.
Array< OneD, Array< OneD, NekDouble > > m_gmat
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
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.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
virtual void v_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.
const boost::shared_ptr< SpatialDomains::GeomFactors > & GetMetricInfo(void) const
virtual void v_Diffuse(const int 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=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
Calculate FR Diffusion for the linear problems using an LDG interface flux.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
boost::shared_ptr< Basis > BasisSharedPtr
Array< OneD, Array< OneD, NekDouble > > m_divFD
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Geometry is curved or has non-constant factors.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
LibUtilities::SessionReaderSharedPtr m_session
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
virtual void v_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".
Array< OneD, Array< OneD, NekDouble > > m_IF2
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.