57 m_session->LoadSolverInfo(
"ShockCaptureType",
62 int nDim = pFields[0]->GetCoordim(0);
63 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
66 for(i = 0; i < nDim; ++i)
74 const int nConvectiveFields,
81 int nBndEdgePts, i, j, k, e;
82 int nDim = fields[0]->GetCoordim(0);
83 int nPts = fields[0]->GetTotPoints();
84 int nCoeffs = fields[0]->GetNcoeffs();
85 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
95 for (j = 0; j < nDim; ++j)
104 for (i = 0; i < nConvectiveFields; ++i)
112 for (k = 0; k < nDim; ++k)
122 for (j = 0; j < nDim; ++j)
124 for (i = 0; i < nConvectiveFields; ++i)
126 fields[i]->IProductWRTDerivBase(j, inarray[i], qcoeffs);
128 fields[i]->AddTraceIntegral (flux[j][i], qcoeffs);
129 fields[i]->SetPhysState (
false);
130 fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);
131 fields[i]->BwdTrans (qcoeffs, qfield[j][i]);
143 int numConvFields = nConvectiveFields;
147 numConvFields = nConvectiveFields - 1;
150 for (j = 0; j < nDim; ++j)
152 for (i = 0; i < numConvFields; ++i)
154 Vmath::Vmul(nPts,qfield[j][i],1,muvar,1,qfield[j][i],1);
161 fields[0]->GetFwdBwdTracePhys(muvar,FwdMuVar,BwdMuVar);
163 int nBndRegions = fields[0]->GetBndCondExpansions().
167 for (i = 0; i < nBndRegions; ++i)
170 int nBndEdges = fields[0]->
171 GetBndCondExpansions()[i]->GetExpSize();
175 for (e = 0; e < nBndEdges ; ++e)
177 nBndEdgePts = fields[0]->GetBndCondExpansions()[i]
178 ->GetExp(e)->GetTotPoints();
180 int id2 = fields[0]->GetTrace()->GetPhys_Offset(
181 fields[0]->GetTraceMap()
182 ->GetBndCondTraceToGlobalTraceMap(cnt++));
184 for (k = 0; k < nBndEdgePts; ++k)
186 BwdMuVar[id2+k] = 0.0;
191 for(i = 0; i < numConvFields; ++i)
193 for(k = 0; k < nTracePts; ++k)
196 0.5 * (FwdMuVar[k] + BwdMuVar[k]) * flux[0][i][k];
204 for (i = 0; i < nConvectiveFields; ++i)
206 for (j = 0; j < nDim; ++j)
208 qdbase[j] = qfield[j][i];
210 fields[i]->IProductWRTDerivBase(qdbase,tmp);
214 fields[i]->AddTraceIntegral (flux[0][i], tmp);
215 fields[i]->SetPhysState (
false);
216 fields[i]->MultiplyByElmtInvMass(tmp, tmp);
217 fields[i]->BwdTrans (tmp, outarray[i]);
229 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
230 int nvariables = fields.num_elements();
231 int nDim = fields[0]->GetCoordim(0);;
239 for(i = 0; i < nDim; ++i)
248 for (i = 0; i < nvariables ; ++i)
254 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
269 fields[i]->GetTrace()->Upwind(Vn,
280 if(fields[0]->GetBndCondExpansions().num_elements())
285 for (j = 0; j < nDim; ++j)
316 int nBndEdgePts, nBndEdges;
318 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
319 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
322 fields[var]->ExtractTracePhys(ufield, uplus);
323 for (i = 0; i < nBndRegions; ++i)
326 nBndEdges = fields[var]->
327 GetBndCondExpansions()[i]->GetExpSize();
330 for (e = 0; e < nBndEdges ; ++e)
332 nBndEdgePts = fields[var]->
333 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
336 GetBndCondExpansions()[i]->GetPhys_Offset(e);
338 id2 = fields[0]->GetTrace()->
339 GetPhys_Offset(fields[0]->GetTraceMap()->
340 GetBndCondTraceToGlobalTraceMap(cnt++));
343 if (fields[var]->GetBndConditions()[i]->
348 GetBndCondExpansions()[i]->
350 &penaltyflux[id2], 1);
353 else if ((fields[var]->GetBndConditions()[i])->
358 &penaltyflux[id2], 1);
373 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
374 int nvariables = fields.num_elements();
375 int nDim = qfield.num_elements();
396 for(i = 0; i < nDim; ++i)
403 for (i = 0; i < nvariables; ++i)
406 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
413 -1.0 * C11, uterm, 1,
417 for (j = 0; j < nDim; ++j)
420 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
434 fields[i]->GetTrace()->Upwind(Vn,
450 if (fields[0]->GetBndCondExpansions().num_elements())
484 int nBndEdges, nBndEdgePts;
485 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
486 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
501 fields[var]->ExtractTracePhys(qfield, qtemp);
503 for (i = 0; i < nBndRegions; ++i)
505 nBndEdges = fields[var]->
506 GetBndCondExpansions()[i]->GetExpSize();
509 for (e = 0; e < nBndEdges ; ++e)
511 nBndEdgePts = fields[var]->
512 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
515 GetBndCondExpansions()[i]->GetPhys_Offset(e);
517 id2 = fields[0]->GetTrace()->
518 GetPhys_Offset(fields[0]->GetTraceMap()->
519 GetBndCondTraceToGlobalTraceMap(cnt++));
523 if(fields[var]->GetBndConditions()[i]->
529 &penaltyflux[id2], 1);
532 else if((fields[var]->GetBndConditions()[i])->
538 GetBndCondExpansions()[i]->
540 &penaltyflux[id2], 1);
std::string m_shockCaptureType
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
virtual void v_Diffuse(const int nConvective, 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)
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
static DiffusionSharedPtr create(std::string diffType)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
LibUtilities::SessionReaderSharedPtr m_session
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)
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)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
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, NekDouble > > m_traceNormals
virtual void v_WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayofArray
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
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_NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)