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