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