Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AdvectionFR.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: AdvectionFR.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: FR advection class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <SolverUtils/Advection/AdvectionFR.h>
37 #include <LibUtilities/Foundations/ManagerAccess.h>
38 #include <LibUtilities/Foundations/Basis.h>
39 #include <LibUtilities/Foundations/Points.h>
40 #include <LibUtilities/Polylib/Polylib.h>
41 #include <StdRegions/StdSegExp.h>
42 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
43 #include <MultiRegions/DisContField1D.h>
44 #include <boost/math/special_functions/gamma.hpp>
45 #include <LocalRegions/Expansion1D.h>
46 #include <LocalRegions/Expansion2D.h>
47 
48 #include <iostream>
49 #include <iomanip>
50 
51 namespace Nektar
52 {
53  namespace SolverUtils
54  {
55  std::string AdvectionFR::type[] = {
56  GetAdvectionFactory().RegisterCreatorFunction(
57  "FRDG", AdvectionFR::create),
58  GetAdvectionFactory().RegisterCreatorFunction(
59  "FRSD", AdvectionFR::create),
60  GetAdvectionFactory().RegisterCreatorFunction(
61  "FRHU", AdvectionFR::create),
62  GetAdvectionFactory().RegisterCreatorFunction(
63  "FRcmin", AdvectionFR::create),
64  GetAdvectionFactory().RegisterCreatorFunction(
65  "FRcinf", AdvectionFR::create)};
66 
67  /**
68  * @brief AdvectionFR uses the Flux Reconstruction (FR) approach to
69  * compute the advection term. The implementation is only for segments,
70  * quadrilaterals and hexahedra at the moment.
71  *
72  * \todo Extension to triangles, tetrahedra and other shapes.
73  * (Long term objective)
74  */
75  AdvectionFR::AdvectionFR(std::string advType):m_advType(advType)
76  {
77  }
78 
79  /**
80  * @brief Initiliase AdvectionFR objects and store them before starting
81  * the time-stepping.
82  *
83  * This routine calls the virtual functions #v_SetupMetrics and
84  * #v_SetupCFunctions to initialise the objects needed by AdvectionFR.
85  *
86  * @param pSession Pointer to session reader.
87  * @param pFields Pointer to fields.
88  */
89  void AdvectionFR::v_InitObject(
90  LibUtilities::SessionReaderSharedPtr pSession,
91  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
92  {
93  Advection::v_InitObject(pSession, pFields);
94  v_SetupMetrics (pSession, pFields);
95  v_SetupCFunctions (pSession, pFields);
96  }
97 
98  /**
99  * @brief Setup the metric terms to compute the contravariant
100  * fluxes. (i.e. this special metric terms transform the fluxes
101  * at the interfaces of each element from the physical space to
102  * the standard space).
103  *
104  * This routine calls the function #GetEdgeQFactors to compute and
105  * store the metric factors following an anticlockwise conventions
106  * along the edges/faces of each element. Note: for 1D problem
107  * the transformation is not needed.
108  *
109  * @param pSession Pointer to session reader.
110  * @param pFields Pointer to fields.
111  *
112  * \todo Add the metric terms for 3D Hexahedra.
113  */
114  void AdvectionFR::v_SetupMetrics(
115  LibUtilities::SessionReaderSharedPtr pSession,
116  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
117  {
118  int i, n;
119  int nquad0, nquad1;
120  int phys_offset;
121  int nLocalSolutionPts;
122  int nElements = pFields[0]->GetExpSize();
123  int nDimensions = pFields[0]->GetCoordim(0);
124  int nSolutionPts = pFields[0]->GetTotPoints();
125  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
126 
127  Array<TwoD, const NekDouble> gmat;
128  Array<OneD, const NekDouble> jac;
129 
130  m_jac = Array<OneD, NekDouble>(nSolutionPts);
131 
132  Array<OneD, NekDouble> auxArray1;
133  Array<OneD, LibUtilities::BasisSharedPtr> base;
134  LibUtilities::PointsKeyVector ptsKeys;
135 
136  switch (nDimensions)
137  {
138  case 1:
139  {
140  for (n = 0; n < nElements; ++n)
141  {
142  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
143  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
144  phys_offset = pFields[0]->GetPhys_Offset(n);
145  jac = pFields[0]->GetExp(n)
146  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
147  ->GetMetricInfo()->GetJac(ptsKeys);
148  for (i = 0; i < nLocalSolutionPts; ++i)
149  {
150  m_jac[i+phys_offset] = jac[0];
151  }
152  }
153  break;
154  }
155  case 2:
156  {
157  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
158  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
159  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
160  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
161  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
162 
163  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
164  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
165  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
166  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
167 
168  LibUtilities::PointsKeyVector ptsKeys;
169 
170  for (n = 0; n < nElements; ++n)
171  {
172  base = pFields[0]->GetExp(n)->GetBase();
173  nquad0 = base[0]->GetNumPoints();
174  nquad1 = base[1]->GetNumPoints();
175 
176  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
177  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
178  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
179  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
180 
181  // Extract the Q factors at each edge point
182  pFields[0]->GetExp(n)->GetEdgeQFactors(
183  0, auxArray1 = m_Q2D_e0[n]);
184  pFields[0]->GetExp(n)->GetEdgeQFactors(
185  1, auxArray1 = m_Q2D_e1[n]);
186  pFields[0]->GetExp(n)->GetEdgeQFactors(
187  2, auxArray1 = m_Q2D_e2[n]);
188  pFields[0]->GetExp(n)->GetEdgeQFactors(
189  3, auxArray1 = m_Q2D_e3[n]);
190 
191  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
192  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
193  phys_offset = pFields[0]->GetPhys_Offset(n);
194 
195  jac = pFields[0]->GetExp(n)
196  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
197  ->GetMetricInfo()->GetJac(ptsKeys);
198  gmat = pFields[0]->GetExp(n)
199  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
200  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
201 
202  if (pFields[0]->GetExp(n)
203  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
204  ->GetMetricInfo()->GetGtype()
205  == SpatialDomains::eDeformed)
206  {
207  for (i = 0; i < nLocalSolutionPts; ++i)
208  {
209  m_jac[i+phys_offset] = jac[i];
210  m_gmat[0][i+phys_offset] = gmat[0][i];
211  m_gmat[1][i+phys_offset] = gmat[1][i];
212  m_gmat[2][i+phys_offset] = gmat[2][i];
213  m_gmat[3][i+phys_offset] = gmat[3][i];
214  }
215  }
216  else
217  {
218  for (i = 0; i < nLocalSolutionPts; ++i)
219  {
220  m_jac[i+phys_offset] = jac[0];
221  m_gmat[0][i+phys_offset] = gmat[0][0];
222  m_gmat[1][i+phys_offset] = gmat[1][0];
223  m_gmat[2][i+phys_offset] = gmat[2][0];
224  m_gmat[3][i+phys_offset] = gmat[3][0];
225  }
226  }
227  }
228 
229  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(
230  nDimensions);
231  for(i = 0; i < nDimensions; ++i)
232  {
233  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
234  }
235  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
236 
237  break;
238  }
239  case 3:
240  {
241  ASSERTL0(false,"3DFR Metric terms not implemented yet");
242  break;
243  }
244  default:
245  {
246  ASSERTL0(false, "Expansion dimension not recognised");
247  break;
248  }
249  }
250  }
251 
252  /**
253  * @brief Setup the derivatives of the correction functions. For more
254  * details see J Sci Comput (2011) 47: 50–72
255  *
256  * This routine calls 3 different bases:
257  * #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme,
258  * #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme,
259  * #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme.
260  * The values of the derivatives of the correction function are then
261  * stored into global variables and reused into the virtual functions
262  * #v_DivCFlux_1D, #v_DivCFlux_2D, #v_DivCFlux_3D to compute the
263  * the divergence of the correction flux for 1D, 2D or 3D problems.
264  *
265  * @param pSession Pointer to session reader.
266  * @param pFields Pointer to fields.
267  */
268  void AdvectionFR::v_SetupCFunctions(
269  LibUtilities::SessionReaderSharedPtr pSession,
270  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
271  {
272  int i, n;
273  NekDouble c0, c1, c2;
274  int nquad0, nquad1, nquad2;
275  int nmodes0, nmodes1, nmodes2;
276  Array<OneD, LibUtilities::BasisSharedPtr> base;
277 
278  int nElements = pFields[0]->GetExpSize();
279  int nDimensions = pFields[0]->GetCoordim(0);
280 
281  switch (nDimensions)
282  {
283  case 1:
284  {
285  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
286  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
287 
288  for (n = 0; n < nElements; ++n)
289  {
290  base = pFields[0]->GetExp(n)->GetBase();
291  nquad0 = base[0]->GetNumPoints();
292  nmodes0 = base[0]->GetNumModes();
293  Array<OneD, const NekDouble> z0;
294  Array<OneD, const NekDouble> w0;
295 
296  base[0]->GetZW(z0, w0);
297 
298  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
299  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
300 
301  // Auxiliary vectors to build up the auxiliary
302  // derivatives of the Legendre polynomials
303  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
304  Array<OneD,NekDouble> dLpp0(nquad0, 0.0);
305  Array<OneD,NekDouble> dLpm0(nquad0, 0.0);
306 
307  // Degree of the correction functions
308  int p0 = nmodes0 - 1;
309 
310  // Function sign to compute left correction function
311  NekDouble sign0 = pow(-1.0, p0);
312 
313  // Factorial factor to build the scheme
314  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
315  / (pow(2.0, p0)
316  * boost::math::tgamma(p0 + 1)
317  * boost::math::tgamma(p0 + 1));
318 
319  // Scalar parameter which recovers the FR schemes
320  if (m_advType == "FRDG")
321  {
322  c0 = 0.0;
323  }
324  else if (m_advType == "FRSD")
325  {
326  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
327  * (ap0 * boost::math::tgamma(p0 + 1))
328  * (ap0 * boost::math::tgamma(p0 + 1)));
329  }
330  else if (m_advType == "FRHU")
331  {
332  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
333  * (ap0 * boost::math::tgamma(p0 + 1))
334  * (ap0 * boost::math::tgamma(p0 + 1)));
335  }
336  else if (m_advType == "FRcmin")
337  {
338  c0 = -2.0 / ((2.0 * p0 + 1.0)
339  * (ap0 * boost::math::tgamma(p0 + 1))
340  * (ap0 * boost::math::tgamma(p0 + 1)));
341  }
342  else if (m_advType == "FRcinf")
343  {
344  c0 = 10000000000000000.0;
345  }
346 
347  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
348  * (ap0 * boost::math::tgamma(p0 + 1))
349  * (ap0 * boost::math::tgamma(p0 + 1));
350 
351  NekDouble overeta0 = 1.0 / (1.0 + etap0);
352 
353  // Derivative of the Legendre polynomials
354  // dLp = derivative of the Legendre polynomial of order p
355  // dLpp = derivative of the Legendre polynomial of order p+1
356  // dLpm = derivative of the Legendre polynomial of order p-1
357  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
358  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
359  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
360 
361  // Building the DG_c_Left
362  for(i = 0; i < nquad0; ++i)
363  {
364  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
365  m_dGL_xi1[n][i] += dLpp0[i];
366  m_dGL_xi1[n][i] *= overeta0;
367  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
368  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
369  }
370 
371  // Building the DG_c_Right
372  for(i = 0; i < nquad0; ++i)
373  {
374  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
375  m_dGR_xi1[n][i] += dLpp0[i];
376  m_dGR_xi1[n][i] *= overeta0;
377  m_dGR_xi1[n][i] += dLp0[i];
378  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
379  }
380  }
381  break;
382  }
383  case 2:
384  {
385  LibUtilities::BasisSharedPtr BasisFR_Left0;
386  LibUtilities::BasisSharedPtr BasisFR_Right0;
387  LibUtilities::BasisSharedPtr BasisFR_Left1;
388  LibUtilities::BasisSharedPtr BasisFR_Right1;
389 
390  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
391  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
392  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
393  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
394 
395  for (n = 0; n < nElements; ++n)
396  {
397  base = pFields[0]->GetExp(n)->GetBase();
398  nquad0 = base[0]->GetNumPoints();
399  nquad1 = base[1]->GetNumPoints();
400  nmodes0 = base[0]->GetNumModes();
401  nmodes1 = base[1]->GetNumModes();
402 
403  Array<OneD, const NekDouble> z0;
404  Array<OneD, const NekDouble> w0;
405  Array<OneD, const NekDouble> z1;
406  Array<OneD, const NekDouble> w1;
407 
408  base[0]->GetZW(z0, w0);
409  base[1]->GetZW(z1, w1);
410 
411  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
412  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
413  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
414  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
415 
416  // Auxiliary vectors to build up the auxiliary
417  // derivatives of the Legendre polynomials
418  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
419  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
420  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
421  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
422  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
423  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
424 
425  // Degree of the correction functions
426  int p0 = nmodes0 - 1;
427  int p1 = nmodes1 - 1;
428 
429  // Function sign to compute left correction function
430  NekDouble sign0 = pow(-1.0, p0);
431  NekDouble sign1 = pow(-1.0, p1);
432 
433  // Factorial factor to build the scheme
434  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
435  / (pow(2.0, p0)
436  * boost::math::tgamma(p0 + 1)
437  * boost::math::tgamma(p0 + 1));
438 
439  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
440  / (pow(2.0, p1)
441  * boost::math::tgamma(p1 + 1)
442  * boost::math::tgamma(p1 + 1));
443 
444  // Scalar parameter which recovers the FR schemes
445  if (m_advType == "FRDG")
446  {
447  c0 = 0.0;
448  c1 = 0.0;
449  }
450  else if (m_advType == "FRSD")
451  {
452  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
453  * (ap0 * boost::math::tgamma(p0 + 1))
454  * (ap0 * boost::math::tgamma(p0 + 1)));
455 
456  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
457  * (ap1 * boost::math::tgamma(p1 + 1))
458  * (ap1 * boost::math::tgamma(p1 + 1)));
459  }
460  else if (m_advType == "FRHU")
461  {
462  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
463  * (ap0 * boost::math::tgamma(p0 + 1))
464  * (ap0 * boost::math::tgamma(p0 + 1)));
465 
466  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
467  * (ap1 * boost::math::tgamma(p1 + 1))
468  * (ap1 * boost::math::tgamma(p1 + 1)));
469  }
470  else if (m_advType == "FRcmin")
471  {
472  c0 = -2.0 / ((2.0 * p0 + 1.0)
473  * (ap0 * boost::math::tgamma(p0 + 1))
474  * (ap0 * boost::math::tgamma(p0 + 1)));
475 
476  c1 = -2.0 / ((2.0 * p1 + 1.0)
477  * (ap1 * boost::math::tgamma(p1 + 1))
478  * (ap1 * boost::math::tgamma(p1 + 1)));
479  }
480  else if (m_advType == "FRcinf")
481  {
482  c0 = 10000000000000000.0;
483  c1 = 10000000000000000.0;
484  }
485 
486  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
487  * (ap0 * boost::math::tgamma(p0 + 1))
488  * (ap0 * boost::math::tgamma(p0 + 1));
489 
490  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
491  * (ap1 * boost::math::tgamma(p1 + 1))
492  * (ap1 * boost::math::tgamma(p1 + 1));
493 
494  NekDouble overeta0 = 1.0 / (1.0 + etap0);
495  NekDouble overeta1 = 1.0 / (1.0 + etap1);
496 
497  // Derivative of the Legendre polynomials
498  // dLp = derivative of the Legendre polynomial of order p
499  // dLpp = derivative of the Legendre polynomial of order p+1
500  // dLpm = derivative of the Legendre polynomial of order p-1
501  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
502  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
503  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
504 
505  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
506  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
507  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
508 
509  // Building the DG_c_Left
510  for(i = 0; i < nquad0; ++i)
511  {
512  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
513  m_dGL_xi1[n][i] += dLpp0[i];
514  m_dGL_xi1[n][i] *= overeta0;
515  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
516  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
517  }
518 
519  // Building the DG_c_Left
520  for(i = 0; i < nquad1; ++i)
521  {
522  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
523  m_dGL_xi2[n][i] += dLpp1[i];
524  m_dGL_xi2[n][i] *= overeta1;
525  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
526  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
527  }
528 
529  // Building the DG_c_Right
530  for(i = 0; i < nquad0; ++i)
531  {
532  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
533  m_dGR_xi1[n][i] += dLpp0[i];
534  m_dGR_xi1[n][i] *= overeta0;
535  m_dGR_xi1[n][i] += dLp0[i];
536  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
537  }
538 
539  // Building the DG_c_Right
540  for(i = 0; i < nquad1; ++i)
541  {
542  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
543  m_dGR_xi2[n][i] += dLpp1[i];
544  m_dGR_xi2[n][i] *= overeta1;
545  m_dGR_xi2[n][i] += dLp1[i];
546  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
547  }
548  }
549  break;
550  }
551  case 3:
552  {
553  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
554  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
555  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
556  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
557  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
558  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
559 
560  for (n = 0; n < nElements; ++n)
561  {
562  base = pFields[0]->GetExp(n)->GetBase();
563  nquad0 = base[0]->GetNumPoints();
564  nquad1 = base[1]->GetNumPoints();
565  nquad2 = base[2]->GetNumPoints();
566  nmodes0 = base[0]->GetNumModes();
567  nmodes1 = base[1]->GetNumModes();
568  nmodes2 = base[2]->GetNumModes();
569 
570  Array<OneD, const NekDouble> z0;
571  Array<OneD, const NekDouble> w0;
572  Array<OneD, const NekDouble> z1;
573  Array<OneD, const NekDouble> w1;
574  Array<OneD, const NekDouble> z2;
575  Array<OneD, const NekDouble> w2;
576 
577  base[0]->GetZW(z0, w0);
578  base[1]->GetZW(z1, w1);
579  base[1]->GetZW(z2, w2);
580 
581  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
582  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
583  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
584  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
585  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
586  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
587 
588  // Auxiliary vectors to build up the auxiliary
589  // derivatives of the Legendre polynomials
590  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
591  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
592  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
593  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
594  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
595  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
596  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
597  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
598  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
599 
600  // Degree of the correction functions
601  int p0 = nmodes0 - 1;
602  int p1 = nmodes1 - 1;
603  int p2 = nmodes2 - 1;
604 
605  // Function sign to compute left correction function
606  NekDouble sign0 = pow(-1.0, p0);
607  NekDouble sign1 = pow(-1.0, p1);
608 
609  // Factorial factor to build the scheme
610  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
611  / (pow(2.0, p0)
612  * boost::math::tgamma(p0 + 1)
613  * boost::math::tgamma(p0 + 1));
614 
615  // Factorial factor to build the scheme
616  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
617  / (pow(2.0, p1)
618  * boost::math::tgamma(p1 + 1)
619  * boost::math::tgamma(p1 + 1));
620 
621  // Factorial factor to build the scheme
622  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
623  / (pow(2.0, p2)
624  * boost::math::tgamma(p2 + 1)
625  * boost::math::tgamma(p2 + 1));
626 
627  // Scalar parameter which recovers the FR schemes
628  if (m_advType == "FRDG")
629  {
630  c0 = 0.0;
631  c1 = 0.0;
632  c2 = 0.0;
633  }
634  else if (m_advType == "FRSD")
635  {
636  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
637  * (ap0 * boost::math::tgamma(p0 + 1))
638  * (ap0 * boost::math::tgamma(p0 + 1)));
639 
640  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
641  * (ap1 * boost::math::tgamma(p1 + 1))
642  * (ap1 * boost::math::tgamma(p1 + 1)));
643 
644  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
645  * (ap2 * boost::math::tgamma(p2 + 1))
646  * (ap2 * boost::math::tgamma(p2 + 1)));
647  }
648  else if (m_advType == "FRHU")
649  {
650  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
651  * (ap0 * boost::math::tgamma(p0 + 1))
652  * (ap0 * boost::math::tgamma(p0 + 1)));
653 
654  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
655  * (ap1 * boost::math::tgamma(p1 + 1))
656  * (ap1 * boost::math::tgamma(p1 + 1)));
657 
658  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
659  * (ap2 * boost::math::tgamma(p2 + 1))
660  * (ap2 * boost::math::tgamma(p2 + 1)));
661  }
662  else if (m_advType == "FRcmin")
663  {
664  c0 = -2.0 / ((2.0 * p0 + 1.0)
665  * (ap0 * boost::math::tgamma(p0 + 1))
666  * (ap0 * boost::math::tgamma(p0 + 1)));
667 
668  c1 = -2.0 / ((2.0 * p1 + 1.0)
669  * (ap1 * boost::math::tgamma(p1 + 1))
670  * (ap1 * boost::math::tgamma(p1 + 1)));
671 
672  c2 = -2.0 / ((2.0 * p2 + 1.0)
673  * (ap2 * boost::math::tgamma(p2 + 1))
674  * (ap2 * boost::math::tgamma(p2 + 1)));
675  }
676  else if (m_advType == "FRcinf")
677  {
678  c0 = 10000000000000000.0;
679  c1 = 10000000000000000.0;
680  c2 = 10000000000000000.0;
681  }
682 
683  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
684  * (ap0 * boost::math::tgamma(p0 + 1))
685  * (ap0 * boost::math::tgamma(p0 + 1));
686 
687  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
688  * (ap1 * boost::math::tgamma(p1 + 1))
689  * (ap1 * boost::math::tgamma(p1 + 1));
690 
691  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
692  * (ap2 * boost::math::tgamma(p2 + 1))
693  * (ap2 * boost::math::tgamma(p2 + 1));
694 
695  NekDouble overeta0 = 1.0 / (1.0 + etap0);
696  NekDouble overeta1 = 1.0 / (1.0 + etap1);
697  NekDouble overeta2 = 1.0 / (1.0 + etap2);
698 
699  // Derivative of the Legendre polynomials
700  // dLp = derivative of the Legendre polynomial of order p
701  // dLpp = derivative of the Legendre polynomial of order p+1
702  // dLpm = derivative of the Legendre polynomial of order p-1
703  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
704  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
705  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
706 
707  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
708  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
709  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
710 
711  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
712  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
713  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
714 
715  // Building the DG_c_Left
716  for(i = 0; i < nquad0; ++i)
717  {
718  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
719  m_dGL_xi1[n][i] += dLpp0[i];
720  m_dGL_xi1[n][i] *= overeta0;
721  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
722  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
723  }
724 
725  // Building the DG_c_Left
726  for(i = 0; i < nquad1; ++i)
727  {
728  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
729  m_dGL_xi2[n][i] += dLpp1[i];
730  m_dGL_xi2[n][i] *= overeta1;
731  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
732  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
733  }
734 
735  // Building the DG_c_Left
736  for(i = 0; i < nquad2; ++i)
737  {
738  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
739  m_dGL_xi3[n][i] += dLpp2[i];
740  m_dGL_xi3[n][i] *= overeta2;
741  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
742  m_dGL_xi3[n][i] = 0.5 * sign1 * m_dGL_xi3[n][i];
743  }
744 
745  // Building the DG_c_Right
746  for(i = 0; i < nquad0; ++i)
747  {
748  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
749  m_dGR_xi1[n][i] += dLpp0[i];
750  m_dGR_xi1[n][i] *= overeta0;
751  m_dGR_xi1[n][i] += dLp0[i];
752  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
753  }
754 
755  // Building the DG_c_Right
756  for(i = 0; i < nquad1; ++i)
757  {
758  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
759  m_dGR_xi2[n][i] += dLpp1[i];
760  m_dGR_xi2[n][i] *= overeta1;
761  m_dGR_xi2[n][i] += dLp1[i];
762  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
763  }
764 
765  // Building the DG_c_Right
766  for(i = 0; i < nquad2; ++i)
767  {
768  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
769  m_dGR_xi3[n][i] += dLpp2[i];
770  m_dGR_xi3[n][i] *= overeta2;
771  m_dGR_xi3[n][i] += dLp2[i];
772  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
773  }
774  }
775  break;
776  }
777  default:
778  {
779  ASSERTL0(false,"Expansion dimension not recognised");
780  break;
781  }
782  }
783  }
784 
785  /**
786  * @brief Compute the advection term at each time-step using the Flux
787  * Reconstruction approach (FR).
788  *
789  * @param nConvectiveFields Number of fields.
790  * @param fields Pointer to fields.
791  * @param advVel Advection velocities.
792  * @param inarray Solution at the previous time-step.
793  * @param outarray Advection term to be passed at the
794  * time integration class.
795  *
796  */
797  void AdvectionFR::v_Advect(
798  const int nConvectiveFields,
799  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
800  const Array<OneD, Array<OneD, NekDouble> > &advVel,
801  const Array<OneD, Array<OneD, NekDouble> > &inarray,
802  Array<OneD, Array<OneD, NekDouble> > &outarray)
803  {
804  int i, j, n;
805  int phys_offset;
806 
807  Array<OneD, NekDouble> auxArray1, auxArray2;
808 
809  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
810  Basis = fields[0]->GetExp(0)->GetBase();
811 
812  int nElements = fields[0]->GetExpSize();
813  int nDimensions = fields[0]->GetCoordim(0);
814  int nSolutionPts = fields[0]->GetTotPoints();
815  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
816  int nCoeffs = fields[0]->GetNcoeffs();
817 
818  // Storage for flux vector.
819  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > fluxvector
820  (nConvectiveFields);
821  // Outarray for Galerkin projection in case of primitive dealising
822  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff
823  (nConvectiveFields);
824  // Store forwards/backwards space along trace space
825  Array<OneD, Array<OneD, NekDouble> > Fwd (nConvectiveFields);
826  Array<OneD, Array<OneD, NekDouble> > Bwd (nConvectiveFields);
827  Array<OneD, Array<OneD, NekDouble> > numflux(nConvectiveFields);
828 
829  // Set up storage for flux vector.
830  for (i = 0; i < nConvectiveFields; ++i)
831  {
832  fluxvector[i] =
833  Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
834 
835  for (j = 0; j < m_spaceDim; ++j)
836  {
837  fluxvector[i][j] = Array<OneD, NekDouble>(nSolutionPts);
838  }
839  }
840 
841  for (i = 0; i < nConvectiveFields; ++i)
842  {
843  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
844  Fwd[i] = Array<OneD, NekDouble>(nTracePts);
845  Bwd[i] = Array<OneD, NekDouble>(nTracePts);
846  numflux[i] = Array<OneD, NekDouble>(nTracePts);
847  fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
848  }
849 
850  // Computing the interface flux at each trace point
851  m_riemann->Solve(m_spaceDim, Fwd, Bwd, numflux);
852 
853  // Divergence of the flux (computing the RHS)
854  switch(nDimensions)
855  {
856  // 1D problems
857  case 1:
858  {
859  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
860  Array<OneD, NekDouble> divFC (nSolutionPts);
861 
862  // Get the advection flux vector (solver specific)
863  m_fluxVector(inarray, fluxvector);
864 
865  // Get the discontinuous flux divergence
866  for(i = 0; i < nConvectiveFields; ++i)
867  {
868  for (n = 0; n < nElements; n++)
869  {
870  phys_offset = fields[0]->GetPhys_Offset(n);
871 
872  fields[i]->GetExp(n)->PhysDeriv(
873  0, fluxvector[i][0] + phys_offset,
874  auxArray2 = DfluxvectorX1 + phys_offset);
875  }
876 
877  // Get the correction flux divergence
878  v_DivCFlux_1D(nConvectiveFields, fields,
879  fluxvector[i][0], numflux[i], divFC);
880 
881  // Back to the physical space using global operations
882  Vmath::Vdiv(nSolutionPts, &divFC[0], 1, &m_jac[0], 1,
883  &outarray[i][0], 1);
884 
885  // Adding the total divergence to outarray (RHS)
886  Vmath::Vadd(nSolutionPts, &outarray[i][0], 1,
887  &DfluxvectorX1[0], 1, &outarray[i][0], 1);
888 
889  // Primitive Dealiasing 1D
890  if (!(Basis[0]->Collocation()))
891  {
892  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
893  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
894  }
895  }
896  break;
897  }
898  // 2D problems
899  case 2:
900  {
901  Array<OneD, NekDouble> DfluxvectorX1(nSolutionPts);
902  Array<OneD, NekDouble> DfluxvectorX2(nSolutionPts);
903  Array<OneD, NekDouble> divFD(nSolutionPts);
904  Array<OneD, NekDouble> divFC(nSolutionPts);
905 
906  // Get the advection flux vector (solver specific)
907  m_fluxVector(inarray, fluxvector);
908 
909  for (i = 0; i < nConvectiveFields; ++i)
910  {
911  // Temporary vectors
912  Array<OneD, NekDouble> f_hat(nSolutionPts);
913  Array<OneD, NekDouble> g_hat(nSolutionPts);
914 
915  Vmath::Vvtvvtp(nSolutionPts,
916  &m_gmat[0][0], 1,
917  &fluxvector[i][0][0], 1,
918  &m_gmat[2][0], 1,
919  &fluxvector[i][1][0], 1,
920  &f_hat[0], 1);
921 
922  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &f_hat[0], 1,
923  &f_hat[0], 1);
924 
925  Vmath::Vvtvvtp(nSolutionPts,
926  &m_gmat[1][0], 1,
927  &fluxvector[i][0][0], 1,
928  &m_gmat[3][0], 1,
929  &fluxvector[i][1][0], 1,
930  &g_hat[0], 1);
931 
932  Vmath::Vmul(nSolutionPts, &m_jac[0], 1, &g_hat[0], 1,
933  &g_hat[0], 1);
934 
935  // Get the discontinuous flux derivatives
936  for (n = 0; n < nElements; n++)
937  {
938  phys_offset = fields[0]->GetPhys_Offset(n);
939  fields[0]->GetExp(n)->StdPhysDeriv(0,
940  auxArray1 = f_hat + phys_offset,
941  auxArray2 = DfluxvectorX1 + phys_offset);
942  fields[0]->GetExp(n)->StdPhysDeriv(1,
943  auxArray1 = g_hat + phys_offset,
944  auxArray2 = DfluxvectorX2 + phys_offset);
945  }
946 
947  // Divergence of the discontinuous flux
948  Vmath::Vadd(nSolutionPts, DfluxvectorX1, 1,
949  DfluxvectorX2, 1, divFD, 1);
950 
951  // Divergence of the correction flux
952  if (Basis[0]->GetPointsType() ==
953  LibUtilities::eGaussGaussLegendre &&
954  Basis[1]->GetPointsType() ==
955  LibUtilities::eGaussGaussLegendre)
956  {
957  v_DivCFlux_2D_Gauss(
958  nConvectiveFields, fields, f_hat, g_hat,
959  numflux[i], divFC);
960  }
961  else
962  {
963  v_DivCFlux_2D(
964  nConvectiveFields, fields,
965  fluxvector[i][0], fluxvector[i][1],
966  numflux[i], divFC);
967 
968  }
969 
970  // Divergence of the final flux
971  Vmath::Vadd(nSolutionPts, divFD, 1, divFC, 1,
972  outarray[i], 1);
973 
974  // Back to the physical space using a global operation
975  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
976  &m_jac[0], 1, &outarray[i][0], 1);
977 
978  // Primitive Dealiasing 2D
979  if (!(Basis[0]->Collocation()))
980  {
981  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
982  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
983  }
984  } // close nConvectiveFields loop
985  break;
986  }
987  // 3D problems
988  case 3:
989  {
990  ASSERTL0(false,"3D FRDG case not implemented yet");
991  break;
992  }
993  }
994  }
995 
996  /**
997  * @brief Compute the divergence of the corrective flux for 1D problems.
998  *
999  * @param nConvectiveFields Number of fields.
1000  * @param fields Pointer to fields.
1001  * @param fluxX1 X1-volumetric flux in physical space.
1002  * @param numericalFlux Interface flux in physical space.
1003  * @param divCFlux Divergence of the corrective flux.
1004  *
1005  */
1006  void AdvectionFR::v_DivCFlux_1D(
1007  const int nConvectiveFields,
1008  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1009  const Array<OneD, const NekDouble> &fluxX1,
1010  const Array<OneD, const NekDouble> &numericalFlux,
1011  Array<OneD, NekDouble> &divCFlux)
1012  {
1013  int n;
1014  int nLocalSolutionPts, phys_offset, t_offset;
1015 
1016  LibUtilities::BasisSharedPtr Basis;
1017  Basis = fields[0]->GetExp(0)->GetBasis(0);
1018 
1019  int nElements = fields[0]->GetExpSize();
1020  int nSolutionPts = fields[0]->GetTotPoints();
1021 
1022 
1023  vector<bool> negatedFluxNormal = (boost::static_pointer_cast<MultiRegions::DisContField1D>(fields[0]))->GetNegatedFluxNormal();
1024 
1025  // Arrays to store the derivatives of the correction flux
1026  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1027  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1028 
1029  // Arrays to store the intercell numerical flux jumps
1030  Array<OneD, NekDouble> JumpL(nElements);
1031  Array<OneD, NekDouble> JumpR(nElements);
1032 
1033  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1034  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1035 
1036  for (n = 0; n < nElements; ++n)
1037  {
1038  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1039  phys_offset = fields[0]->GetPhys_Offset(n);
1040 
1041  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1042  NekDouble tmpFluxVertex = 0;
1043  Vmath::Vcopy(nLocalSolutionPts,
1044  &fluxX1[phys_offset], 1,
1045  &tmparrayX1[0], 1);
1046 
1047  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1048  tmpFluxVertex);
1049 
1050  t_offset = fields[0]->GetTrace()
1051  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1052 
1053  if(negatedFluxNormal[2*n])
1054  {
1055  JumpL[n] = numericalFlux[t_offset] - tmpFluxVertex;
1056  }
1057  else
1058  {
1059  JumpL[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1060  }
1061 
1062  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1063  tmpFluxVertex);
1064 
1065  t_offset = fields[0]->GetTrace()
1066  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1067 
1068  if(negatedFluxNormal[2*n+1])
1069  {
1070  JumpR[n] = -numericalFlux[t_offset] - tmpFluxVertex;
1071  }
1072  else
1073  {
1074  JumpR[n] = numericalFlux[t_offset] - tmpFluxVertex;
1075  }
1076  }
1077 
1078  for (n = 0; n < nElements; ++n)
1079  {
1080  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1081  phys_offset = fields[0]->GetPhys_Offset(n);
1082 
1083  // Left jump multiplied by left derivative of C function
1084  Vmath::Smul(nLocalSolutionPts, JumpL[n], &m_dGL_xi1[n][0], 1,
1085  &DCL[0], 1);
1086 
1087  // Right jump multiplied by right derivative of C function
1088  Vmath::Smul(nLocalSolutionPts, JumpR[n], &m_dGR_xi1[n][0], 1,
1089  &DCR[0], 1);
1090 
1091  // Assembling divergence of the correction flux
1092  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1093  &divCFlux[phys_offset], 1);
1094  }
1095  }
1096 
1097  /**
1098  * @brief Compute the divergence of the corrective flux for 2D problems.
1099  *
1100  * @param nConvectiveFields Number of fields.
1101  * @param fields Pointer to fields.
1102  * @param fluxX1 X1-volumetric flux in physical space.
1103  * @param fluxX2 X2-volumetric flux in physical space.
1104  * @param numericalFlux Interface flux in physical space.
1105  * @param divCFlux Divergence of the corrective flux.
1106  *
1107  * \todo: Switch on shapes eventually here.
1108  */
1109  void AdvectionFR::v_DivCFlux_2D(
1110  const int nConvectiveFields,
1111  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1112  const Array<OneD, const NekDouble> &fluxX1,
1113  const Array<OneD, const NekDouble> &fluxX2,
1114  const Array<OneD, const NekDouble> &numericalFlux,
1115  Array<OneD, NekDouble> &divCFlux)
1116  {
1117  int n, e, i, j, cnt;
1118 
1119  int nElements = fields[0]->GetExpSize();
1120 
1121  int nLocalSolutionPts, nEdgePts;
1122  int trace_offset, phys_offset;
1123  int nquad0, nquad1;
1124 
1125  Array<OneD, NekDouble> auxArray1;
1126  Array<OneD, LibUtilities::BasisSharedPtr> base;
1127 
1128  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1129  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1130 
1131  // Loop on the elements
1132  for(n = 0; n < nElements; ++n)
1133  {
1134  // Offset of the element on the global vector
1135  phys_offset = fields[0]->GetPhys_Offset(n);
1136  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1137 
1138  base = fields[0]->GetExp(n)->GetBase();
1139  nquad0 = base[0]->GetNumPoints();
1140  nquad1 = base[1]->GetNumPoints();
1141 
1142  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1143  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1144  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1145  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1146 
1147  // Loop on the edges
1148  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1149  {
1150  // Number of edge points of edge e
1151  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1152 
1153  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1154  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1155  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1156  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1157  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1158 
1159  // Offset of the trace space correspondent to edge e
1160  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1161  elmtToTrace[n][e]->GetElmtId());
1162 
1163  // Get the normals of edge e
1164  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1165  fields[0]->GetExp(n)->GetEdgeNormal(e);
1166 
1167  // Extract the edge values of flux-x on edge e and order
1168  // them accordingly to the order of the trace space
1169  fields[0]->GetExp(n)->GetEdgePhysVals(
1170  e, elmtToTrace[n][e],
1171  fluxX1 + phys_offset,
1172  auxArray1 = tmparrayX1);
1173 
1174  // Extract the edge values of flux-y on edge e and order
1175  // them accordingly to the order of the trace space
1176  fields[0]->GetExp(n)->GetEdgePhysVals(
1177  e, elmtToTrace[n][e],
1178  fluxX2 + phys_offset,
1179  auxArray1 = tmparrayX2);
1180 
1181  // Multiply the edge components of the flux by the normal
1182  Vmath::Vvtvvtp(nEdgePts, &tmparrayX1[0], 1,
1183  &m_traceNormals[0][trace_offset], 1,
1184  &tmparrayX2[0], 1,
1185  &m_traceNormals[1][trace_offset], 1,
1186  &fluxN[0], 1);
1187 
1188  // Subtract to the Riemann flux the discontinuous flux
1189  Vmath::Vsub(nEdgePts, &numericalFlux[trace_offset], 1,
1190  &fluxN[0], 1, &fluxJumps[0], 1);
1191 
1192  // Check the ordering of the jump vectors
1193  if (fields[0]->GetExp(n)->GetEorient(e) ==
1194  StdRegions::eBackwards)
1195  {
1196  Vmath::Reverse(nEdgePts, &fluxJumps[0], 1,
1197  &fluxJumps[0], 1);
1198  }
1199 
1200  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1201  -1.0 : 1.0;
1202 
1203  for (i = 0; i < nEdgePts; ++i)
1204  {
1205  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1206  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1207  {
1208  fluxJumps[i] = -fluxJumps[i];
1209  }
1210  }
1211 
1212  // Multiply jumps by derivatives of the correction functions
1213  switch (e)
1214  {
1215  case 0:
1216  for (i = 0; i < nquad0; ++i)
1217  {
1218  // Multiply fluxJumps by Q factors
1219  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1220 
1221  for (j = 0; j < nquad1; ++j)
1222  {
1223  cnt = i + j*nquad0;
1224  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1225  }
1226  }
1227  break;
1228  case 1:
1229  for (i = 0; i < nquad1; ++i)
1230  {
1231  // Multiply fluxJumps by Q factors
1232  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1233 
1234  for (j = 0; j < nquad0; ++j)
1235  {
1236  cnt = (nquad0)*i + j;
1237  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1238  }
1239  }
1240  break;
1241  case 2:
1242  for (i = 0; i < nquad0; ++i)
1243  {
1244  // Multiply fluxJumps by Q factors
1245  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1246 
1247  for (j = 0; j < nquad1; ++j)
1248  {
1249  cnt = (nquad0 - 1) + j*nquad0 - i;
1250  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1251  }
1252  }
1253  break;
1254  case 3:
1255  for (i = 0; i < nquad1; ++i)
1256  {
1257  // Multiply fluxJumps by Q factors
1258  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1259  for (j = 0; j < nquad0; ++j)
1260  {
1261  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1262  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1263  }
1264  }
1265  break;
1266  default:
1267  ASSERTL0(false,"edge value (< 3) is out of range");
1268  break;
1269  }
1270  }
1271 
1272  // Sum each edge contribution
1273  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1274  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1275 
1276  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1277  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1278 
1279  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1280  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1281  }
1282  }
1283 
1284  /**
1285  * @brief Compute the divergence of the corrective flux for 2D problems
1286  * where POINTSTYPE="GaussGaussLegendre"
1287  *
1288  * @param nConvectiveFields Number of fields.
1289  * @param fields Pointer to fields.
1290  * @param fluxX1 X1-volumetric flux in physical space.
1291  * @param fluxX2 X2-volumetric flux in physical space.
1292  * @param numericalFlux Interface flux in physical space.
1293  * @param divCFlux Divergence of the corrective flux.
1294  *
1295  * \todo: Switch on shapes eventually here.
1296  */
1297 
1298  void AdvectionFR::v_DivCFlux_2D_Gauss(
1299  const int nConvectiveFields,
1300  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1301  const Array<OneD, const NekDouble> &fluxX1,
1302  const Array<OneD, const NekDouble> &fluxX2,
1303  const Array<OneD, const NekDouble> &numericalFlux,
1304  Array<OneD, NekDouble> &divCFlux)
1305  {
1306  int n, e, i, j, cnt;
1307 
1308  int nElements = fields[0]->GetExpSize();
1309  int nLocalSolutionPts;
1310  int nEdgePts;
1311  int trace_offset;
1312  int phys_offset;
1313  int nquad0;
1314  int nquad1;
1315 
1316  NekDouble fac;
1317 
1318  Array<OneD, NekDouble> auxArray1, auxArray2;
1319  Array<OneD, LibUtilities::BasisSharedPtr> base;
1320 
1321  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1322  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1323 
1324  // Loop on the elements
1325  for(n = 0; n < nElements; ++n)
1326  {
1327  // Offset of the element on the global vector
1328  phys_offset = fields[0]->GetPhys_Offset(n);
1329  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1330 
1331  base = fields[0]->GetExp(n)->GetBase();
1332  nquad0 = base[0]->GetNumPoints();
1333  nquad1 = base[1]->GetNumPoints();
1334 
1335  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1336  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1337  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1338  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1339 
1340  // Loop on the edges
1341  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1342  {
1343  // Number of edge points of edge e
1344  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1345 
1346  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1347  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1348  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
1349  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
1350  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
1351 
1352  // Offset of the trace space correspondent to edge e
1353  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1354  elmtToTrace[n][e]->GetElmtId());
1355 
1356  // Get the normals of edge e
1357  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1358  fields[0]->GetExp(n)->GetEdgeNormal(e);
1359 
1360  // Extract the trasformed normal flux at each edge
1361  switch (e)
1362  {
1363  case 0:
1364  // Extract the edge values of transformed flux-y on
1365  // edge e and order them accordingly to the order of
1366  // the trace space
1367  fields[0]->GetExp(n)->GetEdgePhysVals(
1368  e, elmtToTrace[n][e],
1369  fluxX2 + phys_offset,
1370  auxArray1 = fluxN_D);
1371 
1372  Vmath::Neg (nEdgePts, fluxN_D, 1);
1373 
1374  // Extract the physical Rieamann flux at each edge
1375  Vmath::Vcopy(nEdgePts,
1376  &numericalFlux[trace_offset], 1,
1377  &fluxN[0], 1);
1378 
1379  // Check the ordering of vectors
1380  if (fields[0]->GetExp(n)->GetEorient(e) ==
1381  StdRegions::eBackwards)
1382  {
1383  Vmath::Reverse(nEdgePts,
1384  auxArray1 = fluxN, 1,
1385  auxArray2 = fluxN, 1);
1386 
1387  Vmath::Reverse(nEdgePts,
1388  auxArray1 = fluxN_D, 1,
1389  auxArray2 = fluxN_D, 1);
1390  }
1391 
1392  // Transform Riemann Fluxes in the standard element
1393  for (i = 0; i < nquad0; ++i)
1394  {
1395  // Multiply Riemann Flux by Q factors
1396  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
1397  }
1398 
1399  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1400  -1.0 : 1.0;
1401 
1402  for (i = 0; i < nEdgePts; ++i)
1403  {
1404  if (m_traceNormals[0][trace_offset+i]
1405  != fac*normals[0][i] ||
1406  m_traceNormals[1][trace_offset+i]
1407  != fac*normals[1][i])
1408  {
1409  fluxN_R[i] = -fluxN_R[i];
1410  }
1411  }
1412 
1413  // Subtract to the Riemann flux the discontinuous
1414  // flux
1415  Vmath::Vsub(nEdgePts,
1416  &fluxN_R[0], 1,
1417  &fluxN_D[0], 1, &fluxJumps[0], 1);
1418 
1419  // Multiplicate the flux jumps for the correction
1420  // function
1421  for (i = 0; i < nquad0; ++i)
1422  {
1423  for (j = 0; j < nquad1; ++j)
1424  {
1425  cnt = i + j*nquad0;
1426  divCFluxE0[cnt] =
1427  -fluxJumps[i] * m_dGL_xi2[n][j];
1428  }
1429  }
1430  break;
1431  case 1:
1432  // Extract the edge values of transformed flux-x on
1433  // edge e and order them accordingly to the order of
1434  // the trace space
1435  fields[0]->GetExp(n)->GetEdgePhysVals(
1436  e, elmtToTrace[n][e],
1437  fluxX1 + phys_offset,
1438  auxArray1 = fluxN_D);
1439 
1440  // Extract the physical Rieamann flux at each edge
1441  Vmath::Vcopy(nEdgePts,
1442  &numericalFlux[trace_offset], 1,
1443  &fluxN[0], 1);
1444 
1445  // Check the ordering of vectors
1446  if (fields[0]->GetExp(n)->GetEorient(e) ==
1447  StdRegions::eBackwards)
1448  {
1449  Vmath::Reverse(nEdgePts,
1450  auxArray1 = fluxN, 1,
1451  auxArray2 = fluxN, 1);
1452 
1453  Vmath::Reverse(nEdgePts,
1454  auxArray1 = fluxN_D, 1,
1455  auxArray2 = fluxN_D, 1);
1456  }
1457 
1458  // Transform Riemann Fluxes in the standard element
1459  for (i = 0; i < nquad1; ++i)
1460  {
1461  // Multiply Riemann Flux by Q factors
1462  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
1463  }
1464 
1465  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1466  -1.0 : 1.0;
1467 
1468  for (i = 0; i < nEdgePts; ++i)
1469  {
1470  if (m_traceNormals[0][trace_offset+i]
1471  != fac*normals[0][i] ||
1472  m_traceNormals[1][trace_offset+i]
1473  != fac*normals[1][i])
1474  {
1475  fluxN_R[i] = -fluxN_R[i];
1476  }
1477  }
1478 
1479  // Subtract to the Riemann flux the discontinuous
1480  // flux
1481  Vmath::Vsub(nEdgePts,
1482  &fluxN_R[0], 1,
1483  &fluxN_D[0], 1, &fluxJumps[0], 1);
1484 
1485  // Multiplicate the flux jumps for the correction
1486  // function
1487  for (i = 0; i < nquad1; ++i)
1488  {
1489  for (j = 0; j < nquad0; ++j)
1490  {
1491  cnt = (nquad0)*i + j;
1492  divCFluxE1[cnt] =
1493  fluxJumps[i] * m_dGR_xi1[n][j];
1494  }
1495  }
1496  break;
1497  case 2:
1498 
1499  // Extract the edge values of transformed flux-y on
1500  // edge e and order them accordingly to the order of
1501  // the trace space
1502 
1503  fields[0]->GetExp(n)->GetEdgePhysVals(
1504  e, elmtToTrace[n][e],
1505  fluxX2 + phys_offset,
1506  auxArray1 = fluxN_D);
1507 
1508  // Extract the physical Rieamann flux at each edge
1509  Vmath::Vcopy(nEdgePts,
1510  &numericalFlux[trace_offset], 1,
1511  &fluxN[0], 1);
1512 
1513  // Check the ordering of vectors
1514  if (fields[0]->GetExp(n)->GetEorient(e) ==
1515  StdRegions::eBackwards)
1516  {
1517  Vmath::Reverse(nEdgePts,
1518  auxArray1 = fluxN, 1,
1519  auxArray2 = fluxN, 1);
1520 
1521  Vmath::Reverse(nEdgePts,
1522  auxArray1 = fluxN_D, 1,
1523  auxArray2 = fluxN_D, 1);
1524  }
1525 
1526  // Transform Riemann Fluxes in the standard element
1527  for (i = 0; i < nquad0; ++i)
1528  {
1529  // Multiply Riemann Flux by Q factors
1530  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
1531  }
1532 
1533  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1534  -1.0 : 1.0;
1535 
1536  for (i = 0; i < nEdgePts; ++i)
1537  {
1538  if (m_traceNormals[0][trace_offset+i]
1539  != fac*normals[0][i] ||
1540  m_traceNormals[1][trace_offset+i]
1541  != fac*normals[1][i])
1542  {
1543  fluxN_R[i] = -fluxN_R[i];
1544  }
1545  }
1546 
1547  // Subtract to the Riemann flux the discontinuous
1548  // flux
1549 
1550  Vmath::Vsub(nEdgePts,
1551  &fluxN_R[0], 1,
1552  &fluxN_D[0], 1, &fluxJumps[0], 1);
1553 
1554  // Multiplicate the flux jumps for the correction
1555  // function
1556  for (i = 0; i < nquad0; ++i)
1557  {
1558  for (j = 0; j < nquad1; ++j)
1559  {
1560  cnt = (nquad0 - 1) + j*nquad0 - i;
1561  divCFluxE2[cnt] =
1562  fluxJumps[i] * m_dGR_xi2[n][j];
1563  }
1564  }
1565  break;
1566  case 3:
1567  // Extract the edge values of transformed flux-x on
1568  // edge e and order them accordingly to the order of
1569  // the trace space
1570 
1571  fields[0]->GetExp(n)->GetEdgePhysVals(
1572  e, elmtToTrace[n][e],
1573  fluxX1 + phys_offset,
1574  auxArray1 = fluxN_D);
1575  Vmath::Neg (nEdgePts, fluxN_D, 1);
1576 
1577  // Extract the physical Rieamann flux at each edge
1578  Vmath::Vcopy(nEdgePts,
1579  &numericalFlux[trace_offset], 1,
1580  &fluxN[0], 1);
1581 
1582  // Check the ordering of vectors
1583  if (fields[0]->GetExp(n)->GetEorient(e) ==
1584  StdRegions::eBackwards)
1585  {
1586  Vmath::Reverse(nEdgePts,
1587  auxArray1 = fluxN, 1,
1588  auxArray2 = fluxN, 1);
1589 
1590  Vmath::Reverse(nEdgePts,
1591  auxArray1 = fluxN_D, 1,
1592  auxArray2 = fluxN_D, 1);
1593  }
1594 
1595  // Transform Riemann Fluxes in the standard element
1596  for (i = 0; i < nquad1; ++i)
1597  {
1598  // Multiply Riemann Flux by Q factors
1599  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
1600  }
1601 
1602  fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1603  -1.0 : 1.0;
1604 
1605  for (i = 0; i < nEdgePts; ++i)
1606  {
1607  if (m_traceNormals[0][trace_offset+i]
1608  != fac*normals[0][i] ||
1609  m_traceNormals[1][trace_offset+i]
1610  != fac*normals[1][i])
1611  {
1612  fluxN_R[i] = -fluxN_R[i];
1613  }
1614  }
1615 
1616  // Subtract to the Riemann flux the discontinuous
1617  // flux
1618 
1619  Vmath::Vsub(nEdgePts,
1620  &fluxN_R[0], 1,
1621  &fluxN_D[0], 1, &fluxJumps[0], 1);
1622 
1623  // Multiplicate the flux jumps for the correction
1624  // function
1625  for (i = 0; i < nquad1; ++i)
1626  {
1627  for (j = 0; j < nquad0; ++j)
1628  {
1629  cnt = (nquad0*nquad1 - nquad0) + j
1630  - i*nquad0;
1631  divCFluxE3[cnt] =
1632  -fluxJumps[i] * m_dGL_xi1[n][j];
1633  }
1634  }
1635  break;
1636  default:
1637  ASSERTL0(false,"edge value (< 3) is out of range");
1638  break;
1639  }
1640  }
1641 
1642 
1643  // Sum each edge contribution
1644  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1645  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1646 
1647  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1648  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1649 
1650  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1651  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1652  }
1653  }
1654 
1655 
1656  /**
1657  * @brief Compute the divergence of the corrective flux for 3D problems.
1658  *
1659  * @param nConvectiveFields Number of fields.
1660  * @param fields Pointer to fields.
1661  * @param fluxX1 X1-volumetric flux in physical space.
1662  * @param fluxX2 X2-volumetric flux in physical space.
1663  * @param fluxX3 X3-volumetric flux in physical space.
1664  * @param numericalFlux Interface flux in physical space.
1665  * @param divCFlux Divergence of the corrective flux.
1666  *
1667  * \todo: To be implemented. Switch on shapes eventually here.
1668  */
1669  void AdvectionFR::v_DivCFlux_3D(
1670  const int nConvectiveFields,
1671  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1672  const Array<OneD, const NekDouble> &fluxX1,
1673  const Array<OneD, const NekDouble> &fluxX2,
1674  const Array<OneD, const NekDouble> &fluxX3,
1675  const Array<OneD, const NekDouble> &numericalFlux,
1676  Array<OneD, NekDouble> &divCFlux)
1677  {
1678 
1679  }
1680 
1681  }
1682 }
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2