42 NekDouble Extrapolate::StifflyStable_Betaq_Coeffs[3][3] = {
43 {1.0, 0.0, 0.0}, {2.0, -1.0, 0.0}, {3.0, -3.0, 1.0}};
44 NekDouble Extrapolate::StifflyStable_Alpha_Coeffs[3][3] = {
45 {1.0, 0.0, 0.0}, {2.0, -0.5, 0.0}, {3.0, -1.5, 1.0 / 3.0}};
46 NekDouble Extrapolate::StifflyStable_Gamma0_Coeffs[3] = {1.0, 1.5, 11.0 / 6.0};
59 : m_session(pSession), m_fields(pFields), m_pressure(pPressure),
60 m_velocity(pVel), m_advObject(advObject)
72 "StandardExtrapolate",
"StandardExtrapolate");
139 m_fields[0]->GetBndElmtExpansion(n, BndElmtExp,
false);
141 int nq = BndElmtExp->GetTotPoints();
156 m_fields[0]->ExtractPhysToBndElmt(n, fields[i], Velocity[i]);
168 BndElmtExp->CurlCurl(Velocity, Q);
182 m_fields[0]->ExtractElmtToBndPhys(n, Q[i], BndValues[i]);
185 m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Pvals);
227 m_fields[0]->GetBndElmtExpansion(n, BndElmtExp,
false);
229 int nq = BndElmtExp->GetTotPoints();
246 BndElmtExp->HomogeneousBwdTrans(Velocity[i], Velocity[i]);
248 BndElmtExp->SetWaveSpace(
false);
253 m_fields[0]->GetBoundaryNormals(n, normals);
268 BndElmtExp->PhysDeriv(Velocity[i], grad[0], grad[1]);
272 BndElmtExp->PhysDeriv(Velocity[i], grad[0], grad[1],
278 m_fields[0]->ExtractElmtToBndPhys(n, grad[j], bndVal);
285 Vmath::Vmul(nqb, normals[i], 1, bndVal, 1, bndVal, 1);
296 m_fields[0]->ExtractElmtToBndPhys(n, Velocity[i], u[i]);
304 Vmath::Vvtvp(nqb, normals[i], 1, u[i], 1, un, 1, un, 1);
310 for (
int i = 0; i < nqb; i++)
312 S0[i] = 0.5 * (1.0 - tanh(un[i] / (
m_houtflow->m_U0 *
332 Vmath::Vvtvp(nqb, u[i], 1, bndVal, 1, E[i], 1, E[i], 1);
348 m_houtflow->m_UBndExp[i][n]->GetPhys(), 1, bndVal,
351 Vmath::Vvtvp(nqb, normals[i], 1, bndVal, 1, En, 1, En, 1);
363 m_PBndExp[n]->HomogeneousFwdTrans(pbc, bndVal);
375 int nbcoeffs =
m_PBndExp[n]->GetNcoeffs();
379 m_PBndExp[n]->HomogeneousFwdTrans(pbc, bndVal);
380 m_PBndExp[n]->IProductWRTBase(bndVal, bndCoeffs);
384 m_PBndExp[n]->IProductWRTBase(pbc, bndCoeffs);
389 bndCoeffs, 1,
m_PBndExp[n]->UpdateCoeffs(), 1,
423 Vmath::Smul(nqb, -1.0 * kinvis, divU, 1, bndVal, 1);
431 Vmath::Vvtvp(nqb, normals[i], 1, bndVal, 1, E[i], 1, divU, 1);
434 m_houtflow->m_UBndExp[i][n]->GetPhys(), 1, divU, 1);
441 u[i], 1, divU, 1, divU, 1);
444 if (
m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
446 m_houtflow->m_UBndExp[i][n]->HomogeneousFwdTrans(divU,
451 divU,
m_houtflow->m_UBndExp[i][n]->UpdateCoeffs());
481 m_PBndExp[n]->HomogeneousBwdTrans(pbc, pbc);
489 m_fields[0]->GetBoundaryNormals(n, normals);
498 if (
m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
500 m_houtflow->m_UBndExp[i][n]->HomogeneousFwdTrans(wk, wk);
502 m_houtflow->m_UBndExp[i][n]->IProductWRTBase(wk, wk1);
506 m_houtflow->m_UBndExp[i][n]->UpdateCoeffs(), 1);
527 m_fields[0]->ExtractPhysToBnd(n, Vel[i], velbc[i]);
529 IProdVnTmp = IProdVn + cnt;
530 m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
565 VelBndExp[i][n]->GetTotPoints(), 0.0);
566 VelBndExp[i][n]->SetWaveSpace(
568 VelBndExp[i][n]->BwdTrans(VelBndExp[i][n]->GetCoeffs(),
571 IProdVnTmp = IProdVn + cnt;
572 m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
591 int nlevels = input.size();
595 tmp = input[nlevels - 1];
597 for (
int n = nlevels - 1; n > 0; --n)
599 input[n] = input[n - 1];
624 int outHBCnumber = 0;
625 int numOutHBCPts = 0;
631 if (boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
641 boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
646 numOutHBCPts +=
m_PBndExp[n]->GetTotPoints();
650 else if (boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
653 numOutHBCPts +=
m_PBndExp[n]->GetTotPoints();
694 ASSERTL0(0,
"Dimension not supported");
699 if (numOutHBCPts > 0)
702 numOutHBCPts, outHBCnumber,
m_curl_dim, pSession);
715 if (boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
725 m_fields[0]->GetBndElmtExpansion(n, BndElmtExp,
false);
727 int nq = BndElmtExp->GetTotPoints();
753 if (
m_houtflow->m_pressurePrimCoeff.size() == 0)
768 std::static_pointer_cast<
770 ->m_robinPrimitiveCoeff;
781 std::static_pointer_cast<
784 ->m_robinPrimitiveCoeff;
788 ASSERTL1(UBndConds[n]->GetBoundaryConditionType() ==
791 "conditions to be of Robin type when pressure"
792 "outflow is specticied as Robin Boundary type");
820 std::string primcoeff =
822 boost::lexical_cast<std::string>(
831 m_session, rcond->m_robinFunction.GetExpression(),
832 primcoeff, rcond->GetUserDefined(),
835 UBndConds[n] = bcond;
850 int n_points_0 =
m_fields[0]->GetExp(0)->GetTotPoints();
851 int n_element =
m_fields[0]->GetExpSize();
852 int nvel = inarray.size();
863 for (
int i = 0; i < nvel; ++i)
869 for (
int el = 0; el < n_element; ++el)
871 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
872 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
875 if (n_points != n_points_0)
877 for (
int j = 0; j < nvel; ++j)
881 n_points_0 = n_points;
885 for (
int j = 0; j < nvel; ++j)
895 ->GetDerivFactors(ptsKeys);
897 if (
m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype() ==
900 for (
int j = 0; j < nvel; ++j)
902 for (
int k = 0; k < nvel; ++k)
905 tmp = inarray[k] + cnt, 1, stdVelocity[j], 1,
912 for (
int j = 0; j < nvel; ++j)
914 for (
int k = 0; k < nvel; ++k)
917 tmp = inarray[k] + cnt, 1, stdVelocity[j], 1,
925 Vmath::Vmul(n_points, stdVelocity[0], 1, stdVelocity[0], 1,
927 for (
int k = 1; k < nvel; ++k)
929 Vmath::Vvtvp(n_points, stdVelocity[k], 1, stdVelocity[k], 1,
930 stdVelocity[0], 1, stdVelocity[0], 1);
932 pntVelocity =
Vmath::Vmax(n_points, stdVelocity[0], 1);
933 maxV[el] =
sqrt(pntVelocity);
953 int nlevels = array.size();
954 int nPts = array[0].size();
961 array[nint - 1], 1, array[nlevels - 1], 1);
963 for (
int n = 0; n < nint - 1; ++n)
966 array[nlevels - 1], 1, array[nlevels - 1], 1);
979 int nlevels = array.size();
980 int nPts = array[0].size();
987 array[nint - 1], 1, array[nlevels - 1], 1);
989 for (
int n = 0; n < nint - 1; ++n)
992 array[nlevels - 1], 1, array[nlevels - 1], 1);
1004 int nlevels = array.size();
1005 int nPts = array[0].size();
1018 array[0], 1, accelerationTerm, 1);
1020 for (
int i = 0; i < acc_order; i++)
1024 array[i + 1], 1, accelerationTerm, 1, accelerationTerm, 1);
1027 array[nlevels - 1] = accelerationTerm;
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
std::string GetExpression(void) const
NekDouble Evaluate() const
Provides a generic Factory class.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
An abstract base class encapsulating the concept of advection of a vector field.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::vector< PointsKey > PointsKeyVector
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
The above copyright notice and this permission notice shall be included.
ExtrapolateFactory & GetExtrapolateFactory()
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 Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Svtvm(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 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 Zero(int n, T *x, const int incx)
Zero vector.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
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.
scalarT< T > sqrt(scalarT< T > in)