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