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