Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DiffusionLDG.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: DiffusionLDG.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: LDG diffusion class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <iostream>
38 #include <iomanip>
39 
40 namespace Nektar
41 {
42  namespace SolverUtils
43  {
45  RegisterCreatorFunction("LDG", DiffusionLDG::create);
46 
48  {
49  }
50 
54  {
55  m_session = pSession;
56 
57  m_session->LoadSolverInfo("ShockCaptureType",
58  m_shockCaptureType, "Off");
59 
60  // Setting up the normals
61  int i;
62  int nDim = pFields[0]->GetCoordim(0);
63  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
64 
66  for(i = 0; i < nDim; ++i)
67  {
68  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
69  }
70  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
71  }
72 
74  const int nConvectiveFields,
76  const Array<OneD, Array<OneD, NekDouble> > &inarray,
77  Array<OneD, Array<OneD, NekDouble> > &outarray,
78  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
79  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
80  {
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();
86 
87  Array<OneD, NekDouble> qcoeffs(nCoeffs);
88 
89  Array<OneD, Array<OneD, NekDouble> > fluxvector(nDim);
90 
94 
95  for (j = 0; j < nDim; ++j)
96  {
97  qfield[j] =
98  Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
99  qfieldStd[j] =
100  Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
101  flux[j] =
102  Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
103 
104  for (i = 0; i < nConvectiveFields; ++i)
105  {
106  qfield[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
107  qfieldStd[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
108  flux[j][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
109  }
110  }
111 
112  for (k = 0; k < nDim; ++k)
113  {
114  fluxvector[k] = Array<OneD, NekDouble>(nPts, 0.0);
115  }
116 
117  // Compute q_{\eta} and q_{\xi}
118  // Obtain numerical fluxes
119 
120  v_NumFluxforScalar(fields, inarray, flux, pFwd, pBwd);
121 
122  for (j = 0; j < nDim; ++j)
123  {
124  for (i = 0; i < nConvectiveFields; ++i)
125  {
126  fields[i]->IProductWRTDerivBase(j, inarray[i], qcoeffs);
127  Vmath::Neg (nCoeffs, qcoeffs, 1);
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]);
132  }
133  }
134  // Compute u from q_{\eta} and q_{\xi}
135  // Obtain numerical fluxes
136  v_NumFluxforVector(fields, inarray, qfield, flux[0]);
137 
139  {
140  Array<OneD, NekDouble> muvar(nPts, 0.0);
141  m_ArtificialDiffusionVector(inarray, muvar);
142 
143  int numConvFields = nConvectiveFields;
144 
145  if (m_shockCaptureType == "Smooth")
146  {
147  numConvFields = nConvectiveFields - 1;
148  }
149 
150  for (j = 0; j < nDim; ++j)
151  {
152  for (i = 0; i < numConvFields; ++i)
153  {
154  Vmath::Vmul(nPts,qfield[j][i],1,muvar,1,qfield[j][i],1);
155  }
156  }
157 
158  Array<OneD, NekDouble> FwdMuVar(nTracePts, 0.0);
159  Array<OneD, NekDouble> BwdMuVar(nTracePts, 0.0);
160 
161  fields[0]->GetFwdBwdTracePhys(muvar,FwdMuVar,BwdMuVar);
162 
163  int nBndRegions = fields[0]->GetBndCondExpansions().
164  num_elements();
165  int cnt = 0;
166 
167  for (i = 0; i < nBndRegions; ++i)
168  {
169  // Number of boundary expansion related to that region
170  int nBndEdges = fields[0]->
171  GetBndCondExpansions()[i]->GetExpSize();
172 
173  // Weakly impose boundary conditions by modifying flux
174  // values
175  for (e = 0; e < nBndEdges ; ++e)
176  {
177  nBndEdgePts = fields[0]->GetBndCondExpansions()[i]
178  ->GetExp(e)->GetTotPoints();
179 
180  int id2 = fields[0]->GetTrace()->GetPhys_Offset(
181  fields[0]->GetTraceMap()
182  ->GetBndCondTraceToGlobalTraceMap(cnt++));
183 
184  for (k = 0; k < nBndEdgePts; ++k)
185  {
186  BwdMuVar[id2+k] = 0.0;
187  }
188  }
189  }
190 
191  for(i = 0; i < numConvFields; ++i)
192  {
193  for(k = 0; k < nTracePts; ++k)
194  {
195  flux[0][i][k] =
196  0.5 * (FwdMuVar[k] + BwdMuVar[k]) * flux[0][i][k];
197  }
198  }
199  }
200 
202  Array<OneD, Array<OneD, NekDouble> > qdbase(nDim);
203 
204  for (i = 0; i < nConvectiveFields; ++i)
205  {
206  for (j = 0; j < nDim; ++j)
207  {
208  qdbase[j] = qfield[j][i];
209  }
210  fields[i]->IProductWRTDerivBase(qdbase,tmp);
211 
212  // Evaulate <\phi, \hat{F}\cdot n> - outarray[i]
213  Vmath::Neg (nCoeffs, tmp, 1);
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]);
218  }
219  }
220 
223  const Array<OneD, Array<OneD, NekDouble> > &ufield,
225  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
226  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
227  {
228  int i, j;
229  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
230  int nvariables = fields.num_elements();
231  int nDim = fields[0]->GetCoordim(0);;
232 
233  Array<OneD, NekDouble > Fwd (nTracePts);
234  Array<OneD, NekDouble > Bwd (nTracePts);
235  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
236  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
237 
238  // Get the normal velocity Vn
239  for(i = 0; i < nDim; ++i)
240  {
241  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
242  Vn, 1, Vn, 1);
243  }
244 
245  // Get the sign of (v \cdot n), v = an arbitrary vector
246  // Evaluate upwind flux:
247  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
248  for (i = 0; i < nvariables ; ++i)
249  {
250  // Compute Fwd and Bwd value of ufield of i direction
251  if (pFwd == NullNekDoubleArrayofArray ||
253  {
254  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
255  }
256  else
257  {
258  Fwd = pFwd[i];
259  Bwd = pBwd[i];
260  }
261 
262  // if Vn >= 0, flux = uFwd, i.e.,
263  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
264  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
265 
266  // else if Vn < 0, flux = uBwd, i.e.,
267  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
268  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
269  fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
270  Fwd, Bwd, fluxtemp);
271 
272  // Imposing weak boundary condition with flux
273  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
274  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
275  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
276 
277  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
278  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
279  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
280  if(fields[0]->GetBndCondExpansions().num_elements())
281  {
282  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
283  }
284 
285  for (j = 0; j < nDim; ++j)
286  {
287  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
288  // i.e,
289  // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
290  // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
291 
292  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
293  // i.e,
294  // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
295  // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
296 
297  Vmath::Vmul(nTracePts,
298  m_traceNormals[j], 1,
299  fluxtemp, 1,
300  uflux[j][i], 1);
301  }
302  }
303  }
304 
305 
306 
309  const int var,
310  const Array<OneD, const NekDouble> &ufield,
311  Array<OneD, NekDouble> &penaltyflux)
312  {
313  int i, e, id1, id2;
314 
315  // Number of boundary regions
316  int nBndEdgePts, nBndEdges;
317  int cnt = 0;
318  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
319  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
320  Array<OneD, NekDouble > uplus(nTracePts);
321 
322  fields[var]->ExtractTracePhys(ufield, uplus);
323  for (i = 0; i < nBndRegions; ++i)
324  {
325  // Number of boundary expansion related to that region
326  nBndEdges = fields[var]->
327  GetBndCondExpansions()[i]->GetExpSize();
328 
329  // Weakly impose boundary conditions by modifying flux values
330  for (e = 0; e < nBndEdges ; ++e)
331  {
332  nBndEdgePts = fields[var]->
333  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
334 
335  id1 = fields[var]->
336  GetBndCondExpansions()[i]->GetPhys_Offset(e);
337 
338  id2 = fields[0]->GetTrace()->
339  GetPhys_Offset(fields[0]->GetTraceMap()->
340  GetBndCondTraceToGlobalTraceMap(cnt++));
341 
342  // For Dirichlet boundary condition: uflux = g_D
343  if (fields[var]->GetBndConditions()[i]->
344  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
345  {
346  Vmath::Vcopy(nBndEdgePts,
347  &(fields[var]->
348  GetBndCondExpansions()[i]->
349  GetPhys())[id1], 1,
350  &penaltyflux[id2], 1);
351  }
352  // For Neumann boundary condition: uflux = u+
353  else if ((fields[var]->GetBndConditions()[i])->
354  GetBoundaryConditionType() == SpatialDomains::eNeumann)
355  {
356  Vmath::Vcopy(nBndEdgePts,
357  &uplus[id2], 1,
358  &penaltyflux[id2], 1);
359  }
360  }
361  }
362  }
363 
364 
365 
368  const Array<OneD, Array<OneD, NekDouble> > &ufield,
371  {
372  int i, j;
373  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
374  int nvariables = fields.num_elements();
375  int nDim = qfield.num_elements();
376 
377  NekDouble C11 = 0.0;
378  Array<OneD, NekDouble > Fwd(nTracePts);
379  Array<OneD, NekDouble > Bwd(nTracePts);
380  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
381  Array<OneD, NekDouble > qFwd (nTracePts);
382  Array<OneD, NekDouble > qBwd (nTracePts);
383  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
384  Array<OneD, NekDouble > uterm(nTracePts);
385  /*
386  // Setting up the normals
387  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDim);
388  for(i = 0; i < nDim; ++i)
389  {
390  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
391  }
392  fields[0]->GetTrace()->GetNormals(m_traceNormals);
393  */
394 
395  // Get the normal velocity Vn
396  for(i = 0; i < nDim; ++i)
397  {
398  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
399  Vn, 1, Vn, 1);
400  }
401  // Evaulate upwind flux:
402  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
403  for (i = 0; i < nvariables; ++i)
404  {
405  // Generate Stability term = - C11 ( u- - u+ )
406  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
407 
408  Vmath::Vsub(nTracePts,
409  Fwd, 1, Bwd, 1,
410  uterm, 1);
411 
412  Vmath::Smul(nTracePts,
413  -1.0 * C11, uterm, 1,
414  uterm, 1);
415 
416  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
417  for (j = 0; j < nDim; ++j)
418  {
419  // Compute Fwd and Bwd value of ufield of jth direction
420  fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
421 
422  // if Vn >= 0, flux = uFwd, i.e.,
423  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
424  // qflux = qBwd = q+
425  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
426  // qflux = qBwd = q-
427 
428  // else if Vn < 0, flux = uBwd, i.e.,
429  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
430  // qflux = qFwd = q-
431  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
432  // qflux = qFwd = q+
433 
434  fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
435  qBwd, qFwd,
436  qfluxtemp);
437 
438  Vmath::Vmul(nTracePts,
439  m_traceNormals[j], 1,
440  qfluxtemp, 1,
441  qfluxtemp, 1);
442 
443  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
444  Vmath::Vadd(nTracePts,
445  uterm, 1,
446  qfluxtemp, 1,
447  qfluxtemp, 1);
448 
449  // Imposing weak boundary condition with flux
450  if (fields[0]->GetBndCondExpansions().num_elements())
451  {
452  v_WeakPenaltyforVector(fields, i, j,
453  qfield[j][i],
454  qfluxtemp, C11);
455  }
456 
457  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
458  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
459  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
460  Vmath::Vadd(nTracePts,
461  qfluxtemp, 1,
462  qflux[i], 1,
463  qflux[i], 1);
464  }
465  }
466  }
467 
468 
469 
470  /**
471  * Diffusion: Imposing weak boundary condition for q with flux
472  * uflux = g_D on Dirichlet boundary condition
473  * uflux = u_Fwd on Neumann boundary condition
474  */
477  const int var,
478  const int dir,
479  const Array<OneD, const NekDouble> &qfield,
480  Array<OneD, NekDouble> &penaltyflux,
481  NekDouble C11)
482  {
483  int i, e, id1, id2;
484  int nBndEdges, nBndEdgePts;
485  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
486  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
487 
488  Array<OneD, NekDouble > qtemp(nTracePts);
489  int cnt = 0;
490 
491  /*
492  // Setting up the normals
493  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDim);
494  for(i = 0; i < nDim; ++i)
495  {
496  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
497  }
498  fields[0]->GetTrace()->GetNormals(m_traceNormals);
499  */
500 
501  fields[var]->ExtractTracePhys(qfield, qtemp);
502 
503  for (i = 0; i < nBndRegions; ++i)
504  {
505  nBndEdges = fields[var]->
506  GetBndCondExpansions()[i]->GetExpSize();
507 
508  // Weakly impose boundary conditions by modifying flux values
509  for (e = 0; e < nBndEdges ; ++e)
510  {
511  nBndEdgePts = fields[var]->
512  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
513 
514  id1 = fields[var]->
515  GetBndCondExpansions()[i]->GetPhys_Offset(e);
516 
517  id2 = fields[0]->GetTrace()->
518  GetPhys_Offset(fields[0]->GetTraceMap()->
519  GetBndCondTraceToGlobalTraceMap(cnt++));
520 
521  // For Dirichlet boundary condition:
522  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
523  if(fields[var]->GetBndConditions()[i]->
524  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
525  {
526  Vmath::Vmul(nBndEdgePts,
527  &m_traceNormals[dir][id2], 1,
528  &qtemp[id2], 1,
529  &penaltyflux[id2], 1);
530  }
531  // For Neumann boundary condition: qflux = g_N
532  else if((fields[var]->GetBndConditions()[i])->
533  GetBoundaryConditionType() == SpatialDomains::eNeumann)
534  {
535  Vmath::Vmul(nBndEdgePts,
536  &m_traceNormals[dir][id2], 1,
537  &(fields[var]->
538  GetBndCondExpansions()[i]->
539  GetPhys())[id1], 1,
540  &penaltyflux[id2], 1);
541  }
542  }
543  }
544  }
545 
546  }
547 }
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
Definition: Diffusion.h:131
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()
Definition: Diffusion.cpp:42
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
Definition: Vmath.cpp:485
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLDG.h:48
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLDG.h:61
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.
Definition: Vmath.cpp:213
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
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.
Definition: Vmath.cpp:343
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLDG.h:60
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)
Definition: Vmath.cpp:1061
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.
Definition: Vmath.cpp:299
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.
Definition: Vmath.cpp:183
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)