43 NekDouble Extrapolate::StifflyStable_Betaq_Coeffs[3][3] = {
44 { 1.0, 0.0, 0.0},{ 2.0, -1.0, 0.0},{ 3.0, -3.0, 1.0}};
45 NekDouble Extrapolate::StifflyStable_Alpha_Coeffs[3][3] = {
46 { 1.0, 0.0, 0.0},{ 2.0, -0.5, 0.0},{ 3.0, -1.5, 1.0/3.0}};
47 NekDouble Extrapolate::StifflyStable_Gamma0_Coeffs[3] = {
55 Loki::SingleThreaded > Type;
56 return Type::Instance();
59 Extrapolate::Extrapolate(
65 : m_session(pSession),
67 m_pressure(pPressure),
69 m_advObject(advObject)
81 "StandardExtrapolate",
"StandardExtrapolate");
146 for(n = cnt = 0; n <
m_PBndConds.num_elements(); ++n)
151 m_fields[0]->GetBndElmtExpansion(n, BndElmtExp,
false);
153 int nq = BndElmtExp->GetTotPoints();
168 m_fields[0]->ExtractPhysToBndElmt(n, fields[i],Velocity[i]);
175 m_fields[0]->ExtractPhysToBndElmt(n, N[i], Advection[i]);
180 BndElmtExp->CurlCurl(Velocity, Q);
194 m_fields[0]->ExtractElmtToBndPhys(n, Q[i],BndValues[i]);
197 m_PBndExp[n]->NormVectorIProductWRTBase(BndValues, Pvals);
235 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
240 m_fields[0]->GetBndElmtExpansion(n, BndElmtExp,
false);
242 int nq = BndElmtExp->GetTotPoints();
247 m_fields[0]->ExtractPhysToBndElmt(n, fields[i],
259 BndElmtExp->HomogeneousBwdTrans(Velocity[i],
262 BndElmtExp->SetWaveSpace(
false);
267 m_fields[0]->GetBoundaryNormals(n, normals);
282 BndElmtExp->PhysDeriv(Velocity[i], grad[0], grad[1]);
286 BndElmtExp->PhysDeriv(Velocity[i], grad[0], grad[1],
292 m_fields[0]->ExtractElmtToBndPhys(n, grad[j],bndVal);
302 nGradUn, 1, nGradUn, 1);
311 m_fields[0]->ExtractElmtToBndPhys(n, Velocity[i],u[i]);
327 for(
int i = 0; i < nqb; i++)
348 Vmath::Vvtvp(nqb, u[i], 1, bndVal, 1, E[i], 1, E[i], 1);
364 m_UBndExp[i][n]->GetPhys(),
381 m_PBndExp[n]->HomogeneousFwdTrans(pbc, bndVal);
393 int nbcoeffs =
m_PBndExp[n]->GetNcoeffs();
397 m_PBndExp[n]->HomogeneousFwdTrans(pbc, bndVal);
398 m_PBndExp[n]->IProductWRTBase(bndVal,bndCoeffs);
402 m_PBndExp[n]->IProductWRTBase(pbc,bndCoeffs);
406 bndCoeffs, 1,
m_PBndExp[n]->UpdateCoeffs(),1,
412 m_fields[0]->ExtractElmtToBndPhys(n,
414 m_outflowVel[cnt][i][0],
463 u[i], 1,divU,1,divU,1);
466 if (
m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
468 m_houtflow->m_UBndExp[i][n]->HomogeneousFwdTrans(divU,
472 m_houtflow->m_UBndExp[i][n]->IProductWRTBase(divU,
492 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
506 m_PBndExp[n]->HomogeneousBwdTrans(pbc, pbc);
514 m_fields[0]->GetBoundaryNormals(n, normals);
523 if (
m_houtflow->m_UBndExp[i][n]->GetWaveSpace())
526 HomogeneousFwdTrans(wk, wk);
528 m_houtflow->m_UBndExp[i][n]->IProductWRTBase(wk,wk1);
532 m_houtflow->m_UBndExp[i][n]->UpdateCoeffs(),1);
548 for(n = cnt = 0; n <
m_PBndConds.num_elements(); ++n)
555 m_fields[0]->ExtractPhysToBnd(n, Vel[i], velbc[i]);
557 IProdVnTmp = IProdVn + cnt;
558 m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
584 for(n = cnt = 0; n <
m_PBndConds.num_elements(); ++n)
592 (VelBndExp[i][n]->GetTotPoints(), 0.0);
593 VelBndExp[i][n]->SetWaveSpace(
595 VelBndExp[i][n]->BwdTrans(VelBndExp[i][n]->GetCoeffs(),
598 IProdVnTmp = IProdVn + cnt;
599 m_PBndExp[n]->NormVectorIProductWRTBase(velbc, IProdVnTmp);
618 int nlevels = input.num_elements();
622 tmp = input[nlevels-1];
624 for(
int n = nlevels-1; n > 0; --n)
626 input[n] = input[n-1];
651 int outHBCnumber = 0;
652 int numOutHBCPts = 0;
658 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
662 m_HBCnumber +=
m_PBndExp[n]->GetExpSize();
673 m_HBCnumber +=
m_PBndExp[n]->GetExpSize();
674 numOutHBCPts +=
m_PBndExp[n]->GetTotPoints();
678 else if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
682 numOutHBCPts +=
m_PBndExp[n]->GetTotPoints();
723 ASSERTL0(0,
"Dimension not supported");
741 for(n = 0, cnt = 0; n <
m_PBndConds.num_elements(); ++n)
743 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
753 m_fields[0]->GetBndElmtExpansion(n, BndElmtExp,
false);
755 int nq = BndElmtExp->GetTotPoints();
781 if(
m_houtflow->m_pressurePrimCoeff.num_elements() == 0)
799 boost::static_pointer_cast<
812 boost::static_pointer_cast<
814 >(UBndConds[n])->m_robinPrimitiveCoeff;
819 ASSERTL1(UBndConds[n]->GetBoundaryConditionType()
821 "conditions to be of Robin type when pressure"
822 "outflow is specticied as Robin Boundary type");
842 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
852 std::string primcoeff =
m_houtflow->m_defVelPrimCoeff[i] +
"*" +
857 boost::dynamic_pointer_cast<
862 m_session,rcond->m_robinFunction.GetExpression(),
864 rcond->GetUserDefined(),
867 UBndConds[n] = bcond;
883 int n_points_0 =
m_fields[0]->GetExp(0)->GetTotPoints();
884 int n_element =
m_fields[0]->GetExpSize();
885 int nvel = inarray.num_elements();
896 for (
int i = 0; i < nvel; ++i)
902 for (
int el = 0; el < n_element; ++el)
904 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
905 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
908 if(n_points != n_points_0)
910 for (
int j = 0; j < nvel; ++j)
914 n_points_0 = n_points;
918 for (
int j = 0; j < nvel; ++j)
925 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
927 if (
m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
930 for(
int j = 0; j < nvel; ++j)
932 for(
int k = 0; k < nvel; ++k)
935 tmp = inarray[k] + cnt, 1,
943 for(
int j = 0; j < nvel; ++j)
945 for(
int k = 0; k < nvel; ++k)
948 tmp = inarray[k] + cnt, 1,
957 Vmath::Vmul( n_points, stdVelocity[0], 1, stdVelocity[0], 1,
959 for(
int k = 1; k < nvel; ++k)
966 pntVelocity =
Vmath::Vmax( n_points, stdVelocity[0], 1);
967 maxV[el] = sqrt(pntVelocity);
989 int nlevels = array.num_elements();
990 int nPts = array[0].num_elements();
998 array[nlevels-1], 1);
1000 for(
int n = 0; n < nint-1; ++n)
1003 array[n],1, array[nlevels-1],1,
1004 array[nlevels-1],1);
1019 int nlevels = array.num_elements();
1020 int nPts = array[0].num_elements();
1028 array[nlevels-1], 1);
1030 for(
int n = 0; n < nint-1; ++n)
1033 array[n],1, array[nlevels-1],1,
1034 array[nlevels-1],1);
1047 int nlevels = array.num_elements();
1048 int nPts = array[0].num_elements();
1064 accelerationTerm, 1);
1066 for(
int i = 0; i < acc_order; i++)
1071 accelerationTerm, 1,
1072 accelerationTerm, 1);
1075 array[nlevels-1] = accelerationTerm;
1082 for(cnt = n = 0; n <
m_PBndConds.num_elements(); ++n)
#define ASSERTL0(condition, msg)
std::vector< PointsKey > PointsKeyVector
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
ExtrapolateFactory & GetExtrapolateFactory()
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
std::string GetExpression(void) const
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
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
NekDouble Evaluate() const
LibUtilities::NekFactory< std::string, Extrapolate, const LibUtilities::SessionReaderSharedPtr &, Array< OneD, MultiRegions::ExpListSharedPtr > &, MultiRegions::ExpListSharedPtr &, const Array< OneD, int > &, const SolverUtils::AdvectionSharedPtr & > ExtrapolateFactory
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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 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.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
void Zero(int n, T *x, const int incx)
Zero vector.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Geometry is curved or has non-constant factors.
boost::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
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.
Defines a callback function which evaluates the flux vector.
Provides a generic Factory class.