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,
79 int nBndEdgePts, i, j, k, e;
80 int nDim = fields[0]->GetCoordim(0);
81 int nPts = fields[0]->GetTotPoints();
82 int nCoeffs = fields[0]->GetNcoeffs();
83 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
94 for (j = 0; j < nDim; ++j)
103 for (i = 0; i < nConvectiveFields; ++i)
111 for (k = 0; k < nDim; ++k)
121 for (j = 0; j < nDim; ++j)
123 for (i = 0; i < nConvectiveFields; ++i)
125 fields[i]->IProductWRTDerivBase(j, inarray[i], qcoeffs);
127 fields[i]->AddTraceIntegral (flux[j][i], qcoeffs);
128 fields[i]->SetPhysState (
false);
129 fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);
130 fields[i]->BwdTrans (qcoeffs, qfield[j][i]);
142 int numConvFields = nConvectiveFields;
146 numConvFields = nConvectiveFields - 1;
149 for (j = 0; j < nDim; ++j)
151 for (i = 0; i < numConvFields; ++i)
153 Vmath::Vmul(nPts,qfield[j][i],1,muvar,1,qfield[j][i],1);
160 fields[0]->GetFwdBwdTracePhys(muvar,FwdMuVar,BwdMuVar);
162 int nBndRegions = fields[0]->GetBndCondExpansions().
166 for (i = 0; i < nBndRegions; ++i)
169 int nBndEdges = fields[0]->
170 GetBndCondExpansions()[i]->GetExpSize();
174 for (e = 0; e < nBndEdges ; ++e)
176 nBndEdgePts = fields[0]->GetBndCondExpansions()[i]
177 ->GetExp(e)->GetTotPoints();
179 int id2 = fields[0]->GetTrace()->GetPhys_Offset(
180 fields[0]->GetTraceMap()
181 ->GetBndCondTraceToGlobalTraceMap(cnt++));
183 for (k = 0; k < nBndEdgePts; ++k)
185 BwdMuVar[id2+k] = 0.0;
190 for(i = 0; i < numConvFields; ++i)
192 for(k = 0; k < nTracePts; ++k)
195 0.5 * (FwdMuVar[k] + BwdMuVar[k]) * flux[0][i][k];
203 for (i = 0; i < nConvectiveFields; ++i)
205 for (j = 0; j < nDim; ++j)
207 qdbase[j] = qfield[j][i];
209 fields[i]->IProductWRTDerivBase(qdbase,tmp);
213 fields[i]->AddTraceIntegral (flux[0][i], tmp);
214 fields[i]->SetPhysState (
false);
215 fields[i]->MultiplyByElmtInvMass(tmp, tmp);
216 fields[i]->BwdTrans (tmp, outarray[i]);
226 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
227 int nvariables = fields.num_elements();
228 int nDim = uflux.num_elements();
236 for(i = 0; i < nDim; ++i)
245 for (j = 0; j < nDim; ++j)
247 for (i = 0; i < nvariables ; ++i)
250 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
260 fields[i]->GetTrace()->Upwind(Vn,
272 if(fields[0]->GetBndCondExpansions().num_elements())
306 int nBndEdgePts, nBndEdges;
308 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
309 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
312 fields[var]->ExtractTracePhys(ufield, uplus);
313 for (i = 0; i < nBndRegions; ++i)
316 nBndEdges = fields[var]->
317 GetBndCondExpansions()[i]->GetExpSize();
320 for (e = 0; e < nBndEdges ; ++e)
322 nBndEdgePts = fields[var]->
323 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
326 GetBndCondExpansions()[i]->GetPhys_Offset(e);
328 id2 = fields[0]->GetTrace()->
329 GetPhys_Offset(fields[0]->GetTraceMap()->
330 GetBndCondTraceToGlobalTraceMap(cnt++));
333 if (fields[var]->GetBndConditions()[i]->
338 GetBndCondExpansions()[i]->
340 &penaltyflux[id2], 1);
343 else if ((fields[var]->GetBndConditions()[i])->
348 &penaltyflux[id2], 1);
363 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
364 int nvariables = fields.num_elements();
365 int nDim = qfield.num_elements();
386 for(i = 0; i < nDim; ++i)
393 for (i = 0; i < nvariables; ++i)
396 for (j = 0; j < nDim; ++j)
399 fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
413 fields[i]->GetTrace()->Upwind(Vn,
423 fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
430 -1.0 * C11, uterm, 1,
440 if (fields[0]->GetBndCondExpansions().num_elements())
474 int nBndEdges, nBndEdgePts;
475 int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
476 int nTracePts = fields[0]->GetTrace()->GetTotPoints();
492 fields[var]->ExtractTracePhys(qfield, qtemp);
494 for (i = 0; i < nBndRegions; ++i)
496 nBndEdges = fields[var]->
497 GetBndCondExpansions()[i]->GetExpSize();
500 for (e = 0; e < nBndEdges ; ++e)
502 nBndEdgePts = fields[var]->
503 GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
506 GetBndCondExpansions()[i]->GetPhys_Offset(e);
508 id2 = fields[0]->GetTrace()->
509 GetPhys_Offset(fields[0]->GetTraceMap()->
510 GetBndCondTraceToGlobalTraceMap(cnt++));
514 if(fields[var]->GetBndConditions()[i]->
520 &penaltyflux[id2], 1);
523 else if((fields[var]->GetBndConditions()[i])->
529 GetBndCondExpansions()[i]->
531 &penaltyflux[id2], 1);
std::string m_shockCaptureType
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
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
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)
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.
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)
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)
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.