Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DiffusionLFR.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: DiffusionLFR.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: LFR diffusion class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LocalRegions/Expansion1D.h>
37 #include <LocalRegions/Expansion2D.h>
38 #include <SolverUtils/Diffusion/DiffusionLFR.h>
39 #include <LibUtilities/Polylib/Polylib.h>
40 #include <MultiRegions/DisContField1D.h>
41 #include <boost/math/special_functions/gamma.hpp>
42 #include <iostream>
43 #include <iomanip>
44 
45 namespace Nektar
46 {
47  namespace SolverUtils
48  {
49  std::string DiffusionLFR::type[] = {
50  GetDiffusionFactory().RegisterCreatorFunction(
51  "LFRDG", DiffusionLFR::create),
52  GetDiffusionFactory().RegisterCreatorFunction(
53  "LFRSD", DiffusionLFR::create),
54  GetDiffusionFactory().RegisterCreatorFunction(
55  "LFRHU", DiffusionLFR::create),
56  GetDiffusionFactory().RegisterCreatorFunction(
57  "LFRcmin", DiffusionLFR::create),
58  GetDiffusionFactory().RegisterCreatorFunction(
59  "LFRcinf", DiffusionLFR::create)};
60 
61  /**
62  * @brief DiffusionLFR uses the Flux Reconstruction (FR) approach to
63  * compute the diffusion term. The implementation is only for segments,
64  * quadrilaterals and hexahedra at the moment.
65  *
66  * \todo Extension to triangles, tetrahedra and other shapes.
67  * (Long term objective)
68  */
69  DiffusionLFR::DiffusionLFR(std::string diffType):m_diffType(diffType)
70  {
71  }
72 
73  /**
74  * @brief Initiliase DiffusionLFR objects and store them before starting
75  * the time-stepping.
76  *
77  * This routine calls the virtual functions #v_SetupMetrics and
78  * #v_SetupCFunctions to initialise the objects needed by DiffusionLFR.
79  *
80  * @param pSession Pointer to session reader.
81  * @param pFields Pointer to fields.
82  */
83  void DiffusionLFR::v_InitObject(
84  LibUtilities::SessionReaderSharedPtr pSession,
85  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
86  {
87  m_session = pSession;
88  v_SetupMetrics(pSession, pFields);
89  v_SetupCFunctions(pSession, pFields);
90 
91  // Initialising arrays
92  int i, j;
93  int nConvectiveFields = pFields.num_elements();
94  int nDim = pFields[0]->GetCoordim(0);
95  int nSolutionPts = pFields[0]->GetTotPoints();
96  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
97 
98  m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
99  nConvectiveFields);
100  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
101  nConvectiveFields);
102  m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
103  nConvectiveFields);
104  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
105  nConvectiveFields);
106  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
107  nConvectiveFields);
108  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
109  nConvectiveFields);
110  m_IF2 = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
111  m_DFC2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
112  nConvectiveFields);
113  m_divFD = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
114  m_divFC = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
115 
116  m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
117  nConvectiveFields);
118  m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
119  nConvectiveFields);
120 
121 
122 
123  for (i = 0; i < nConvectiveFields; ++i)
124  {
125  m_IF1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
126  m_DU1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
127  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
128  m_BD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
129  m_D1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
130  m_DD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
131  m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
132  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
133  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
134  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
135 
136  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
137  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
138 
139  for (j = 0; j < nDim; ++j)
140  {
141  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
142  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
143  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
144  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
147  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
148 
149  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
150  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
151  }
152  }
153 
154  }
155 
156  /**
157  * @brief Setup the metric terms to compute the contravariant
158  * fluxes. (i.e. this special metric terms transform the fluxes
159  * at the interfaces of each element from the physical space to
160  * the standard space).
161  *
162  * This routine calls the function #GetEdgeQFactors to compute and
163  * store the metric factors following an anticlockwise conventions
164  * along the edges/faces of each element. Note: for 1D problem
165  * the transformation is not needed.
166  *
167  * @param pSession Pointer to session reader.
168  * @param pFields Pointer to fields.
169  *
170  * \todo Add the metric terms for 3D Hexahedra.
171  */
172  void DiffusionLFR::v_SetupMetrics(
173  LibUtilities::SessionReaderSharedPtr pSession,
174  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
175  {
176  int i, n;
177  int nquad0, nquad1;
178  int phys_offset;
179  int nLocalSolutionPts;
180  int nElements = pFields[0]->GetExpSize();
181  int nDimensions = pFields[0]->GetCoordim(0);
182  int nSolutionPts = pFields[0]->GetTotPoints();
183  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
184 
185  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
186  for (i = 0; i < nDimensions; ++i)
187  {
188  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
189  }
190  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
191 
192  Array<TwoD, const NekDouble> gmat;
193  Array<OneD, const NekDouble> jac;
194 
195  m_jac = Array<OneD, NekDouble>(nSolutionPts);
196 
197  Array<OneD, NekDouble> auxArray1;
198  Array<OneD, LibUtilities::BasisSharedPtr> base;
199  LibUtilities::PointsKeyVector ptsKeys;
200 
201  switch (nDimensions)
202  {
203  case 1:
204  {
205  for (n = 0; n < nElements; ++n)
206  {
207  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
208  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
209  phys_offset = pFields[0]->GetPhys_Offset(n);
210  jac = pFields[0]->GetExp(n)
211  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
212  ->GetMetricInfo()->GetJac(ptsKeys);
213  for (i = 0; i < nLocalSolutionPts; ++i)
214  {
215  m_jac[i+phys_offset] = jac[0];
216  }
217  }
218  break;
219  }
220  case 2:
221  {
222  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
223  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
224  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
225  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
226  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
227 
228  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
229  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
230  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
231  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
232 
233  for (n = 0; n < nElements; ++n)
234  {
235  base = pFields[0]->GetExp(n)->GetBase();
236  nquad0 = base[0]->GetNumPoints();
237  nquad1 = base[1]->GetNumPoints();
238 
239  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
240  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
241  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
242  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
243 
244  // Extract the Q factors at each edge point
245  pFields[0]->GetExp(n)->GetEdgeQFactors(
246  0, auxArray1 = m_Q2D_e0[n]);
247  pFields[0]->GetExp(n)->GetEdgeQFactors(
248  1, auxArray1 = m_Q2D_e1[n]);
249  pFields[0]->GetExp(n)->GetEdgeQFactors(
250  2, auxArray1 = m_Q2D_e2[n]);
251  pFields[0]->GetExp(n)->GetEdgeQFactors(
252  3, auxArray1 = m_Q2D_e3[n]);
253 
254  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
255  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
256  phys_offset = pFields[0]->GetPhys_Offset(n);
257 
258  jac = pFields[0]->GetExp(n)
259  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
260  ->GetMetricInfo()->GetJac(ptsKeys);
261  gmat = pFields[0]->GetExp(n)
262  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
263  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
264 
265  if (pFields[0]->GetExp(n)
266  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
267  ->GetMetricInfo()->GetGtype()
268  == SpatialDomains::eDeformed)
269  {
270  for (i = 0; i < nLocalSolutionPts; ++i)
271  {
272  m_jac[i+phys_offset] = jac[i];
273  m_gmat[0][i+phys_offset] = gmat[0][i];
274  m_gmat[1][i+phys_offset] = gmat[1][i];
275  m_gmat[2][i+phys_offset] = gmat[2][i];
276  m_gmat[3][i+phys_offset] = gmat[3][i];
277  }
278  }
279  else
280  {
281  for (i = 0; i < nLocalSolutionPts; ++i)
282  {
283  m_jac[i+phys_offset] = jac[0];
284  m_gmat[0][i+phys_offset] = gmat[0][0];
285  m_gmat[1][i+phys_offset] = gmat[1][0];
286  m_gmat[2][i+phys_offset] = gmat[2][0];
287  m_gmat[3][i+phys_offset] = gmat[3][0];
288  }
289  }
290  }
291  break;
292  }
293  case 3:
294  {
295  ASSERTL0(false,"3DFR Metric terms not implemented yet");
296  break;
297  }
298  default:
299  {
300  ASSERTL0(false, "Expansion dimension not recognised");
301  break;
302  }
303  }
304  }
305 
306  /**
307  * @brief Setup the derivatives of the correction functions. For more
308  * details see J Sci Comput (2011) 47: 50–72
309  *
310  * This routine calls 3 different bases:
311  * #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme,
312  * #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme,
313  * #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme.
314  * The values of the derivatives of the correction function are then
315  * stored into global variables and reused into the virtual functions
316  * #v_DivCFlux_1D, #v_DivCFlux_2D, #v_DivCFlux_3D to compute the
317  * the divergence of the correction flux for 1D, 2D or 3D problems.
318  *
319  * @param pSession Pointer to session reader.
320  * @param pFields Pointer to fields.
321  */
322  void DiffusionLFR::v_SetupCFunctions(
323  LibUtilities::SessionReaderSharedPtr pSession,
324  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
325  {
326  int i, n;
327  NekDouble c0, c1, c2;
328  int nquad0, nquad1, nquad2;
329  int nmodes0, nmodes1, nmodes2;
330  Array<OneD, LibUtilities::BasisSharedPtr> base;
331 
332  int nElements = pFields[0]->GetExpSize();
333  int nDim = pFields[0]->GetCoordim(0);
334 
335  switch (nDim)
336  {
337  case 1:
338  {
339  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
340  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
341 
342  for (n = 0; n < nElements; ++n)
343  {
344  base = pFields[0]->GetExp(n)->GetBase();
345  nquad0 = base[0]->GetNumPoints();
346  nmodes0 = base[0]->GetNumModes();
347  Array<OneD, const NekDouble> z0;
348  Array<OneD, const NekDouble> w0;
349 
350  base[0]->GetZW(z0, w0);
351 
352  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
353  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
354 
355  // Auxiliary vectors to build up the auxiliary
356  // derivatives of the Legendre polynomials
357  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
358  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
359  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
360 
361  // Degree of the correction functions
362  int p0 = nmodes0 - 1;
363 
364  // Function sign to compute left correction function
365  NekDouble sign0 = pow(-1.0, p0);
366 
367  // Factorial factor to build the scheme
368  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
369  / (pow(2.0, p0)
370  * boost::math::tgamma(p0 + 1)
371  * boost::math::tgamma(p0 + 1));
372 
373  // Scalar parameter which recovers the FR schemes
374  if (m_diffType == "LFRDG")
375  {
376  c0 = 0.0;
377  }
378  else if (m_diffType == "LFRSD")
379  {
380  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
381  * (ap0 * boost::math::tgamma(p0 + 1))
382  * (ap0 * boost::math::tgamma(p0 + 1)));
383  }
384  else if (m_diffType == "LFRHU")
385  {
386  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
387  * (ap0 * boost::math::tgamma(p0 + 1))
388  * (ap0 * boost::math::tgamma(p0 + 1)));
389  }
390  else if (m_diffType == "LFRcmin")
391  {
392  c0 = -2.0 / ((2.0 * p0 + 1.0)
393  * (ap0 * boost::math::tgamma(p0 + 1))
394  * (ap0 * boost::math::tgamma(p0 + 1)));
395  }
396  else if (m_diffType == "LFRinf")
397  {
398  c0 = 10000000000000000.0;
399  }
400 
401  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
402  * (ap0 * boost::math::tgamma(p0 + 1))
403  * (ap0 * boost::math::tgamma(p0 + 1));
404 
405  NekDouble overeta0 = 1.0 / (1.0 + etap0);
406 
407  // Derivative of the Legendre polynomials
408  // dLp = derivative of the Legendre polynomial of order p
409  // dLpp = derivative of the Legendre polynomial of order p+1
410  // dLpm = derivative of the Legendre polynomial of order p-1
411  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
412  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
413  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
414 
415  // Building the DG_c_Left
416  for(i = 0; i < nquad0; ++i)
417  {
418  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
419  m_dGL_xi1[n][i] += dLpp0[i];
420  m_dGL_xi1[n][i] *= overeta0;
421  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
422  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
423  }
424 
425  // Building the DG_c_Right
426  for(i = 0; i < nquad0; ++i)
427  {
428  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
429  m_dGR_xi1[n][i] += dLpp0[i];
430  m_dGR_xi1[n][i] *= overeta0;
431  m_dGR_xi1[n][i] += dLp0[i];
432  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
433  }
434  }
435  break;
436  }
437  case 2:
438  {
439  LibUtilities::BasisSharedPtr BasisFR_Left0;
440  LibUtilities::BasisSharedPtr BasisFR_Right0;
441  LibUtilities::BasisSharedPtr BasisFR_Left1;
442  LibUtilities::BasisSharedPtr BasisFR_Right1;
443 
444  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
445  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
446  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
447  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
448 
449  for (n = 0; n < nElements; ++n)
450  {
451  base = pFields[0]->GetExp(n)->GetBase();
452  nquad0 = base[0]->GetNumPoints();
453  nquad1 = base[1]->GetNumPoints();
454  nmodes0 = base[0]->GetNumModes();
455  nmodes1 = base[1]->GetNumModes();
456 
457  Array<OneD, const NekDouble> z0;
458  Array<OneD, const NekDouble> w0;
459  Array<OneD, const NekDouble> z1;
460  Array<OneD, const NekDouble> w1;
461 
462  base[0]->GetZW(z0, w0);
463  base[1]->GetZW(z1, w1);
464 
465  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
466  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
467  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
468  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
469 
470  // Auxiliary vectors to build up the auxiliary
471  // derivatives of the Legendre polynomials
472  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
473  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
474  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
475  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
476  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
477  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
478 
479  // Degree of the correction functions
480  int p0 = nmodes0 - 1;
481  int p1 = nmodes1 - 1;
482 
483  // Function sign to compute left correction function
484  NekDouble sign0 = pow(-1.0, p0);
485  NekDouble sign1 = pow(-1.0, p1);
486 
487  // Factorial factor to build the scheme
488  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
489  / (pow(2.0, p0)
490  * boost::math::tgamma(p0 + 1)
491  * boost::math::tgamma(p0 + 1));
492 
493  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
494  / (pow(2.0, p1)
495  * boost::math::tgamma(p1 + 1)
496  * boost::math::tgamma(p1 + 1));
497 
498  // Scalar parameter which recovers the FR schemes
499  if (m_diffType == "LFRDG")
500  {
501  c0 = 0.0;
502  c1 = 0.0;
503  }
504  else if (m_diffType == "LFRSD")
505  {
506  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
507  * (ap0 * boost::math::tgamma(p0 + 1))
508  * (ap0 * boost::math::tgamma(p0 + 1)));
509 
510  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
511  * (ap1 * boost::math::tgamma(p1 + 1))
512  * (ap1 * boost::math::tgamma(p1 + 1)));
513  }
514  else if (m_diffType == "LFRHU")
515  {
516  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
517  * (ap0 * boost::math::tgamma(p0 + 1))
518  * (ap0 * boost::math::tgamma(p0 + 1)));
519 
520  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
521  * (ap1 * boost::math::tgamma(p1 + 1))
522  * (ap1 * boost::math::tgamma(p1 + 1)));
523  }
524  else if (m_diffType == "LFRcmin")
525  {
526  c0 = -2.0 / ((2.0 * p0 + 1.0)
527  * (ap0 * boost::math::tgamma(p0 + 1))
528  * (ap0 * boost::math::tgamma(p0 + 1)));
529 
530  c1 = -2.0 / ((2.0 * p1 + 1.0)
531  * (ap1 * boost::math::tgamma(p1 + 1))
532  * (ap1 * boost::math::tgamma(p1 + 1)));
533  }
534  else if (m_diffType == "LFRinf")
535  {
536  c0 = 10000000000000000.0;
537  c1 = 10000000000000000.0;
538  }
539 
540  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
541  * (ap0 * boost::math::tgamma(p0 + 1))
542  * (ap0 * boost::math::tgamma(p0 + 1));
543 
544  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
545  * (ap1 * boost::math::tgamma(p1 + 1))
546  * (ap1 * boost::math::tgamma(p1 + 1));
547 
548  NekDouble overeta0 = 1.0 / (1.0 + etap0);
549  NekDouble overeta1 = 1.0 / (1.0 + etap1);
550 
551  // Derivative of the Legendre polynomials
552  // dLp = derivative of the Legendre polynomial of order p
553  // dLpp = derivative of the Legendre polynomial of order p+1
554  // dLpm = derivative of the Legendre polynomial of order p-1
555  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
556  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
557  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
558 
559  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
560  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
561  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
562 
563  // Building the DG_c_Left
564  for(i = 0; i < nquad0; ++i)
565  {
566  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
567  m_dGL_xi1[n][i] += dLpp0[i];
568  m_dGL_xi1[n][i] *= overeta0;
569  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
570  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
571  }
572 
573  // Building the DG_c_Left
574  for(i = 0; i < nquad1; ++i)
575  {
576  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
577  m_dGL_xi2[n][i] += dLpp1[i];
578  m_dGL_xi2[n][i] *= overeta1;
579  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
580  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
581  }
582 
583  // Building the DG_c_Right
584  for(i = 0; i < nquad0; ++i)
585  {
586  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
587  m_dGR_xi1[n][i] += dLpp0[i];
588  m_dGR_xi1[n][i] *= overeta0;
589  m_dGR_xi1[n][i] += dLp0[i];
590  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
591  }
592 
593  // Building the DG_c_Right
594  for(i = 0; i < nquad1; ++i)
595  {
596  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
597  m_dGR_xi2[n][i] += dLpp1[i];
598  m_dGR_xi2[n][i] *= overeta1;
599  m_dGR_xi2[n][i] += dLp1[i];
600  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
601  }
602  }
603  break;
604  }
605  case 3:
606  {
607  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
608  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
609  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
610  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
611  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
612  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
613 
614  for (n = 0; n < nElements; ++n)
615  {
616  base = pFields[0]->GetExp(n)->GetBase();
617  nquad0 = base[0]->GetNumPoints();
618  nquad1 = base[1]->GetNumPoints();
619  nquad2 = base[2]->GetNumPoints();
620  nmodes0 = base[0]->GetNumModes();
621  nmodes1 = base[1]->GetNumModes();
622  nmodes2 = base[2]->GetNumModes();
623 
624  Array<OneD, const NekDouble> z0;
625  Array<OneD, const NekDouble> w0;
626  Array<OneD, const NekDouble> z1;
627  Array<OneD, const NekDouble> w1;
628  Array<OneD, const NekDouble> z2;
629  Array<OneD, const NekDouble> w2;
630 
631  base[0]->GetZW(z0, w0);
632  base[1]->GetZW(z1, w1);
633  base[1]->GetZW(z2, w2);
634 
635  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
636  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
637  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
638  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
639  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
640  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
641 
642  // Auxiliary vectors to build up the auxiliary
643  // derivatives of the Legendre polynomials
644  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
645  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
646  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
647  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
648  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
649  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
650  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
651  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
652  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
653 
654  // Degree of the correction functions
655  int p0 = nmodes0 - 1;
656  int p1 = nmodes1 - 1;
657  int p2 = nmodes2 - 1;
658 
659  // Function sign to compute left correction function
660  NekDouble sign0 = pow(-1.0, p0);
661  NekDouble sign1 = pow(-1.0, p1);
662  NekDouble sign2 = pow(-1.0, p2);
663 
664  // Factorial factor to build the scheme
665  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
666  / (pow(2.0, p0)
667  * boost::math::tgamma(p0 + 1)
668  * boost::math::tgamma(p0 + 1));
669 
670  // Factorial factor to build the scheme
671  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
672  / (pow(2.0, p1)
673  * boost::math::tgamma(p1 + 1)
674  * boost::math::tgamma(p1 + 1));
675 
676  // Factorial factor to build the scheme
677  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
678  / (pow(2.0, p2)
679  * boost::math::tgamma(p2 + 1)
680  * boost::math::tgamma(p2 + 1));
681 
682  // Scalar parameter which recovers the FR schemes
683  if (m_diffType == "LFRDG")
684  {
685  c0 = 0.0;
686  c1 = 0.0;
687  c2 = 0.0;
688  }
689  else if (m_diffType == "LFRSD")
690  {
691  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
692  * (ap0 * boost::math::tgamma(p0 + 1))
693  * (ap0 * boost::math::tgamma(p0 + 1)));
694 
695  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
696  * (ap1 * boost::math::tgamma(p1 + 1))
697  * (ap1 * boost::math::tgamma(p1 + 1)));
698 
699  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
700  * (ap2 * boost::math::tgamma(p2 + 1))
701  * (ap2 * boost::math::tgamma(p2 + 1)));
702  }
703  else if (m_diffType == "LFRHU")
704  {
705  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
706  * (ap0 * boost::math::tgamma(p0 + 1))
707  * (ap0 * boost::math::tgamma(p0 + 1)));
708 
709  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
710  * (ap1 * boost::math::tgamma(p1 + 1))
711  * (ap1 * boost::math::tgamma(p1 + 1)));
712 
713  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
714  * (ap2 * boost::math::tgamma(p2 + 1))
715  * (ap2 * boost::math::tgamma(p2 + 1)));
716  }
717  else if (m_diffType == "LFRcmin")
718  {
719  c0 = -2.0 / ((2.0 * p0 + 1.0)
720  * (ap0 * boost::math::tgamma(p0 + 1))
721  * (ap0 * boost::math::tgamma(p0 + 1)));
722 
723  c1 = -2.0 / ((2.0 * p1 + 1.0)
724  * (ap1 * boost::math::tgamma(p1 + 1))
725  * (ap1 * boost::math::tgamma(p1 + 1)));
726 
727  c2 = -2.0 / ((2.0 * p2 + 1.0)
728  * (ap2 * boost::math::tgamma(p2 + 1))
729  * (ap2 * boost::math::tgamma(p2 + 1)));
730  }
731  else if (m_diffType == "LFRinf")
732  {
733  c0 = 10000000000000000.0;
734  c1 = 10000000000000000.0;
735  c2 = 10000000000000000.0;
736  }
737 
738  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
739  * (ap0 * boost::math::tgamma(p0 + 1))
740  * (ap0 * boost::math::tgamma(p0 + 1));
741 
742  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
743  * (ap1 * boost::math::tgamma(p1 + 1))
744  * (ap1 * boost::math::tgamma(p1 + 1));
745 
746  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
747  * (ap2 * boost::math::tgamma(p2 + 1))
748  * (ap2 * boost::math::tgamma(p2 + 1));
749 
750  NekDouble overeta0 = 1.0 / (1.0 + etap0);
751  NekDouble overeta1 = 1.0 / (1.0 + etap1);
752  NekDouble overeta2 = 1.0 / (1.0 + etap2);
753 
754  // Derivative of the Legendre polynomials
755  // dLp = derivative of the Legendre polynomial of order p
756  // dLpp = derivative of the Legendre polynomial of order p+1
757  // dLpm = derivative of the Legendre polynomial of order p-1
758  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
759  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
760  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
761 
762  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
763  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
764  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
765 
766  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
767  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
768  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
769 
770  // Building the DG_c_Left
771  for(i = 0; i < nquad0; ++i)
772  {
773  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
774  m_dGL_xi1[n][i] += dLpp0[i];
775  m_dGL_xi1[n][i] *= overeta0;
776  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
777  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
778  }
779 
780  // Building the DG_c_Left
781  for(i = 0; i < nquad1; ++i)
782  {
783  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
784  m_dGL_xi2[n][i] += dLpp1[i];
785  m_dGL_xi2[n][i] *= overeta1;
786  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
787  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
788  }
789 
790  // Building the DG_c_Left
791  for(i = 0; i < nquad2; ++i)
792  {
793  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
794  m_dGL_xi3[n][i] += dLpp2[i];
795  m_dGL_xi3[n][i] *= overeta2;
796  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
797  m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
798  }
799 
800  // Building the DG_c_Right
801  for(i = 0; i < nquad0; ++i)
802  {
803  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
804  m_dGR_xi1[n][i] += dLpp0[i];
805  m_dGR_xi1[n][i] *= overeta0;
806  m_dGR_xi1[n][i] += dLp0[i];
807  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
808  }
809 
810  // Building the DG_c_Right
811  for(i = 0; i < nquad1; ++i)
812  {
813  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
814  m_dGR_xi2[n][i] += dLpp1[i];
815  m_dGR_xi2[n][i] *= overeta1;
816  m_dGR_xi2[n][i] += dLp1[i];
817  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
818  }
819 
820  // Building the DG_c_Right
821  for(i = 0; i < nquad2; ++i)
822  {
823  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
824  m_dGR_xi3[n][i] += dLpp2[i];
825  m_dGR_xi3[n][i] *= overeta2;
826  m_dGR_xi3[n][i] += dLp2[i];
827  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
828  }
829  }
830  break;
831  }
832  default:
833  {
834  ASSERTL0(false,"Expansion dimension not recognised");
835  break;
836  }
837  }
838  }
839 
840  /**
841  * @brief Calculate FR Diffusion for the linear problems
842  * using an LDG interface flux.
843  *
844  */
845  void DiffusionLFR::v_Diffuse(
846  const int nConvectiveFields,
847  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
848  const Array<OneD, Array<OneD, NekDouble> > &inarray,
849  Array<OneD, Array<OneD, NekDouble> > &outarray)
850  {
851  int i, j, n;
852  int phys_offset;
853  //Array<TwoD, const NekDouble> gmat;
854  //Array<OneD, const NekDouble> jac;
855  Array<OneD, NekDouble> auxArray1, auxArray2;
856 
857  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
858  Basis = fields[0]->GetExp(0)->GetBase();
859 
860  int nElements = fields[0]->GetExpSize();
861  int nDim = fields[0]->GetCoordim(0);
862  int nSolutionPts = fields[0]->GetTotPoints();
863  int nCoeffs = fields[0]->GetNcoeffs();
864 
865  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
866  for (i = 0; i < nConvectiveFields; ++i)
867  {
868  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
869  }
870 
871  // Compute interface numerical fluxes for inarray in physical space
872  v_NumFluxforScalar(fields, inarray, m_IF1);
873 
874  switch(nDim)
875  {
876  // 1D problems
877  case 1:
878  {
879  for (i = 0; i < nConvectiveFields; ++i)
880  {
881  // Computing the physical first-order discountinuous
882  // derivative
883  for (n = 0; n < nElements; n++)
884  {
885  phys_offset = fields[0]->GetPhys_Offset(n);
886 
887  fields[i]->GetExp(n)->PhysDeriv(0,
888  auxArray1 = inarray[i] + phys_offset,
889  auxArray2 = m_DU1[i][0] + phys_offset);
890  }
891 
892  // Computing the standard first-order correction
893  // derivative
894  v_DerCFlux_1D(nConvectiveFields, fields, inarray[i],
895  m_IF1[i][0], m_DFC1[i][0]);
896 
897  // Back to the physical space using global operations
898  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
899  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
900  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
901  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
902 
903  // Computing total first order derivatives
904  Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1,
905  &m_DU1[i][0][0], 1, &m_D1[i][0][0], 1);
906 
907  Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1,
908  &m_tmp1[i][0][0], 1);
909  }
910 
911  // Computing interface numerical fluxes for m_D1
912  // in physical space
913  v_NumFluxforVector(fields, inarray, m_tmp1, m_IF2);
914 
915  for (i = 0; i < nConvectiveFields; ++i)
916  {
917  // Computing the physical second-order discountinuous
918  // derivative
919  for (n = 0; n < nElements; n++)
920  {
921  phys_offset = fields[0]->GetPhys_Offset(n);
922 
923  fields[i]->GetExp(n)->PhysDeriv(0,
924  auxArray1 = m_D1[i][0] + phys_offset,
925  auxArray2 = m_DD1[i][0] + phys_offset);
926  }
927 
928  // Computing the standard second-order correction
929  // derivative
930  v_DerCFlux_1D(nConvectiveFields, fields, m_D1[i][0],
931  m_IF2[i], m_DFC2[i][0]);
932 
933  // Back to the physical space using global operations
934  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
935  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
936  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
937  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
938 
939  // Adding the total divergence to outarray (RHS)
940  Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1,
941  &m_DD1[i][0][0], 1, &outarray[i][0], 1);
942 
943  // Primitive Dealiasing 1D
944  if (!(Basis[0]->Collocation()))
945  {
946  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
947  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
948  }
949  }
950  break;
951  }
952  // 2D problems
953  case 2:
954  {
955  for(i = 0; i < nConvectiveFields; ++i)
956  {
957  for (j = 0; j < nDim; ++j)
958  {
959  // Temporary vectors
960  Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
961  Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
962 
963  if (j == 0)
964  {
965  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
966  &m_gmat[0][0], 1, &u1_hat[0], 1);
967 
968  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
969  &m_jac[0], 1, &u1_hat[0], 1);
970 
971  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
972  &m_gmat[1][0], 1, &u2_hat[0], 1);
973 
974  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
975  &m_jac[0], 1, &u2_hat[0], 1);
976  }
977  else if (j == 1)
978  {
979  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
980  &m_gmat[2][0], 1, &u1_hat[0], 1);
981 
982  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
983  &m_jac[0], 1, &u1_hat[0], 1);
984 
985  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
986  &m_gmat[3][0], 1, &u2_hat[0], 1);
987 
988  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
989  &m_jac[0], 1, &u2_hat[0], 1);
990  }
991 
992  for (n = 0; n < nElements; n++)
993  {
994  phys_offset = fields[0]->GetPhys_Offset(n);
995 
996  fields[i]->GetExp(n)->StdPhysDeriv(0,
997  auxArray1 = u1_hat + phys_offset,
998  auxArray2 = m_tmp1[i][j] + phys_offset);
999 
1000  fields[i]->GetExp(n)->StdPhysDeriv(1,
1001  auxArray1 = u2_hat + phys_offset,
1002  auxArray2 = m_tmp2[i][j] + phys_offset);
1003  }
1004 
1005  Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1006  &m_tmp2[i][j][0], 1,
1007  &m_DU1[i][j][0], 1);
1008 
1009  // Divide by the metric jacobian
1010  Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1,
1011  &m_jac[0], 1, &m_DU1[i][j][0], 1);
1012 
1013  // Computing the standard first-order correction
1014  // derivatives
1015  v_DerCFlux_2D(nConvectiveFields, j, fields,
1016  inarray[i], m_IF1[i][j],
1017  m_DFC1[i][j]);
1018  }
1019 
1020  // Multiplying derivatives by B matrix to get auxiliary
1021  // variables
1022  for (j = 0; j < nSolutionPts; ++j)
1023  {
1024  // std(q1)
1025  m_BD1[i][0][j] =
1026  (m_gmat[0][j]*m_gmat[0][j] +
1027  m_gmat[2][j]*m_gmat[2][j]) *
1028  m_DFC1[i][0][j] +
1029  (m_gmat[1][j]*m_gmat[0][j] +
1030  m_gmat[3][j]*m_gmat[2][j]) *
1031  m_DFC1[i][1][j];
1032 
1033  // std(q2)
1034  m_BD1[i][1][j] =
1035  (m_gmat[1][j]*m_gmat[0][j] +
1036  m_gmat[3][j]*m_gmat[2][j]) *
1037  m_DFC1[i][0][j] +
1038  (m_gmat[1][j]*m_gmat[1][j] +
1039  m_gmat[3][j]*m_gmat[3][j]) *
1040  m_DFC1[i][1][j];
1041  }
1042 
1043  // Multiplying derivatives by A^(-1) to get back
1044  // into the physical space
1045  for (j = 0; j < nSolutionPts; j++)
1046  {
1047  // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1048  m_DFC1[i][0][j] =
1049  m_gmat[3][j] * m_BD1[i][0][j] -
1050  m_gmat[2][j] * m_BD1[i][1][j];
1051 
1052  // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1053  m_DFC1[i][1][j] =
1054  m_gmat[0][j] * m_BD1[i][1][j] -
1055  m_gmat[1][j] * m_BD1[i][0][j];
1056  }
1057 
1058  // Computing the physical first-order derivatives
1059  for (j = 0; j < nDim; ++j)
1060  {
1061  Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1062  &m_DFC1[i][j][0], 1,
1063  &m_D1[i][j][0], 1);
1064  }
1065  }
1066 
1067  // Computing interface numerical fluxes for m_D1
1068  // in physical space
1069  v_NumFluxforVector(fields, inarray, m_D1, m_IF2);
1070 
1071  // Computing the standard second-order discontinuous
1072  // derivatives
1073  for (i = 0; i < nConvectiveFields; ++i)
1074  {
1075  Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1076  Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1077 
1078  for (j = 0; j < nSolutionPts; j++)
1079  {
1080  f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1081  m_D1[i][1][j] * m_gmat[2][j])*m_jac[j];
1082 
1083  g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1084  m_D1[i][1][j] * m_gmat[3][j])*m_jac[j];
1085  }
1086 
1087  for (n = 0; n < nElements; n++)
1088  {
1089  phys_offset = fields[0]->GetPhys_Offset(n);
1090 
1091  fields[0]->GetExp(n)->StdPhysDeriv(0,
1092  auxArray1 = f_hat + phys_offset,
1093  auxArray2 = m_DD1[i][0] + phys_offset);
1094 
1095  fields[0]->GetExp(n)->StdPhysDeriv(1,
1096  auxArray1 = g_hat + phys_offset,
1097  auxArray2 = m_DD1[i][1] + phys_offset);
1098  }
1099 
1100  // Divergence of the standard discontinuous flux
1101  Vmath::Vadd(nSolutionPts,
1102  &m_DD1[i][0][0], 1,
1103  &m_DD1[i][1][0], 1,
1104  &m_divFD[i][0], 1);
1105 
1106  // Divergence of the standard correction flux
1107  if (Basis[0]->GetPointsType() ==
1108  LibUtilities::eGaussGaussLegendre &&
1109  Basis[1]->GetPointsType() ==
1110  LibUtilities::eGaussGaussLegendre)
1111  {
1112  v_DivCFlux_2D_Gauss(nConvectiveFields, fields,
1113  f_hat, g_hat,
1114  m_IF2[i], m_divFC[i]);
1115  }
1116  else
1117  {
1118  v_DivCFlux_2D(nConvectiveFields, fields,
1119  m_D1[i][0], m_D1[i][1],
1120  m_IF2[i], m_divFC[i]);
1121  }
1122 
1123  // Divergence of the standard final flux
1124  Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1,
1125  &m_divFC[i][0], 1, &outarray[i][0], 1);
1126 
1127  // Dividing by the jacobian to get back into
1128  // physical space
1129  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
1130  &m_jac[0], 1, &outarray[i][0], 1);
1131 
1132  // Primitive Dealiasing 2D
1133  if (!(Basis[0]->Collocation()))
1134  {
1135  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1136  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1137  }
1138  }//Close loop on nConvectiveFields
1139  break;
1140  }
1141  // 3D problems
1142  case 3:
1143  {
1144  ASSERTL0(false, "3D FRDG case not implemented yet");
1145  break;
1146  }
1147  }
1148  }
1149 
1150  /**
1151  * @brief Builds the numerical flux for the 1st order derivatives
1152  *
1153  */
1154  void DiffusionLFR::v_NumFluxforScalar(
1155  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1156  const Array<OneD, Array<OneD, NekDouble> > &ufield,
1157  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
1158  {
1159  int i, j;
1160  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1161  int nvariables = fields.num_elements();
1162  int nDim = fields[0]->GetCoordim(0);
1163 
1164  Array<OneD, NekDouble > Fwd (nTracePts);
1165  Array<OneD, NekDouble > Bwd (nTracePts);
1166  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1167  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1168 
1169  // Get the normal velocity Vn
1170  for (i = 0; i < nDim; ++i)
1171  {
1172  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1173  Vn, 1, Vn, 1);
1174  }
1175 
1176  // Get the sign of (v \cdot n), v = an arbitrary vector
1177  // Evaluate upwind flux:
1178  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1179  for (j = 0; j < nDim; ++j)
1180  {
1181  for (i = 0; i < nvariables ; ++i)
1182  {
1183  // Compute Fwd and Bwd value of ufield of i direction
1184  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1185 
1186  // if Vn >= 0, flux = uFwd, i.e.,
1187  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1188  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1189 
1190  // else if Vn < 0, flux = uBwd, i.e.,
1191  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1192  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1193 
1194  fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1195 
1196  // Imposing weak boundary condition with flux
1197  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1198  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1199  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1200 
1201  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1202  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1203  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1204 
1205  if(fields[0]->GetBndCondExpansions().num_elements())
1206  {
1207  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1208  }
1209 
1210  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1211  // i.e,
1212  // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1213  // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1214 
1215  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1216  // i.e,
1217  // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1218  // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1219 
1220  Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1221  }
1222  }
1223  }
1224 
1225  /**
1226  * @brief Imposes appropriate bcs for the 1st order derivatives
1227  *
1228  */
1229  void DiffusionLFR::v_WeakPenaltyforScalar(
1230  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1231  const int var,
1232  const Array<OneD, const NekDouble> &ufield,
1233  Array<OneD, NekDouble> &penaltyflux)
1234  {
1235  // Variable initialisation
1236  int i, e, id1, id2;
1237  int nBndEdgePts, nBndEdges;
1238  int cnt = 0;
1239  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1240  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1241 
1242  // Extract physical values of the fields at the trace space
1243  Array<OneD, NekDouble > uplus(nTracePts);
1244  fields[var]->ExtractTracePhys(ufield, uplus);
1245 
1246  // Impose boundary conditions
1247  for (i = 0; i < nBndRegions; ++i)
1248  {
1249  // Number of boundary edges related to bcRegion 'i'
1250  nBndEdges = fields[var]->
1251  GetBndCondExpansions()[i]->GetExpSize();
1252 
1253  // Weakly impose boundary conditions by modifying flux values
1254  for (e = 0; e < nBndEdges ; ++e)
1255  {
1256  // Number of points on the expansion
1257  nBndEdgePts = fields[var]->
1258  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1259 
1260  // Offset of the boundary expansion
1261  id1 = fields[var]->
1262  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1263 
1264  // Offset of the trace space related to boundary expansion
1265  id2 = fields[0]->GetTrace()->
1266  GetPhys_Offset(fields[0]->GetTraceMap()->
1267  GetBndCondTraceToGlobalTraceMap(cnt++));
1268 
1269  // Dirichlet bcs ==> uflux = gD
1270  if (fields[var]->GetBndConditions()[i]->
1271  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1272  {
1273  Vmath::Vcopy(nBndEdgePts,
1274  &(fields[var]->
1275  GetBndCondExpansions()[i]->
1276  GetPhys())[id1], 1,
1277  &penaltyflux[id2], 1);
1278  }
1279  // Neumann bcs ==> uflux = u+
1280  else if ((fields[var]->GetBndConditions()[i])->
1281  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1282  {
1283  Vmath::Vcopy(nBndEdgePts,
1284  &uplus[id2], 1,
1285  &penaltyflux[id2], 1);
1286  }
1287  }
1288  }
1289  }
1290 
1291  /**
1292  * @brief Build the numerical flux for the 2nd order derivatives
1293  *
1294  */
1295  void DiffusionLFR::v_NumFluxforVector(
1296  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1297  const Array<OneD, Array<OneD, NekDouble> > &ufield,
1298  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
1299  Array<OneD, Array<OneD, NekDouble> > &qflux)
1300  {
1301  int i, j;
1302  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1303  int nvariables = fields.num_elements();
1304  int nDim = fields[0]->GetCoordim(0);
1305 
1306  NekDouble C11 = 0.0;
1307  Array<OneD, NekDouble > Fwd(nTracePts);
1308  Array<OneD, NekDouble > Bwd(nTracePts);
1309  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1310 
1311  Array<OneD, NekDouble > qFwd (nTracePts);
1312  Array<OneD, NekDouble > qBwd (nTracePts);
1313  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1314  Array<OneD, NekDouble > uterm(nTracePts);
1315 
1316  // Get the normal velocity Vn
1317  for(i = 0; i < nDim; ++i)
1318  {
1319  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1320  Vn, 1, Vn, 1);
1321  }
1322 
1323  // Evaulate upwind flux:
1324  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1325  for (i = 0; i < nvariables; ++i)
1326  {
1327  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1328  for (j = 0; j < nDim; ++j)
1329  {
1330  // Compute Fwd and Bwd value of ufield of jth direction
1331  fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1332 
1333  // if Vn >= 0, flux = uFwd, i.e.,
1334  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1335  // qflux = qBwd = q+
1336  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1337  // qflux = qBwd = q-
1338 
1339  // else if Vn < 0, flux = uBwd, i.e.,
1340  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1341  // qflux = qFwd = q-
1342  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1343  // qflux = qFwd = q+
1344 
1345  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1346 
1347  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1348  qfluxtemp, 1);
1349 
1350  /*
1351  // Generate Stability term = - C11 ( u- - u+ )
1352  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1353 
1354  Vmath::Vsub(nTracePts,
1355  Fwd, 1, Bwd, 1,
1356  uterm, 1);
1357 
1358  Vmath::Smul(nTracePts,
1359  -1.0 * C11, uterm, 1,
1360  uterm, 1);
1361 
1362  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
1363  Vmath::Vadd(nTracePts,
1364  uterm, 1,
1365  qfluxtemp, 1,
1366  qfluxtemp, 1);
1367  */
1368 
1369  // Imposing weak boundary condition with flux
1370  if (fields[0]->GetBndCondExpansions().num_elements())
1371  {
1372  v_WeakPenaltyforVector(fields, i, j, qfield[i][j],
1373  qfluxtemp, C11);
1374  }
1375 
1376  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1377  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1378  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1379  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1380  qflux[i], 1);
1381  }
1382  }
1383  }
1384 
1385 
1386 
1387  /**
1388  * @brief Imposes appropriate bcs for the 2nd order derivatives
1389  *
1390  */
1391  void DiffusionLFR::v_WeakPenaltyforVector(
1392  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1393  const int var,
1394  const int dir,
1395  const Array<OneD, const NekDouble> &qfield,
1396  Array<OneD, NekDouble> &penaltyflux,
1397  NekDouble C11)
1398  {
1399  int i, e, id1, id2;
1400  int nBndEdges, nBndEdgePts;
1401  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1402  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1403 
1404  Array<OneD, NekDouble > uterm(nTracePts);
1405  Array<OneD, NekDouble > qtemp(nTracePts);
1406  int cnt = 0;
1407 
1408  fields[var]->ExtractTracePhys(qfield, qtemp);
1409 
1410  for (i = 0; i < nBndRegions; ++i)
1411  {
1412  nBndEdges = fields[var]->
1413  GetBndCondExpansions()[i]->GetExpSize();
1414 
1415  // Weakly impose boundary conditions by modifying flux values
1416  for (e = 0; e < nBndEdges ; ++e)
1417  {
1418  nBndEdgePts = fields[var]->
1419  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1420 
1421  id1 = fields[var]->
1422  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1423 
1424  id2 = fields[0]->GetTrace()->
1425  GetPhys_Offset(fields[0]->GetTraceMap()->
1426  GetBndCondTraceToGlobalTraceMap(cnt++));
1427 
1428  // For Dirichlet boundary condition:
1429  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1430  if (fields[var]->GetBndConditions()[i]->
1431  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1432  {
1433  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1434  &qtemp[id2], 1, &penaltyflux[id2], 1);
1435  }
1436  // For Neumann boundary condition: qflux = g_N
1437  else if ((fields[var]->GetBndConditions()[i])->
1438  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1439  {
1440  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1441  &(fields[var]->GetBndCondExpansions()[i]->
1442  GetPhys())[id1], 1, &penaltyflux[id2], 1);
1443  }
1444  }
1445  }
1446  }
1447 
1448  /**
1449  * @brief Compute the derivative of the corrective flux for 1D problems.
1450  *
1451  * @param nConvectiveFields Number of fields.
1452  * @param fields Pointer to fields.
1453  * @param flux Volumetric flux in the physical space.
1454  * @param iFlux Numerical interface flux in physical space.
1455  * @param derCFlux Derivative of the corrective flux.
1456  *
1457  */
1458  void DiffusionLFR::v_DerCFlux_1D(
1459  const int nConvectiveFields,
1460  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1461  const Array<OneD, const NekDouble> &flux,
1462  const Array<OneD, const NekDouble> &iFlux,
1463  Array<OneD, NekDouble> &derCFlux)
1464  {
1465  int n;
1466  int nLocalSolutionPts, phys_offset, t_offset;
1467 
1468  Array<OneD, NekDouble> auxArray1, auxArray2;
1469  Array<TwoD, const NekDouble> gmat;
1470  Array<OneD, const NekDouble> jac;
1471 
1472  LibUtilities::PointsKeyVector ptsKeys;
1473  LibUtilities::BasisSharedPtr Basis;
1474  Basis = fields[0]->GetExp(0)->GetBasis(0);
1475 
1476  int nElements = fields[0]->GetExpSize();
1477  int nSolutionPts = fields[0]->GetTotPoints();
1478 
1479  vector<bool> negatedFluxNormal = (boost::static_pointer_cast<MultiRegions::DisContField1D>(fields[0]))->GetNegatedFluxNormal();
1480 
1481  // Arrays to store the derivatives of the correction flux
1482  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1483  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1484 
1485  // Arrays to store the intercell numerical flux jumps
1486  Array<OneD, NekDouble> JumpL(nElements);
1487  Array<OneD, NekDouble> JumpR(nElements);
1488 
1489  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1490  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1491 
1492  for (n = 0; n < nElements; ++n)
1493  {
1494  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1495  phys_offset = fields[0]->GetPhys_Offset(n);
1496 
1497  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1498  NekDouble tmpFluxVertex = 0;
1499  Vmath::Vcopy(nLocalSolutionPts,
1500  &flux[phys_offset], 1,
1501  &tmparrayX1[0], 1);
1502 
1503  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1504  tmpFluxVertex);
1505 
1506  t_offset = fields[0]->GetTrace()
1507  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1508 
1509  if(negatedFluxNormal[2*n])
1510  {
1511  JumpL[n] = iFlux[t_offset] - tmpFluxVertex;
1512  }
1513  else
1514  {
1515  JumpL[n] = -iFlux[t_offset] - tmpFluxVertex;
1516  }
1517 
1518  t_offset = fields[0]->GetTrace()
1519  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1520 
1521  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1522  tmpFluxVertex);
1523  if(negatedFluxNormal[2*n+1])
1524  {
1525  JumpR[n] = -iFlux[t_offset] - tmpFluxVertex;
1526  }
1527  else
1528  {
1529  JumpR[n] = iFlux[t_offset] - tmpFluxVertex;
1530  }
1531  }
1532 
1533  for (n = 0; n < nElements; ++n)
1534  {
1535  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1536  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1537  phys_offset = fields[0]->GetPhys_Offset(n);
1538  jac = fields[0]->GetExp(n)
1539  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1540  ->GetMetricInfo()->GetJac(ptsKeys);
1541 
1542  JumpL[n] = JumpL[n] * jac[0];
1543  JumpR[n] = JumpR[n] * jac[0];
1544 
1545  // Left jump multiplied by left derivative of C function
1546  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1547  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1548 
1549  // Right jump multiplied by right derivative of C function
1550  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1551  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1552 
1553  // Assembling derivative of the correction flux
1554  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1555  &derCFlux[phys_offset], 1);
1556  }
1557  }
1558 
1559  /**
1560  * @brief Compute the derivative of the corrective flux wrt a given
1561  * coordinate for 2D problems.
1562  *
1563  * There could be a bug for deformed elements since first derivatives
1564  * are not exactly the same results of DiffusionLDG as expected
1565  *
1566  * @param nConvectiveFields Number of fields.
1567  * @param fields Pointer to fields.
1568  * @param flux Volumetric flux in the physical space.
1569  * @param numericalFlux Numerical interface flux in physical space.
1570  * @param derCFlux Derivative of the corrective flux.
1571  *
1572  * \todo: Switch on shapes eventually here.
1573  */
1574  void DiffusionLFR::v_DerCFlux_2D(
1575  const int nConvectiveFields,
1576  const int direction,
1577  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1578  const Array<OneD, const NekDouble> &flux,
1579  const Array<OneD, NekDouble> &iFlux,
1580  Array<OneD, NekDouble> &derCFlux)
1581  {
1582  int n, e, i, j, cnt;
1583 
1584  Array<TwoD, const NekDouble> gmat;
1585  Array<OneD, const NekDouble> jac;
1586 
1587  int nElements = fields[0]->GetExpSize();
1588 
1589  int trace_offset, phys_offset;
1590  int nLocalSolutionPts;
1591  int nquad0, nquad1;
1592  int nEdgePts;
1593 
1594  LibUtilities::PointsKeyVector ptsKeys;
1595  Array<OneD, NekDouble> auxArray1, auxArray2;
1596  Array<OneD, LibUtilities::BasisSharedPtr> base;
1597  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1598  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1599 
1600  // Loop on the elements
1601  for (n = 0; n < nElements; ++n)
1602  {
1603  // Offset of the element on the global vector
1604  phys_offset = fields[0]->GetPhys_Offset(n);
1605  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1606  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1607 
1608  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1609  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1610 
1611  base = fields[0]->GetExp(n)->GetBase();
1612  nquad0 = base[0]->GetNumPoints();
1613  nquad1 = base[1]->GetNumPoints();
1614 
1615  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1616  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1617  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1618  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1619 
1620  // Loop on the edges
1621  for (e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1622  {
1623  // Number of edge points of edge 'e'
1624  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1625 
1626  // Array for storing volumetric fluxes on each edge
1627  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1628 
1629  // Offset of the trace space correspondent to edge 'e'
1630  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1631  elmtToTrace[n][e]->GetElmtId());
1632 
1633  // Get the normals of edge 'e'
1634  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1635  //fields[0]->GetExp(n)->GetEdgeNormal(e);
1636 
1637  // Extract the edge values of the volumetric fluxes
1638  // on edge 'e' and order them accordingly to the order
1639  // of the trace space
1640  fields[0]->GetExp(n)->GetEdgePhysVals(e, elmtToTrace[n][e],
1641  flux + phys_offset,
1642  auxArray1 = tmparray);
1643 
1644  // Compute the fluxJumps per each edge 'e' and each
1645  // flux point
1646  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1647  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1648  &tmparray[0], 1, &fluxJumps[0], 1);
1649 
1650  // Check the ordering of the fluxJumps and reverse
1651  // it in case of backward definition of edge 'e'
1652  if (fields[0]->GetExp(n)->GetEorient(e) ==
1653  StdRegions::eBackwards)
1654  {
1655  Vmath::Reverse(nEdgePts,
1656  &fluxJumps[0], 1,
1657  &fluxJumps[0], 1);
1658  }
1659 
1660  // Deformed elements
1661  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1662  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1663  == SpatialDomains::eDeformed)
1664  {
1665  // Extract the Jacobians along edge 'e'
1666  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1667 
1668  fields[0]->GetExp(n)->GetEdgePhysVals(
1669  e, elmtToTrace[n][e],
1670  jac, auxArray1 = jacEdge);
1671 
1672  // Check the ordering of the fluxJumps and reverse
1673  // it in case of backward definition of edge 'e'
1674  if (fields[0]->GetExp(n)->GetEorient(e) ==
1675  StdRegions::eBackwards)
1676  {
1677  Vmath::Reverse(nEdgePts,
1678  &jacEdge[0], 1,
1679  &jacEdge[0], 1);
1680  }
1681 
1682  // Multiply the fluxJumps by the edge 'e' Jacobians
1683  // to bring the fluxJumps into the standard space
1684  for (j = 0; j < nEdgePts; j++)
1685  {
1686  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1687  }
1688  }
1689  // Non-deformed elements
1690  else
1691  {
1692  // Multiply the fluxJumps by the edge 'e' Jacobians
1693  // to bring the fluxJumps into the standard space
1694  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1695  fluxJumps, 1);
1696  }
1697 
1698  // Multiply jumps by derivatives of the correction functions
1699  // All the quntities at this level should be defined into
1700  // the standard space
1701  switch (e)
1702  {
1703  case 0:
1704  for (i = 0; i < nquad0; ++i)
1705  {
1706  // Multiply fluxJumps by Q factors
1707  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1708 
1709  for (j = 0; j < nquad1; ++j)
1710  {
1711  cnt = i + j*nquad0;
1712  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1713  }
1714  }
1715  break;
1716  case 1:
1717  for (i = 0; i < nquad1; ++i)
1718  {
1719  // Multiply fluxJumps by Q factors
1720  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1721 
1722  for (j = 0; j < nquad0; ++j)
1723  {
1724  cnt = (nquad0)*i + j;
1725  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1726  }
1727  }
1728  break;
1729  case 2:
1730  for (i = 0; i < nquad0; ++i)
1731  {
1732  // Multiply fluxJumps by Q factors
1733  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1734 
1735  for (j = 0; j < nquad1; ++j)
1736  {
1737  cnt = (nquad0 - 1) + j*nquad0 - i;
1738  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1739  }
1740  }
1741  break;
1742  case 3:
1743  for (i = 0; i < nquad1; ++i)
1744  {
1745  // Multiply fluxJumps by Q factors
1746  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1747  for (j = 0; j < nquad0; ++j)
1748  {
1749  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1750  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1751  }
1752  }
1753  break;
1754 
1755  default:
1756  ASSERTL0(false, "edge value (< 3) is out of range");
1757  break;
1758  }
1759  }
1760 
1761  // Summing the component of the flux for calculating the
1762  // derivatives per each direction
1763  if (direction == 0)
1764  {
1765  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1766  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1767  }
1768  else if (direction == 1)
1769  {
1770  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1771  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1772  }
1773  }
1774  }
1775 
1776  /**
1777  * @brief Compute the divergence of the corrective flux for 2D problems.
1778  *
1779  * @param nConvectiveFields Number of fields.
1780  * @param fields Pointer to fields.
1781  * @param fluxX1 X1 - volumetric flux in physical space.
1782  * @param fluxX2 X2 - volumetric flux in physical space.
1783  * @param numericalFlux Interface flux in physical space.
1784  * @param divCFlux Divergence of the corrective flux for
1785  * 2D problems.
1786  *
1787  * \todo: Switch on shapes eventually here.
1788  */
1789  void DiffusionLFR::v_DivCFlux_2D(
1790  const int nConvectiveFields,
1791  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1792  const Array<OneD, const NekDouble> &fluxX1,
1793  const Array<OneD, const NekDouble> &fluxX2,
1794  const Array<OneD, const NekDouble> &numericalFlux,
1795  Array<OneD, NekDouble> &divCFlux)
1796  {
1797  int n, e, i, j, cnt;
1798 
1799  int nElements = fields[0]->GetExpSize();
1800 
1801  int nLocalSolutionPts;
1802  int nEdgePts;
1803  int trace_offset;
1804  int phys_offset;
1805  int nquad0;
1806  int nquad1;
1807 
1808  Array<OneD, NekDouble> auxArray1, auxArray2;
1809  Array<OneD, LibUtilities::BasisSharedPtr> base;
1810 
1811  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1812  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1813 
1814  // Loop on the elements
1815  for(n = 0; n < nElements; ++n)
1816  {
1817  // Offset of the element on the global vector
1818  phys_offset = fields[0]->GetPhys_Offset(n);
1819  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1820 
1821  base = fields[0]->GetExp(n)->GetBase();
1822  nquad0 = base[0]->GetNumPoints();
1823  nquad1 = base[1]->GetNumPoints();
1824 
1825  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1826  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1827  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1828  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1829 
1830  // Loop on the edges
1831  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1832  {
1833  // Number of edge points of edge e
1834  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1835 
1836  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1837  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1838  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1839  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1840  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1841 
1842  // Offset of the trace space correspondent to edge e
1843  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1844  elmtToTrace[n][e]->GetElmtId());
1845 
1846  // Get the normals of edge e
1847  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1848  fields[0]->GetExp(n)->GetEdgeNormal(e);
1849 
1850  // Extract the edge values of flux-x on edge e and order
1851  // them accordingly to the order of the trace space
1852  fields[0]->GetExp(n)->GetEdgePhysVals(
1853  e, elmtToTrace[n][e],
1854  fluxX1 + phys_offset,
1855  auxArray1 = tmparrayX1);
1856 
1857  // Extract the edge values of flux-y on edge e and order
1858  // them accordingly to the order of the trace space
1859  fields[0]->GetExp(n)->GetEdgePhysVals(
1860  e, elmtToTrace[n][e],
1861  fluxX2 + phys_offset,
1862  auxArray1 = tmparrayX2);
1863 
1864  // Multiply the edge components of the flux by the normal
1865  for (i = 0; i < nEdgePts; ++i)
1866  {
1867  fluxN[i] =
1868  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
1869  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
1870  }
1871 
1872  // Subtract to the Riemann flux the discontinuous flux
1873  Vmath::Vsub(nEdgePts,
1874  &numericalFlux[trace_offset], 1,
1875  &fluxN[0], 1, &fluxJumps[0], 1);
1876 
1877  // Check the ordering of the jump vectors
1878  if (fields[0]->GetExp(n)->GetEorient(e) ==
1879  StdRegions::eBackwards)
1880  {
1881  Vmath::Reverse(nEdgePts,
1882  auxArray1 = fluxJumps, 1,
1883  auxArray2 = fluxJumps, 1);
1884  }
1885 
1886  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1887  -1.0 : 1.0;
1888 
1889  for (i = 0; i < nEdgePts; ++i)
1890  {
1891  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1892  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1893  {
1894  fluxJumps[i] = -fluxJumps[i];
1895  }
1896  }
1897 
1898  // Multiply jumps by derivatives of the correction functions
1899  switch (e)
1900  {
1901  case 0:
1902  for (i = 0; i < nquad0; ++i)
1903  {
1904  // Multiply fluxJumps by Q factors
1905  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1906 
1907  for (j = 0; j < nquad1; ++j)
1908  {
1909  cnt = i + j*nquad0;
1910  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1911  }
1912  }
1913  break;
1914  case 1:
1915  for (i = 0; i < nquad1; ++i)
1916  {
1917  // Multiply fluxJumps by Q factors
1918  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1919 
1920  for (j = 0; j < nquad0; ++j)
1921  {
1922  cnt = (nquad0)*i + j;
1923  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1924  }
1925  }
1926  break;
1927  case 2:
1928  for (i = 0; i < nquad0; ++i)
1929  {
1930  // Multiply fluxJumps by Q factors
1931  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1932 
1933  for (j = 0; j < nquad1; ++j)
1934  {
1935  cnt = (nquad0 - 1) + j*nquad0 - i;
1936  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1937  }
1938  }
1939  break;
1940  case 3:
1941  for (i = 0; i < nquad1; ++i)
1942  {
1943  // Multiply fluxJumps by Q factors
1944  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1945  for (j = 0; j < nquad0; ++j)
1946  {
1947  cnt = (nquad0*nquad1 - nquad0) + j - i*nquad0;
1948  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1949  }
1950  }
1951  break;
1952 
1953  default:
1954  ASSERTL0(false,"edge value (< 3) is out of range");
1955  break;
1956  }
1957  }
1958 
1959  // Sum each edge contribution
1960  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1961  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1962 
1963  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1964  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1965 
1966  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1967  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1968  }
1969  }
1970 
1971  /**
1972  * @brief Compute the divergence of the corrective flux for 2D problems
1973  * where POINTSTYPE="GaussGaussLegendre"
1974  *
1975  * @param nConvectiveFields Number of fields.
1976  * @param fields Pointer to fields.
1977  * @param fluxX1 X1-volumetric flux in physical space.
1978  * @param fluxX2 X2-volumetric flux in physical space.
1979  * @param numericalFlux Interface flux in physical space.
1980  * @param divCFlux Divergence of the corrective flux.
1981  *
1982  * \todo: Switch on shapes eventually here.
1983  */
1984 
1985  void DiffusionLFR::v_DivCFlux_2D_Gauss(
1986  const int nConvectiveFields,
1987  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
1988  const Array<OneD, const NekDouble> &fluxX1,
1989  const Array<OneD, const NekDouble> &fluxX2,
1990  const Array<OneD, const NekDouble> &numericalFlux,
1991  Array<OneD, NekDouble> &divCFlux)
1992  {
1993  int n, e, i, j, cnt;
1994 
1995  int nElements = fields[0]->GetExpSize();
1996  int nLocalSolutionPts;
1997  int nEdgePts;
1998  int trace_offset;
1999  int phys_offset;
2000  int nquad0;
2001  int nquad1;
2002 
2003  Array<OneD, NekDouble> auxArray1, auxArray2;
2004  Array<OneD, LibUtilities::BasisSharedPtr> base;
2005 
2006  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
2007  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2008 
2009  // Loop on the elements
2010  for(n = 0; n < nElements; ++n)
2011  {
2012  // Offset of the element on the global vector
2013  phys_offset = fields[0]->GetPhys_Offset(n);
2014  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2015 
2016  base = fields[0]->GetExp(n)->GetBase();
2017  nquad0 = base[0]->GetNumPoints();
2018  nquad1 = base[1]->GetNumPoints();
2019 
2020  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2021  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2022  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2023  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2024 
2025  // Loop on the edges
2026  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
2027  {
2028  // Number of edge points of edge e
2029  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
2030 
2031  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2032  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2033  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2034  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2035  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2036 
2037  // Offset of the trace space correspondent to edge e
2038  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2039  elmtToTrace[n][e]->GetElmtId());
2040 
2041  // Get the normals of edge e
2042  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2043  fields[0]->GetExp(n)->GetEdgeNormal(e);
2044 
2045  // Extract the trasformed normal flux at each edge
2046  switch (e)
2047  {
2048  case 0:
2049  // Extract the edge values of transformed flux-y on
2050  // edge e and order them accordingly to the order of
2051  // the trace space
2052  fields[0]->GetExp(n)->GetEdgePhysVals(
2053  e, elmtToTrace[n][e],
2054  fluxX2 + phys_offset,
2055  auxArray1 = fluxN_D);
2056 
2057  Vmath::Neg (nEdgePts, fluxN_D, 1);
2058 
2059  // Extract the physical Rieamann flux at each edge
2060  Vmath::Vcopy(nEdgePts,
2061  &numericalFlux[trace_offset], 1,
2062  &fluxN[0], 1);
2063 
2064  // Check the ordering of vectors
2065  if (fields[0]->GetExp(n)->GetEorient(e) ==
2066  StdRegions::eBackwards)
2067  {
2068  Vmath::Reverse(nEdgePts,
2069  auxArray1 = fluxN, 1,
2070  auxArray2 = fluxN, 1);
2071 
2072  Vmath::Reverse(nEdgePts,
2073  auxArray1 = fluxN_D, 1,
2074  auxArray2 = fluxN_D, 1);
2075  }
2076 
2077  // Transform Riemann Fluxes in the standard element
2078  for (i = 0; i < nquad0; ++i)
2079  {
2080  // Multiply Riemann Flux by Q factors
2081  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2082  }
2083 
2084  for (i = 0; i < nEdgePts; ++i)
2085  {
2086  if (m_traceNormals[0][trace_offset+i] !=
2087  normals[0][i] ||
2088  m_traceNormals[1][trace_offset+i] !=
2089  normals[1][i])
2090  {
2091  fluxN_R[i] = -fluxN_R[i];
2092  }
2093  }
2094 
2095  // Subtract to the Riemann flux the discontinuous
2096  // flux
2097  Vmath::Vsub(nEdgePts,
2098  &fluxN_R[0], 1,
2099  &fluxN_D[0], 1, &fluxJumps[0], 1);
2100 
2101  // Multiplicate the flux jumps for the correction
2102  // function
2103  for (i = 0; i < nquad0; ++i)
2104  {
2105  for (j = 0; j < nquad1; ++j)
2106  {
2107  cnt = i + j*nquad0;
2108  divCFluxE0[cnt] =
2109  -fluxJumps[i] * m_dGL_xi2[n][j];
2110  }
2111  }
2112  break;
2113  case 1:
2114  // Extract the edge values of transformed flux-x on
2115  // edge e and order them accordingly to the order of
2116  // the trace space
2117  fields[0]->GetExp(n)->GetEdgePhysVals(
2118  e, elmtToTrace[n][e],
2119  fluxX1 + phys_offset,
2120  auxArray1 = fluxN_D);
2121 
2122  // Extract the physical Rieamann flux at each edge
2123  Vmath::Vcopy(nEdgePts,
2124  &numericalFlux[trace_offset], 1,
2125  &fluxN[0], 1);
2126 
2127  // Check the ordering of vectors
2128  if (fields[0]->GetExp(n)->GetEorient(e) ==
2129  StdRegions::eBackwards)
2130  {
2131  Vmath::Reverse(nEdgePts,
2132  auxArray1 = fluxN, 1,
2133  auxArray2 = fluxN, 1);
2134 
2135  Vmath::Reverse(nEdgePts,
2136  auxArray1 = fluxN_D, 1,
2137  auxArray2 = fluxN_D, 1);
2138  }
2139 
2140  // Transform Riemann Fluxes in the standard element
2141  for (i = 0; i < nquad1; ++i)
2142  {
2143  // Multiply Riemann Flux by Q factors
2144  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2145  }
2146 
2147  for (i = 0; i < nEdgePts; ++i)
2148  {
2149  if (m_traceNormals[0][trace_offset+i] !=
2150  normals[0][i] ||
2151  m_traceNormals[1][trace_offset+i] !=
2152  normals[1][i])
2153  {
2154  fluxN_R[i] = -fluxN_R[i];
2155  }
2156  }
2157 
2158  // Subtract to the Riemann flux the discontinuous
2159  // flux
2160  Vmath::Vsub(nEdgePts,
2161  &fluxN_R[0], 1,
2162  &fluxN_D[0], 1, &fluxJumps[0], 1);
2163 
2164  // Multiplicate the flux jumps for the correction
2165  // function
2166  for (i = 0; i < nquad1; ++i)
2167  {
2168  for (j = 0; j < nquad0; ++j)
2169  {
2170  cnt = (nquad0)*i + j;
2171  divCFluxE1[cnt] =
2172  fluxJumps[i] * m_dGR_xi1[n][j];
2173  }
2174  }
2175  break;
2176  case 2:
2177 
2178  // Extract the edge values of transformed flux-y on
2179  // edge e and order them accordingly to the order of
2180  // the trace space
2181 
2182  fields[0]->GetExp(n)->GetEdgePhysVals(
2183  e, elmtToTrace[n][e],
2184  fluxX2 + phys_offset,
2185  auxArray1 = fluxN_D);
2186 
2187  // Extract the physical Rieamann flux at each edge
2188  Vmath::Vcopy(nEdgePts,
2189  &numericalFlux[trace_offset], 1,
2190  &fluxN[0], 1);
2191 
2192  // Check the ordering of vectors
2193  if (fields[0]->GetExp(n)->GetEorient(e) ==
2194  StdRegions::eBackwards)
2195  {
2196  Vmath::Reverse(nEdgePts,
2197  auxArray1 = fluxN, 1,
2198  auxArray2 = fluxN, 1);
2199 
2200  Vmath::Reverse(nEdgePts,
2201  auxArray1 = fluxN_D, 1,
2202  auxArray2 = fluxN_D, 1);
2203  }
2204 
2205  // Transform Riemann Fluxes in the standard element
2206  for (i = 0; i < nquad0; ++i)
2207  {
2208  // Multiply Riemann Flux by Q factors
2209  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2210  }
2211 
2212  for (i = 0; i < nEdgePts; ++i)
2213  {
2214  if (m_traceNormals[0][trace_offset+i] !=
2215  normals[0][i] ||
2216  m_traceNormals[1][trace_offset+i] !=
2217  normals[1][i])
2218  {
2219  fluxN_R[i] = -fluxN_R[i];
2220  }
2221  }
2222 
2223  // Subtract to the Riemann flux the discontinuous
2224  // flux
2225 
2226  Vmath::Vsub(nEdgePts,
2227  &fluxN_R[0], 1,
2228  &fluxN_D[0], 1, &fluxJumps[0], 1);
2229 
2230  // Multiplicate the flux jumps for the correction
2231  // function
2232  for (i = 0; i < nquad0; ++i)
2233  {
2234  for (j = 0; j < nquad1; ++j)
2235  {
2236  cnt = (nquad0 - 1) + j*nquad0 - i;
2237  divCFluxE2[cnt] =
2238  fluxJumps[i] * m_dGR_xi2[n][j];
2239  }
2240  }
2241  break;
2242  case 3:
2243  // Extract the edge values of transformed flux-x on
2244  // edge e and order them accordingly to the order of
2245  // the trace space
2246 
2247  fields[0]->GetExp(n)->GetEdgePhysVals(
2248  e, elmtToTrace[n][e],
2249  fluxX1 + phys_offset,
2250  auxArray1 = fluxN_D);
2251  Vmath::Neg (nEdgePts, fluxN_D, 1);
2252 
2253  // Extract the physical Rieamann flux at each edge
2254  Vmath::Vcopy(nEdgePts,
2255  &numericalFlux[trace_offset], 1,
2256  &fluxN[0], 1);
2257 
2258  // Check the ordering of vectors
2259  if (fields[0]->GetExp(n)->GetEorient(e) ==
2260  StdRegions::eBackwards)
2261  {
2262  Vmath::Reverse(nEdgePts,
2263  auxArray1 = fluxN, 1,
2264  auxArray2 = fluxN, 1);
2265 
2266  Vmath::Reverse(nEdgePts,
2267  auxArray1 = fluxN_D, 1,
2268  auxArray2 = fluxN_D, 1);
2269  }
2270 
2271  // Transform Riemann Fluxes in the standard element
2272  for (i = 0; i < nquad1; ++i)
2273  {
2274  // Multiply Riemann Flux by Q factors
2275  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2276  }
2277 
2278  for (i = 0; i < nEdgePts; ++i)
2279  {
2280  if (m_traceNormals[0][trace_offset+i] !=
2281  normals[0][i] ||
2282  m_traceNormals[1][trace_offset+i] !=
2283  normals[1][i])
2284  {
2285  fluxN_R[i] = -fluxN_R[i];
2286  }
2287  }
2288 
2289  // Subtract to the Riemann flux the discontinuous
2290  // flux
2291 
2292  Vmath::Vsub(nEdgePts,
2293  &fluxN_R[0], 1,
2294  &fluxN_D[0], 1, &fluxJumps[0], 1);
2295 
2296  // Multiplicate the flux jumps for the correction
2297  // function
2298  for (i = 0; i < nquad1; ++i)
2299  {
2300  for (j = 0; j < nquad0; ++j)
2301  {
2302  cnt = (nquad0*nquad1 - nquad0) + j
2303  - i*nquad0;
2304  divCFluxE3[cnt] =
2305  -fluxJumps[i] * m_dGL_xi1[n][j];
2306  }
2307  }
2308  break;
2309  default:
2310  ASSERTL0(false,"edge value (< 3) is out of range");
2311  break;
2312  }
2313  }
2314 
2315 
2316  // Sum each edge contribution
2317  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2318  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2319 
2320  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2321  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2322 
2323  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2324  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2325  }
2326  }
2327 
2328  }// close namespace SolverUtils
2329 }// close namespace nektar++
Generated on Sun Mar 15 2015 00:11:30 for Nektar++ by   doxygen 1.8.1.2