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