Nektar++
Static Public Member Functions | Public Attributes | Static Public Attributes | Protected Member Functions | Protected Attributes | List of all members
Nektar::SolverUtils::DiffusionLFR Class Reference

#include <DiffusionLFR.h>

Inheritance diagram for Nektar::SolverUtils::DiffusionLFR:
[legend]

Static Public Member Functions

static DiffusionSharedPtr create (std::string diffType)
 

Public Attributes

Array< OneD, NekDouble > m_jac
 
Array< OneD, Array< OneD, NekDouble > > m_gmat
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
 
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
 
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
 
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
 
DNekMatSharedPtr m_Ixm
 
DNekMatSharedPtr m_Ixp
 

Static Public Attributes

static std::string type []
 

Protected Member Functions

 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. More...
 
virtual void v_InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 Initiliase DiffusionLFR objects and store them before starting the time-stepping. More...
 
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 the fluxes at the interfaces of each element from the physical space to the standard space). More...
 
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–72. More...
 
virtual void v_Diffuse (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 Calculate FR Diffusion for the linear problems using an LDG interface flux. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 
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". More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::Diffusion
virtual void v_SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual Array< OneD, Array< OneD, Array< OneD, NekDouble > > > & v_GetFluxTensor ()
 

Protected Attributes

Array< OneD, Array< OneD, NekDouble > > m_traceNormals
 
LibUtilities::SessionReaderSharedPtr m_session
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
 
Array< OneD, Array< OneD, NekDouble > > m_IF2
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
 
Array< OneD, Array< OneD, NekDouble > > m_divFD
 
Array< OneD, Array< OneD, NekDouble > > m_divFC
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
 
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
 
std::string m_diffType
 
- Protected Attributes inherited from Nektar::SolverUtils::Diffusion
DiffusionFluxVecCB m_fluxVector
 
DiffusionFluxVecCBNS m_fluxVectorNS
 
DiffusionFluxPenaltyNS m_fluxPenaltyNS
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Diffusion
virtual SOLVER_UTILS_EXPORT ~Diffusion ()
 
SOLVER_UTILS_EXPORT void InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 
SOLVER_UTILS_EXPORT void Diffuse (const std::size_t nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 
SOLVER_UTILS_EXPORT void FluxVec (Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &fluxvector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVector (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxVector (DiffusionFluxVecCB fluxVector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVectorNS (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxVectorNS (DiffusionFluxVecCBNS fluxVector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxPenaltyNS (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxPenaltyNS (DiffusionFluxPenaltyNS flux)
 
void SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual Array< OneD, Array< OneD, Array< OneD, NekDouble > > > & GetFluxTensor ()
 

Detailed Description

Definition at line 44 of file DiffusionLFR.h.

Constructor & Destructor Documentation

◆ DiffusionLFR()

Nektar::SolverUtils::DiffusionLFR::DiffusionLFR ( std::string  diffType)
protected

DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.

Todo:
Extension to triangles, tetrahedra and other shapes. (Long term objective)

Definition at line 71 of file DiffusionLFR.cpp.

Referenced by create().

71  :m_diffType(diffType)
72  {
73  }

Member Function Documentation

◆ create()

static DiffusionSharedPtr Nektar::SolverUtils::DiffusionLFR::create ( std::string  diffType)
inlinestatic

Definition at line 47 of file DiffusionLFR.h.

References DiffusionLFR().

48  {
49  return DiffusionSharedPtr(new DiffusionLFR(diffType));
50  }
Nektar::SolverUtils::DiffusionSharedPtr
std::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:183
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:71

◆ v_DerCFlux_1D()

void Nektar::SolverUtils::DiffusionLFR::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 
)
protectedvirtual

Compute the derivative of the corrective flux for 1D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxVolumetric flux in the physical space.
iFluxNumerical interface flux in physical space.
derCFluxDerivative of the corrective flux.

Definition at line 1474 of file DiffusionLFR.cpp.

References Nektar::LocalRegions::Expansion::GetMetricInfo(), m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by v_Diffuse().

1480  {
1481  boost::ignore_unused(nConvectiveFields);
1482 
1483  int n;
1484  int nLocalSolutionPts, phys_offset, t_offset;
1485 
1486  Array<OneD, NekDouble> auxArray1, auxArray2;
1487  Array<TwoD, const NekDouble> gmat;
1488  Array<OneD, const NekDouble> jac;
1489 
1490  LibUtilities::PointsKeyVector ptsKeys;
1491  LibUtilities::BasisSharedPtr Basis;
1492  Basis = fields[0]->GetExp(0)->GetBasis(0);
1493 
1494  int nElements = fields[0]->GetExpSize();
1495  int nSolutionPts = fields[0]->GetTotPoints();
1496 
1497  vector<bool> negatedFluxNormal =
1498  std::static_pointer_cast<MultiRegions::DisContField1D>(
1499  fields[0])->GetNegatedFluxNormal();
1500 
1501  // Arrays to store the derivatives of the correction flux
1502  Array<OneD, NekDouble> DCL(nSolutionPts/nElements, 0.0);
1503  Array<OneD, NekDouble> DCR(nSolutionPts/nElements, 0.0);
1504 
1505  // Arrays to store the intercell numerical flux jumps
1506  Array<OneD, NekDouble> JumpL(nElements);
1507  Array<OneD, NekDouble> JumpR(nElements);
1508 
1509  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1510  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1511 
1512  for (n = 0; n < nElements; ++n)
1513  {
1514  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1515  phys_offset = fields[0]->GetPhys_Offset(n);
1516 
1517  Array<OneD, NekDouble> tmparrayX1(nLocalSolutionPts, 0.0);
1518  NekDouble tmpFluxVertex = 0;
1519  Vmath::Vcopy(nLocalSolutionPts,
1520  &flux[phys_offset], 1,
1521  &tmparrayX1[0], 1);
1522 
1523  fields[0]->GetExp(n)->GetVertexPhysVals(0, tmparrayX1,
1524  tmpFluxVertex);
1525 
1526  t_offset = fields[0]->GetTrace()
1527  ->GetPhys_Offset(elmtToTrace[n][0]->GetElmtId());
1528 
1529  if(negatedFluxNormal[2*n])
1530  {
1531  JumpL[n] = iFlux[t_offset] - tmpFluxVertex;
1532  }
1533  else
1534  {
1535  JumpL[n] = -iFlux[t_offset] - tmpFluxVertex;
1536  }
1537 
1538  t_offset = fields[0]->GetTrace()
1539  ->GetPhys_Offset(elmtToTrace[n][1]->GetElmtId());
1540 
1541  fields[0]->GetExp(n)->GetVertexPhysVals(1, tmparrayX1,
1542  tmpFluxVertex);
1543  if(negatedFluxNormal[2*n+1])
1544  {
1545  JumpR[n] = -iFlux[t_offset] - tmpFluxVertex;
1546  }
1547  else
1548  {
1549  JumpR[n] = iFlux[t_offset] - tmpFluxVertex;
1550  }
1551  }
1552 
1553  for (n = 0; n < nElements; ++n)
1554  {
1555  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1556  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1557  phys_offset = fields[0]->GetPhys_Offset(n);
1558  jac = fields[0]->GetExp(n)
1559  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
1560  ->GetMetricInfo()->GetJac(ptsKeys);
1561 
1562  JumpL[n] = JumpL[n] * jac[0];
1563  JumpR[n] = JumpR[n] * jac[0];
1564 
1565  // Left jump multiplied by left derivative of C function
1566  Vmath::Smul(nLocalSolutionPts, JumpL[n],
1567  &m_dGL_xi1[n][0], 1, &DCL[0], 1);
1568 
1569  // Right jump multiplied by right derivative of C function
1570  Vmath::Smul(nLocalSolutionPts, JumpR[n],
1571  &m_dGR_xi1[n][0], 1, &DCR[0], 1);
1572 
1573  // Assembling derivative of the correction flux
1574  Vmath::Vadd(nLocalSolutionPts, &DCL[0], 1, &DCR[0], 1,
1575  &derCFlux[phys_offset], 1);
1576  }
1577  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:74
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

◆ v_DerCFlux_2D()

void Nektar::SolverUtils::DiffusionLFR::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 
)
protectedvirtual

Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.

There could be a bug for deformed elements since first derivatives are not exactly the same results of DiffusionLDG as expected

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxVolumetric flux in the physical space.
numericalFluxNumerical interface flux in physical space.
derCFluxDerivative of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 1594 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, Nektar::SpatialDomains::eDeformed, Nektar::LocalRegions::Expansion::GetMetricInfo(), m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, Vmath::Reverse(), Vmath::Smul(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by v_Diffuse().

1601  {
1602  boost::ignore_unused(nConvectiveFields);
1603 
1604  int n, e, i, j, cnt;
1605 
1606  Array<TwoD, const NekDouble> gmat;
1607  Array<OneD, const NekDouble> jac;
1608 
1609  int nElements = fields[0]->GetExpSize();
1610 
1611  int trace_offset, phys_offset;
1612  int nLocalSolutionPts;
1613  int nquad0, nquad1;
1614  int nEdgePts;
1615 
1616  LibUtilities::PointsKeyVector ptsKeys;
1617  Array<OneD, NekDouble> auxArray1, auxArray2;
1618  Array<OneD, LibUtilities::BasisSharedPtr> base;
1619  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1620  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1621 
1622  // Loop on the elements
1623  for (n = 0; n < nElements; ++n)
1624  {
1625  // Offset of the element on the global vector
1626  phys_offset = fields[0]->GetPhys_Offset(n);
1627  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1628  ptsKeys = fields[0]->GetExp(n)->GetPointsKeys();
1629 
1630  jac = fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1631  ->GetGeom2D()->GetMetricInfo()->GetJac(ptsKeys);
1632 
1633  base = fields[0]->GetExp(n)->GetBase();
1634  nquad0 = base[0]->GetNumPoints();
1635  nquad1 = base[1]->GetNumPoints();
1636 
1637  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1638  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1639  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1640  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1641 
1642  // Loop on the edges
1643  for (e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1644  {
1645  // Number of edge points of edge 'e'
1646  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1647 
1648  // Array for storing volumetric fluxes on each edge
1649  Array<OneD, NekDouble> tmparray(nEdgePts, 0.0);
1650 
1651  // Offset of the trace space correspondent to edge 'e'
1652  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1653  elmtToTrace[n][e]->GetElmtId());
1654 
1655  // Get the normals of edge 'e'
1656  //const Array<OneD, const Array<OneD, NekDouble> > &normals =
1657  //fields[0]->GetExp(n)->GetEdgeNormal(e);
1658 
1659  // Extract the edge values of the volumetric fluxes
1660  // on edge 'e' and order them accordingly to the order
1661  // of the trace space
1662  fields[0]->GetExp(n)->GetEdgePhysVals(e, elmtToTrace[n][e],
1663  flux + phys_offset,
1664  auxArray1 = tmparray);
1665 
1666  // Compute the fluxJumps per each edge 'e' and each
1667  // flux point
1668  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1669  Vmath::Vsub(nEdgePts, &iFlux[trace_offset], 1,
1670  &tmparray[0], 1, &fluxJumps[0], 1);
1671 
1672  // Check the ordering of the fluxJumps and reverse
1673  // it in case of backward definition of edge 'e'
1674  if (fields[0]->GetExp(n)->GetEorient(e) ==
1675  StdRegions::eBackwards)
1676  {
1677  Vmath::Reverse(nEdgePts,
1678  &fluxJumps[0], 1,
1679  &fluxJumps[0], 1);
1680  }
1681 
1682  // Deformed elements
1683  if (fields[0]->GetExp(n)->as<LocalRegions::Expansion2D>()
1684  ->GetGeom2D()->GetMetricInfo()->GetGtype()
1685  == SpatialDomains::eDeformed)
1686  {
1687  // Extract the Jacobians along edge 'e'
1688  Array<OneD, NekDouble> jacEdge(nEdgePts, 0.0);
1689 
1690  fields[0]->GetExp(n)->GetEdgePhysVals(
1691  e, elmtToTrace[n][e],
1692  jac, auxArray1 = jacEdge);
1693 
1694  // Check the ordering of the fluxJumps and reverse
1695  // it in case of backward definition of edge 'e'
1696  if (fields[0]->GetExp(n)->GetEorient(e) ==
1697  StdRegions::eBackwards)
1698  {
1699  Vmath::Reverse(nEdgePts,
1700  &jacEdge[0], 1,
1701  &jacEdge[0], 1);
1702  }
1703 
1704  // Multiply the fluxJumps by the edge 'e' Jacobians
1705  // to bring the fluxJumps into the standard space
1706  for (j = 0; j < nEdgePts; j++)
1707  {
1708  fluxJumps[j] = fluxJumps[j] * jacEdge[j];
1709  }
1710  }
1711  // Non-deformed elements
1712  else
1713  {
1714  // Multiply the fluxJumps by the edge 'e' Jacobians
1715  // to bring the fluxJumps into the standard space
1716  Vmath::Smul(nEdgePts, jac[0], fluxJumps, 1,
1717  fluxJumps, 1);
1718  }
1719 
1720  // Multiply jumps by derivatives of the correction functions
1721  // All the quntities at this level should be defined into
1722  // the standard space
1723  switch (e)
1724  {
1725  case 0:
1726  for (i = 0; i < nquad0; ++i)
1727  {
1728  // Multiply fluxJumps by Q factors
1729  //fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1730 
1731  for (j = 0; j < nquad1; ++j)
1732  {
1733  cnt = i + j*nquad0;
1734  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1735  }
1736  }
1737  break;
1738  case 1:
1739  for (i = 0; i < nquad1; ++i)
1740  {
1741  // Multiply fluxJumps by Q factors
1742  //fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1743 
1744  for (j = 0; j < nquad0; ++j)
1745  {
1746  cnt = (nquad0)*i + j;
1747  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1748  }
1749  }
1750  break;
1751  case 2:
1752  for (i = 0; i < nquad0; ++i)
1753  {
1754  // Multiply fluxJumps by Q factors
1755  //fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1756 
1757  for (j = 0; j < nquad1; ++j)
1758  {
1759  cnt = j*nquad0 + i;
1760  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1761  }
1762  }
1763  break;
1764  case 3:
1765  for (i = 0; i < nquad1; ++i)
1766  {
1767  // Multiply fluxJumps by Q factors
1768  //fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1769  for (j = 0; j < nquad0; ++j)
1770  {
1771  cnt = j + i*nquad0;
1772  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1773  }
1774  }
1775  break;
1776 
1777  default:
1778  ASSERTL0(false, "edge value (< 3) is out of range");
1779  break;
1780  }
1781  }
1782 
1783  // Summing the component of the flux for calculating the
1784  // derivatives per each direction
1785  if (direction == 0)
1786  {
1787  Vmath::Vadd(nLocalSolutionPts, &divCFluxE1[0], 1,
1788  &divCFluxE3[0], 1, &derCFlux[phys_offset], 1);
1789  }
1790  else if (direction == 1)
1791  {
1792  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1793  &divCFluxE2[0], 1, &derCFlux[phys_offset], 1);
1794  }
1795  }
1796  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1088
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

◆ v_Diffuse()

void Nektar::SolverUtils::DiffusionLFR::v_Diffuse ( const std::size_t  nConvectiveFields,
const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  inarray,
Array< OneD, Array< OneD, NekDouble > > &  outarray,
const Array< OneD, Array< OneD, NekDouble > > &  pFwd = NullNekDoubleArrayofArray,
const Array< OneD, Array< OneD, NekDouble > > &  pBwd = NullNekDoubleArrayofArray 
)
protectedvirtual

Calculate FR Diffusion for the linear problems using an LDG interface flux.

Implements Nektar::SolverUtils::Diffusion.

Definition at line 853 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::LibUtilities::eGaussGaussLegendre, m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_gmat, m_IF1, m_IF2, m_jac, m_tmp1, m_tmp2, v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumFluxforScalar(), v_NumFluxforVector(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vmul().

860  {
861  boost::ignore_unused(pFwd, pBwd);
862 
863  int i, j, n;
864  int phys_offset;
865  //Array<TwoD, const NekDouble> gmat;
866  //Array<OneD, const NekDouble> jac;
867  Array<OneD, NekDouble> auxArray1, auxArray2;
868 
869  Array<OneD, LibUtilities::BasisSharedPtr> Basis;
870  Basis = fields[0]->GetExp(0)->GetBase();
871 
872  int nElements = fields[0]->GetExpSize();
873  int nDim = fields[0]->GetCoordim(0);
874  int nSolutionPts = fields[0]->GetTotPoints();
875  int nCoeffs = fields[0]->GetNcoeffs();
876 
877  Array<OneD, Array<OneD, NekDouble> > outarrayCoeff(nConvectiveFields);
878  for (i = 0; i < nConvectiveFields; ++i)
879  {
880  outarrayCoeff[i] = Array<OneD, NekDouble>(nCoeffs);
881  }
882 
883  // Compute interface numerical fluxes for inarray in physical space
884  v_NumFluxforScalar(fields, inarray, m_IF1);
885 
886  switch(nDim)
887  {
888  // 1D problems
889  case 1:
890  {
891  for (i = 0; i < nConvectiveFields; ++i)
892  {
893  // Computing the physical first-order discountinuous
894  // derivative
895  for (n = 0; n < nElements; n++)
896  {
897  phys_offset = fields[0]->GetPhys_Offset(n);
898 
899  fields[i]->GetExp(n)->PhysDeriv(0,
900  auxArray1 = inarray[i] + phys_offset,
901  auxArray2 = m_DU1[i][0] + phys_offset);
902  }
903 
904  // Computing the standard first-order correction
905  // derivative
906  v_DerCFlux_1D(nConvectiveFields, fields, inarray[i],
907  m_IF1[i][0], m_DFC1[i][0]);
908 
909  // Back to the physical space using global operations
910  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
911  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
912  Vmath::Vdiv(nSolutionPts, &m_DFC1[i][0][0], 1,
913  &m_jac[0], 1, &m_DFC1[i][0][0], 1);
914 
915  // Computing total first order derivatives
916  Vmath::Vadd(nSolutionPts, &m_DFC1[i][0][0], 1,
917  &m_DU1[i][0][0], 1, &m_D1[i][0][0], 1);
918 
919  Vmath::Vcopy(nSolutionPts, &m_D1[i][0][0], 1,
920  &m_tmp1[i][0][0], 1);
921  }
922 
923  // Computing interface numerical fluxes for m_D1
924  // in physical space
925  v_NumFluxforVector(fields, inarray, m_tmp1, m_IF2);
926 
927  for (i = 0; i < nConvectiveFields; ++i)
928  {
929  // Computing the physical second-order discountinuous
930  // derivative
931  for (n = 0; n < nElements; n++)
932  {
933  phys_offset = fields[0]->GetPhys_Offset(n);
934 
935  fields[i]->GetExp(n)->PhysDeriv(0,
936  auxArray1 = m_D1[i][0] + phys_offset,
937  auxArray2 = m_DD1[i][0] + phys_offset);
938  }
939 
940  // Computing the standard second-order correction
941  // derivative
942  v_DerCFlux_1D(nConvectiveFields, fields, m_D1[i][0],
943  m_IF2[i], m_DFC2[i][0]);
944 
945  // Back to the physical space using global operations
946  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
947  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
948  Vmath::Vdiv(nSolutionPts, &m_DFC2[i][0][0], 1,
949  &m_jac[0], 1, &m_DFC2[i][0][0], 1);
950 
951  // Adding the total divergence to outarray (RHS)
952  Vmath::Vadd(nSolutionPts, &m_DFC2[i][0][0], 1,
953  &m_DD1[i][0][0], 1, &outarray[i][0], 1);
954 
955  // Primitive Dealiasing 1D
956  if (!(Basis[0]->Collocation()))
957  {
958  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
959  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
960  }
961  }
962  break;
963  }
964  // 2D problems
965  case 2:
966  {
967  for(i = 0; i < nConvectiveFields; ++i)
968  {
969  for (j = 0; j < nDim; ++j)
970  {
971  // Temporary vectors
972  Array<OneD, NekDouble> u1_hat(nSolutionPts, 0.0);
973  Array<OneD, NekDouble> u2_hat(nSolutionPts, 0.0);
974 
975  if (j == 0)
976  {
977  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
978  &m_gmat[0][0], 1, &u1_hat[0], 1);
979 
980  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
981  &m_jac[0], 1, &u1_hat[0], 1);
982 
983  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
984  &m_gmat[1][0], 1, &u2_hat[0], 1);
985 
986  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
987  &m_jac[0], 1, &u2_hat[0], 1);
988  }
989  else if (j == 1)
990  {
991  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
992  &m_gmat[2][0], 1, &u1_hat[0], 1);
993 
994  Vmath::Vmul(nSolutionPts, &u1_hat[0], 1,
995  &m_jac[0], 1, &u1_hat[0], 1);
996 
997  Vmath::Vmul(nSolutionPts, &inarray[i][0], 1,
998  &m_gmat[3][0], 1, &u2_hat[0], 1);
999 
1000  Vmath::Vmul(nSolutionPts, &u2_hat[0], 1,
1001  &m_jac[0], 1, &u2_hat[0], 1);
1002  }
1003 
1004  for (n = 0; n < nElements; n++)
1005  {
1006  phys_offset = fields[0]->GetPhys_Offset(n);
1007 
1008  fields[i]->GetExp(n)->StdPhysDeriv(0,
1009  auxArray1 = u1_hat + phys_offset,
1010  auxArray2 = m_tmp1[i][j] + phys_offset);
1011 
1012  fields[i]->GetExp(n)->StdPhysDeriv(1,
1013  auxArray1 = u2_hat + phys_offset,
1014  auxArray2 = m_tmp2[i][j] + phys_offset);
1015  }
1016 
1017  Vmath::Vadd(nSolutionPts, &m_tmp1[i][j][0], 1,
1018  &m_tmp2[i][j][0], 1,
1019  &m_DU1[i][j][0], 1);
1020 
1021  // Divide by the metric jacobian
1022  Vmath::Vdiv(nSolutionPts, &m_DU1[i][j][0], 1,
1023  &m_jac[0], 1, &m_DU1[i][j][0], 1);
1024 
1025  // Computing the standard first-order correction
1026  // derivatives
1027  v_DerCFlux_2D(nConvectiveFields, j, fields,
1028  inarray[i], m_IF1[i][j],
1029  m_DFC1[i][j]);
1030  }
1031 
1032  // Multiplying derivatives by B matrix to get auxiliary
1033  // variables
1034  for (j = 0; j < nSolutionPts; ++j)
1035  {
1036  // std(q1)
1037  m_BD1[i][0][j] =
1038  (m_gmat[0][j]*m_gmat[0][j] +
1039  m_gmat[2][j]*m_gmat[2][j]) *
1040  m_DFC1[i][0][j] +
1041  (m_gmat[1][j]*m_gmat[0][j] +
1042  m_gmat[3][j]*m_gmat[2][j]) *
1043  m_DFC1[i][1][j];
1044 
1045  // std(q2)
1046  m_BD1[i][1][j] =
1047  (m_gmat[1][j]*m_gmat[0][j] +
1048  m_gmat[3][j]*m_gmat[2][j]) *
1049  m_DFC1[i][0][j] +
1050  (m_gmat[1][j]*m_gmat[1][j] +
1051  m_gmat[3][j]*m_gmat[3][j]) *
1052  m_DFC1[i][1][j];
1053  }
1054 
1055  // Multiplying derivatives by A^(-1) to get back
1056  // into the physical space
1057  for (j = 0; j < nSolutionPts; j++)
1058  {
1059  // q1 = A11^(-1)*std(q1) + A12^(-1)*std(q2)
1060  m_DFC1[i][0][j] =
1061  m_gmat[3][j] * m_BD1[i][0][j] -
1062  m_gmat[2][j] * m_BD1[i][1][j];
1063 
1064  // q2 = A21^(-1)*std(q1) + A22^(-1)*std(q2)
1065  m_DFC1[i][1][j] =
1066  m_gmat[0][j] * m_BD1[i][1][j] -
1067  m_gmat[1][j] * m_BD1[i][0][j];
1068  }
1069 
1070  // Computing the physical first-order derivatives
1071  for (j = 0; j < nDim; ++j)
1072  {
1073  Vmath::Vadd(nSolutionPts, &m_DU1[i][j][0], 1,
1074  &m_DFC1[i][j][0], 1,
1075  &m_D1[i][j][0], 1);
1076  }
1077  }
1078 
1079  // Computing interface numerical fluxes for m_D1
1080  // in physical space
1081  v_NumFluxforVector(fields, inarray, m_D1, m_IF2);
1082 
1083  // Computing the standard second-order discontinuous
1084  // derivatives
1085  for (i = 0; i < nConvectiveFields; ++i)
1086  {
1087  Array<OneD, NekDouble> f_hat(nSolutionPts, 0.0);
1088  Array<OneD, NekDouble> g_hat(nSolutionPts, 0.0);
1089 
1090  for (j = 0; j < nSolutionPts; j++)
1091  {
1092  f_hat[j] = (m_D1[i][0][j] * m_gmat[0][j] +
1093  m_D1[i][1][j] * m_gmat[2][j])*m_jac[j];
1094 
1095  g_hat[j] = (m_D1[i][0][j] * m_gmat[1][j] +
1096  m_D1[i][1][j] * m_gmat[3][j])*m_jac[j];
1097  }
1098 
1099  for (n = 0; n < nElements; n++)
1100  {
1101  phys_offset = fields[0]->GetPhys_Offset(n);
1102 
1103  fields[0]->GetExp(n)->StdPhysDeriv(0,
1104  auxArray1 = f_hat + phys_offset,
1105  auxArray2 = m_DD1[i][0] + phys_offset);
1106 
1107  fields[0]->GetExp(n)->StdPhysDeriv(1,
1108  auxArray1 = g_hat + phys_offset,
1109  auxArray2 = m_DD1[i][1] + phys_offset);
1110  }
1111 
1112  // Divergence of the standard discontinuous flux
1113  Vmath::Vadd(nSolutionPts,
1114  &m_DD1[i][0][0], 1,
1115  &m_DD1[i][1][0], 1,
1116  &m_divFD[i][0], 1);
1117 
1118  // Divergence of the standard correction flux
1119  if (Basis[0]->GetPointsType() ==
1120  LibUtilities::eGaussGaussLegendre &&
1121  Basis[1]->GetPointsType() ==
1122  LibUtilities::eGaussGaussLegendre)
1123  {
1124  v_DivCFlux_2D_Gauss(nConvectiveFields, fields,
1125  f_hat, g_hat,
1126  m_IF2[i], m_divFC[i]);
1127  }
1128  else
1129  {
1130  v_DivCFlux_2D(nConvectiveFields, fields,
1131  m_D1[i][0], m_D1[i][1],
1132  m_IF2[i], m_divFC[i]);
1133  }
1134 
1135  // Divergence of the standard final flux
1136  Vmath::Vadd(nSolutionPts, &m_divFD[i][0], 1,
1137  &m_divFC[i][0], 1, &outarray[i][0], 1);
1138 
1139  // Dividing by the jacobian to get back into
1140  // physical space
1141  Vmath::Vdiv(nSolutionPts, &outarray[i][0], 1,
1142  &m_jac[0], 1, &outarray[i][0], 1);
1143 
1144  // Primitive Dealiasing 2D
1145  if (!(Basis[0]->Collocation()))
1146  {
1147  fields[i]->FwdTrans(outarray[i], outarrayCoeff[i]);
1148  fields[i]->BwdTrans(outarrayCoeff[i], outarray[i]);
1149  }
1150  }//Close loop on nConvectiveFields
1151  break;
1152  }
1153  // 3D problems
1154  case 3:
1155  {
1156  ASSERTL0(false, "3D FRDG case not implemented yet");
1157  break;
1158  }
1159  }
1160  }
Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D
virtual void v_DivCFlux_2D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems.
Definition: DiffusionLFR.cpp:1811
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforVector
virtual void v_NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
Build the numerical flux for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1307
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:78
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:244
Nektar::LibUtilities::eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:54
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:86
Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar
virtual void v_NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
Builds the numerical flux for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1166
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:79
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_1D
virtual void v_DerCFlux_1D(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, const NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux for 1D problems.
Definition: DiffusionLFR.cpp:1474
Nektar::SolverUtils::DiffusionLFR::v_DerCFlux_2D
virtual void v_DerCFlux_2D(const int nConvectiveFields, const int direction, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &flux, const Array< OneD, NekDouble > &iFlux, Array< OneD, NekDouble > &derCFlux)
Compute the derivative of the corrective flux wrt a given coordinate for 2D problems.
Definition: DiffusionLFR.cpp:1594
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:84
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:88
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:77
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:81
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:80
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186
Nektar::SolverUtils::DiffusionLFR::v_DivCFlux_2D_Gauss
virtual void v_DivCFlux_2D_Gauss(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const NekDouble > &fluxX1, const Array< OneD, const NekDouble > &fluxX2, const Array< OneD, const NekDouble > &numericalFlux, Array< OneD, NekDouble > &divCFlux)
Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".
Definition: DiffusionLFR.cpp:2009
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:83

◆ v_DivCFlux_2D()

void Nektar::SolverUtils::DiffusionLFR::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 
)
protectedvirtual

Compute the divergence of the corrective flux for 2D problems.

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1 - volumetric flux in physical space.
fluxX2X2 - volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux for 2D problems.
Todo:
: Switch on shapes eventually here.

Definition at line 1811 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, m_traceNormals, Vmath::Reverse(), Vmath::Vadd(), and Vmath::Vsub().

Referenced by v_Diffuse().

1818  {
1819  boost::ignore_unused(nConvectiveFields);
1820 
1821  int n, e, i, j, cnt;
1822 
1823  int nElements = fields[0]->GetExpSize();
1824 
1825  int nLocalSolutionPts;
1826  int nEdgePts;
1827  int trace_offset;
1828  int phys_offset;
1829  int nquad0;
1830  int nquad1;
1831 
1832  Array<OneD, NekDouble> auxArray1, auxArray2;
1833  Array<OneD, LibUtilities::BasisSharedPtr> base;
1834 
1835  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1836  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
1837 
1838  // Loop on the elements
1839  for(n = 0; n < nElements; ++n)
1840  {
1841  // Offset of the element on the global vector
1842  phys_offset = fields[0]->GetPhys_Offset(n);
1843  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
1844 
1845  base = fields[0]->GetExp(n)->GetBase();
1846  nquad0 = base[0]->GetNumPoints();
1847  nquad1 = base[1]->GetNumPoints();
1848 
1849  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
1850  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
1851  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
1852  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
1853 
1854  // Loop on the edges
1855  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
1856  {
1857  // Number of edge points of edge e
1858  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
1859 
1860  Array<OneD, NekDouble> tmparrayX1(nEdgePts, 0.0);
1861  Array<OneD, NekDouble> tmparrayX2(nEdgePts, 0.0);
1862  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
1863  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
1864  Array<OneD, NekDouble> fluxJumps(nEdgePts, 0.0);
1865 
1866  // Offset of the trace space correspondent to edge e
1867  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
1868  elmtToTrace[n][e]->GetElmtId());
1869 
1870  // Get the normals of edge e
1871  const Array<OneD, const Array<OneD, NekDouble> > &normals =
1872  fields[0]->GetExp(n)->GetEdgeNormal(e);
1873 
1874  // Extract the edge values of flux-x on edge e and order
1875  // them accordingly to the order of the trace space
1876  fields[0]->GetExp(n)->GetEdgePhysVals(
1877  e, elmtToTrace[n][e],
1878  fluxX1 + phys_offset,
1879  auxArray1 = tmparrayX1);
1880 
1881  // Extract the edge values of flux-y on edge e and order
1882  // them accordingly to the order of the trace space
1883  fields[0]->GetExp(n)->GetEdgePhysVals(
1884  e, elmtToTrace[n][e],
1885  fluxX2 + phys_offset,
1886  auxArray1 = tmparrayX2);
1887 
1888  // Multiply the edge components of the flux by the normal
1889  for (i = 0; i < nEdgePts; ++i)
1890  {
1891  fluxN[i] =
1892  tmparrayX1[i]*m_traceNormals[0][trace_offset+i] +
1893  tmparrayX2[i]*m_traceNormals[1][trace_offset+i];
1894  }
1895 
1896  // Subtract to the Riemann flux the discontinuous flux
1897  Vmath::Vsub(nEdgePts,
1898  &numericalFlux[trace_offset], 1,
1899  &fluxN[0], 1, &fluxJumps[0], 1);
1900 
1901  // Check the ordering of the jump vectors
1902  if (fields[0]->GetExp(n)->GetEorient(e) ==
1903  StdRegions::eBackwards)
1904  {
1905  Vmath::Reverse(nEdgePts,
1906  auxArray1 = fluxJumps, 1,
1907  auxArray2 = fluxJumps, 1);
1908  }
1909 
1910  NekDouble fac = fields[0]->GetExp(n)->EdgeNormalNegated(e) ?
1911  -1.0 : 1.0;
1912 
1913  for (i = 0; i < nEdgePts; ++i)
1914  {
1915  if (m_traceNormals[0][trace_offset+i] != fac*normals[0][i]
1916  || m_traceNormals[1][trace_offset+i] != fac*normals[1][i])
1917  {
1918  fluxJumps[i] = -fluxJumps[i];
1919  }
1920  }
1921 
1922  // Multiply jumps by derivatives of the correction functions
1923  switch (e)
1924  {
1925  case 0:
1926  for (i = 0; i < nquad0; ++i)
1927  {
1928  // Multiply fluxJumps by Q factors
1929  fluxJumps[i] = -(m_Q2D_e0[n][i]) * fluxJumps[i];
1930 
1931  for (j = 0; j < nquad1; ++j)
1932  {
1933  cnt = i + j*nquad0;
1934  divCFluxE0[cnt] = fluxJumps[i] * m_dGL_xi2[n][j];
1935  }
1936  }
1937  break;
1938  case 1:
1939  for (i = 0; i < nquad1; ++i)
1940  {
1941  // Multiply fluxJumps by Q factors
1942  fluxJumps[i] = (m_Q2D_e1[n][i]) * fluxJumps[i];
1943 
1944  for (j = 0; j < nquad0; ++j)
1945  {
1946  cnt = (nquad0)*i + j;
1947  divCFluxE1[cnt] = fluxJumps[i] * m_dGR_xi1[n][j];
1948  }
1949  }
1950  break;
1951  case 2:
1952  for (i = 0; i < nquad0; ++i)
1953  {
1954  // Multiply fluxJumps by Q factors
1955  fluxJumps[i] = (m_Q2D_e2[n][i]) * fluxJumps[i];
1956 
1957  for (j = 0; j < nquad1; ++j)
1958  {
1959  cnt = j*nquad0 + i;
1960  divCFluxE2[cnt] = fluxJumps[i] * m_dGR_xi2[n][j];
1961  }
1962  }
1963  break;
1964  case 3:
1965  for (i = 0; i < nquad1; ++i)
1966  {
1967  // Multiply fluxJumps by Q factors
1968  fluxJumps[i] = -(m_Q2D_e3[n][i]) * fluxJumps[i];
1969  for (j = 0; j < nquad0; ++j)
1970  {
1971  cnt = j + i*nquad0;
1972  divCFluxE3[cnt] = fluxJumps[i] * m_dGL_xi1[n][j];
1973  }
1974  }
1975  break;
1976 
1977  default:
1978  ASSERTL0(false,"edge value (< 3) is out of range");
1979  break;
1980  }
1981  }
1982 
1983  // Sum each edge contribution
1984  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
1985  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
1986 
1987  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1988  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
1989 
1990  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
1991  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
1992  }
1993  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:59
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1088
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:57
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

◆ v_DivCFlux_2D_Gauss()

void Nektar::SolverUtils::DiffusionLFR::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 
)
protectedvirtual

Compute the divergence of the corrective flux for 2D problems where POINTSTYPE="GaussGaussLegendre".

Parameters
nConvectiveFieldsNumber of fields.
fieldsPointer to fields.
fluxX1X1-volumetric flux in physical space.
fluxX2X2-volumetric flux in physical space.
numericalFluxInterface flux in physical space.
divCFluxDivergence of the corrective flux.
Todo:
: Switch on shapes eventually here.

Definition at line 2009 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::StdRegions::eBackwards, m_dGL_xi1, m_dGL_xi2, m_dGR_xi1, m_dGR_xi2, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, m_traceNormals, Vmath::Neg(), Vmath::Reverse(), Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vsub().

Referenced by v_Diffuse().

2016  {
2017  boost::ignore_unused(nConvectiveFields);
2018 
2019  int n, e, i, j, cnt;
2020 
2021  int nElements = fields[0]->GetExpSize();
2022  int nLocalSolutionPts;
2023  int nEdgePts;
2024  int trace_offset;
2025  int phys_offset;
2026  int nquad0;
2027  int nquad1;
2028 
2029  Array<OneD, NekDouble> auxArray1, auxArray2;
2030  Array<OneD, LibUtilities::BasisSharedPtr> base;
2031 
2032  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
2033  &elmtToTrace = fields[0]->GetTraceMap()->GetElmtToTrace();
2034 
2035  // Loop on the elements
2036  for(n = 0; n < nElements; ++n)
2037  {
2038  // Offset of the element on the global vector
2039  phys_offset = fields[0]->GetPhys_Offset(n);
2040  nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
2041 
2042  base = fields[0]->GetExp(n)->GetBase();
2043  nquad0 = base[0]->GetNumPoints();
2044  nquad1 = base[1]->GetNumPoints();
2045 
2046  Array<OneD, NekDouble> divCFluxE0(nLocalSolutionPts, 0.0);
2047  Array<OneD, NekDouble> divCFluxE1(nLocalSolutionPts, 0.0);
2048  Array<OneD, NekDouble> divCFluxE2(nLocalSolutionPts, 0.0);
2049  Array<OneD, NekDouble> divCFluxE3(nLocalSolutionPts, 0.0);
2050 
2051  // Loop on the edges
2052  for(e = 0; e < fields[0]->GetExp(n)->GetNedges(); ++e)
2053  {
2054  // Number of edge points of edge e
2055  nEdgePts = fields[0]->GetExp(n)->GetEdgeNumPoints(e);
2056 
2057  Array<OneD, NekDouble> fluxN (nEdgePts, 0.0);
2058  Array<OneD, NekDouble> fluxT (nEdgePts, 0.0);
2059  Array<OneD, NekDouble> fluxN_R (nEdgePts, 0.0);
2060  Array<OneD, NekDouble> fluxN_D (nEdgePts, 0.0);
2061  Array<OneD, NekDouble> fluxJumps (nEdgePts, 0.0);
2062 
2063  // Offset of the trace space correspondent to edge e
2064  trace_offset = fields[0]->GetTrace()->GetPhys_Offset(
2065  elmtToTrace[n][e]->GetElmtId());
2066 
2067  // Get the normals of edge e
2068  const Array<OneD, const Array<OneD, NekDouble> > &normals =
2069  fields[0]->GetExp(n)->GetEdgeNormal(e);
2070 
2071  // Extract the trasformed normal flux at each edge
2072  switch (e)
2073  {
2074  case 0:
2075  // Extract the edge values of transformed flux-y on
2076  // edge e and order them accordingly to the order of
2077  // the trace space
2078  fields[0]->GetExp(n)->GetEdgePhysVals(
2079  e, elmtToTrace[n][e],
2080  fluxX2 + phys_offset,
2081  auxArray1 = fluxN_D);
2082 
2083  Vmath::Neg (nEdgePts, fluxN_D, 1);
2084 
2085  // Extract the physical Rieamann flux at each edge
2086  Vmath::Vcopy(nEdgePts,
2087  &numericalFlux[trace_offset], 1,
2088  &fluxN[0], 1);
2089 
2090  // Check the ordering of vectors
2091  if (fields[0]->GetExp(n)->GetEorient(e) ==
2092  StdRegions::eBackwards)
2093  {
2094  Vmath::Reverse(nEdgePts,
2095  auxArray1 = fluxN, 1,
2096  auxArray2 = fluxN, 1);
2097 
2098  Vmath::Reverse(nEdgePts,
2099  auxArray1 = fluxN_D, 1,
2100  auxArray2 = fluxN_D, 1);
2101  }
2102 
2103  // Transform Riemann Fluxes in the standard element
2104  for (i = 0; i < nquad0; ++i)
2105  {
2106  // Multiply Riemann Flux by Q factors
2107  fluxN_R[i] = (m_Q2D_e0[n][i]) * fluxN[i];
2108  }
2109 
2110  for (i = 0; i < nEdgePts; ++i)
2111  {
2112  if (m_traceNormals[0][trace_offset+i] !=
2113  normals[0][i] ||
2114  m_traceNormals[1][trace_offset+i] !=
2115  normals[1][i])
2116  {
2117  fluxN_R[i] = -fluxN_R[i];
2118  }
2119  }
2120 
2121  // Subtract to the Riemann flux the discontinuous
2122  // flux
2123  Vmath::Vsub(nEdgePts,
2124  &fluxN_R[0], 1,
2125  &fluxN_D[0], 1, &fluxJumps[0], 1);
2126 
2127  // Multiplicate the flux jumps for the correction
2128  // function
2129  for (i = 0; i < nquad0; ++i)
2130  {
2131  for (j = 0; j < nquad1; ++j)
2132  {
2133  cnt = i + j*nquad0;
2134  divCFluxE0[cnt] =
2135  -fluxJumps[i] * m_dGL_xi2[n][j];
2136  }
2137  }
2138  break;
2139  case 1:
2140  // Extract the edge values of transformed flux-x on
2141  // edge e and order them accordingly to the order of
2142  // the trace space
2143  fields[0]->GetExp(n)->GetEdgePhysVals(
2144  e, elmtToTrace[n][e],
2145  fluxX1 + phys_offset,
2146  auxArray1 = fluxN_D);
2147 
2148  // Extract the physical Rieamann flux at each edge
2149  Vmath::Vcopy(nEdgePts,
2150  &numericalFlux[trace_offset], 1,
2151  &fluxN[0], 1);
2152 
2153  // Check the ordering of vectors
2154  if (fields[0]->GetExp(n)->GetEorient(e) ==
2155  StdRegions::eBackwards)
2156  {
2157  Vmath::Reverse(nEdgePts,
2158  auxArray1 = fluxN, 1,
2159  auxArray2 = fluxN, 1);
2160 
2161  Vmath::Reverse(nEdgePts,
2162  auxArray1 = fluxN_D, 1,
2163  auxArray2 = fluxN_D, 1);
2164  }
2165 
2166  // Transform Riemann Fluxes in the standard element
2167  for (i = 0; i < nquad1; ++i)
2168  {
2169  // Multiply Riemann Flux by Q factors
2170  fluxN_R[i] = (m_Q2D_e1[n][i]) * fluxN[i];
2171  }
2172 
2173  for (i = 0; i < nEdgePts; ++i)
2174  {
2175  if (m_traceNormals[0][trace_offset+i] !=
2176  normals[0][i] ||
2177  m_traceNormals[1][trace_offset+i] !=
2178  normals[1][i])
2179  {
2180  fluxN_R[i] = -fluxN_R[i];
2181  }
2182  }
2183 
2184  // Subtract to the Riemann flux the discontinuous
2185  // flux
2186  Vmath::Vsub(nEdgePts,
2187  &fluxN_R[0], 1,
2188  &fluxN_D[0], 1, &fluxJumps[0], 1);
2189 
2190  // Multiplicate the flux jumps for the correction
2191  // function
2192  for (i = 0; i < nquad1; ++i)
2193  {
2194  for (j = 0; j < nquad0; ++j)
2195  {
2196  cnt = (nquad0)*i + j;
2197  divCFluxE1[cnt] =
2198  fluxJumps[i] * m_dGR_xi1[n][j];
2199  }
2200  }
2201  break;
2202  case 2:
2203 
2204  // Extract the edge values of transformed flux-y on
2205  // edge e and order them accordingly to the order of
2206  // the trace space
2207 
2208  fields[0]->GetExp(n)->GetEdgePhysVals(
2209  e, elmtToTrace[n][e],
2210  fluxX2 + phys_offset,
2211  auxArray1 = fluxN_D);
2212 
2213  // Extract the physical Rieamann flux at each edge
2214  Vmath::Vcopy(nEdgePts,
2215  &numericalFlux[trace_offset], 1,
2216  &fluxN[0], 1);
2217 
2218  // Check the ordering of vectors
2219  if (fields[0]->GetExp(n)->GetEorient(e) ==
2220  StdRegions::eBackwards)
2221  {
2222  Vmath::Reverse(nEdgePts,
2223  auxArray1 = fluxN, 1,
2224  auxArray2 = fluxN, 1);
2225 
2226  Vmath::Reverse(nEdgePts,
2227  auxArray1 = fluxN_D, 1,
2228  auxArray2 = fluxN_D, 1);
2229  }
2230 
2231  // Transform Riemann Fluxes in the standard element
2232  for (i = 0; i < nquad0; ++i)
2233  {
2234  // Multiply Riemann Flux by Q factors
2235  fluxN_R[i] = (m_Q2D_e2[n][i]) * fluxN[i];
2236  }
2237 
2238  for (i = 0; i < nEdgePts; ++i)
2239  {
2240  if (m_traceNormals[0][trace_offset+i] !=
2241  normals[0][i] ||
2242  m_traceNormals[1][trace_offset+i] !=
2243  normals[1][i])
2244  {
2245  fluxN_R[i] = -fluxN_R[i];
2246  }
2247  }
2248 
2249  // Subtract to the Riemann flux the discontinuous
2250  // flux
2251 
2252  Vmath::Vsub(nEdgePts,
2253  &fluxN_R[0], 1,
2254  &fluxN_D[0], 1, &fluxJumps[0], 1);
2255 
2256  // Multiplicate the flux jumps for the correction
2257  // function
2258  for (i = 0; i < nquad0; ++i)
2259  {
2260  for (j = 0; j < nquad1; ++j)
2261  {
2262  cnt = j*nquad0 + i;
2263  divCFluxE2[cnt] =
2264  fluxJumps[i] * m_dGR_xi2[n][j];
2265  }
2266  }
2267  break;
2268  case 3:
2269  // Extract the edge values of transformed flux-x on
2270  // edge e and order them accordingly to the order of
2271  // the trace space
2272 
2273  fields[0]->GetExp(n)->GetEdgePhysVals(
2274  e, elmtToTrace[n][e],
2275  fluxX1 + phys_offset,
2276  auxArray1 = fluxN_D);
2277  Vmath::Neg (nEdgePts, fluxN_D, 1);
2278 
2279  // Extract the physical Rieamann flux at each edge
2280  Vmath::Vcopy(nEdgePts,
2281  &numericalFlux[trace_offset], 1,
2282  &fluxN[0], 1);
2283 
2284  // Check the ordering of vectors
2285  if (fields[0]->GetExp(n)->GetEorient(e) ==
2286  StdRegions::eBackwards)
2287  {
2288  Vmath::Reverse(nEdgePts,
2289  auxArray1 = fluxN, 1,
2290  auxArray2 = fluxN, 1);
2291 
2292  Vmath::Reverse(nEdgePts,
2293  auxArray1 = fluxN_D, 1,
2294  auxArray2 = fluxN_D, 1);
2295  }
2296 
2297  // Transform Riemann Fluxes in the standard element
2298  for (i = 0; i < nquad1; ++i)
2299  {
2300  // Multiply Riemann Flux by Q factors
2301  fluxN_R[i] = (m_Q2D_e3[n][i]) * fluxN[i];
2302  }
2303 
2304  for (i = 0; i < nEdgePts; ++i)
2305  {
2306  if (m_traceNormals[0][trace_offset+i] !=
2307  normals[0][i] ||
2308  m_traceNormals[1][trace_offset+i] !=
2309  normals[1][i])
2310  {
2311  fluxN_R[i] = -fluxN_R[i];
2312  }
2313  }
2314 
2315  // Subtract to the Riemann flux the discontinuous
2316  // flux
2317 
2318  Vmath::Vsub(nEdgePts,
2319  &fluxN_R[0], 1,
2320  &fluxN_D[0], 1, &fluxJumps[0], 1);
2321 
2322  // Multiplicate the flux jumps for the correction
2323  // function
2324  for (i = 0; i < nquad1; ++i)
2325  {
2326  for (j = 0; j < nquad0; ++j)
2327  {
2328  cnt = j + i*nquad0;
2329  divCFluxE3[cnt] =
2330  -fluxJumps[i] * m_dGL_xi1[n][j];
2331  }
2332  }
2333  break;
2334  default:
2335  ASSERTL0(false,"edge value (< 3) is out of range");
2336  break;
2337  }
2338  }
2339 
2340 
2341  // Sum each edge contribution
2342  Vmath::Vadd(nLocalSolutionPts, &divCFluxE0[0], 1,
2343  &divCFluxE1[0], 1, &divCFlux[phys_offset], 1);
2344 
2345  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2346  &divCFluxE2[0], 1, &divCFlux[phys_offset], 1);
2347 
2348  Vmath::Vadd(nLocalSolutionPts, &divCFlux[phys_offset], 1,
2349  &divCFluxE3[0], 1, &divCFlux[phys_offset], 1);
2350  }
2351  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:59
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
Vmath::Reverse
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1088
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:58
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:57
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:346
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302

◆ v_InitObject()

void Nektar::SolverUtils::DiffusionLFR::v_InitObject ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Initiliase DiffusionLFR objects and store them before starting the time-stepping.

This routine calls the virtual functions v_SetupMetrics and v_SetupCFunctions to initialise the objects needed by DiffusionLFR.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Reimplemented from Nektar::SolverUtils::Diffusion.

Definition at line 85 of file DiffusionLFR.cpp.

References m_BD1, m_D1, m_DD1, m_DFC1, m_DFC2, m_divFC, m_divFD, m_DU1, m_IF1, m_IF2, m_session, m_tmp1, m_tmp2, v_SetupCFunctions(), and v_SetupMetrics().

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  }
Nektar::SolverUtils::DiffusionLFR::m_DU1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DU1
Definition: DiffusionLFR.h:78
Nektar::SolverUtils::DiffusionLFR::m_divFC
Array< OneD, Array< OneD, NekDouble > > m_divFC
Definition: DiffusionLFR.h:86
Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions
virtual void v_SetupCFunctions(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–7...
Definition: DiffusionLFR.cpp:326
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:79
Nektar::SolverUtils::DiffusionLFR::v_SetupMetrics
virtual void v_SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
Definition: DiffusionLFR.cpp:174
Nektar::SolverUtils::DiffusionLFR::m_DFC2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC2
Definition: DiffusionLFR.h:84
Nektar::SolverUtils::DiffusionLFR::m_tmp1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp1
Definition: DiffusionLFR.h:88
Nektar::SolverUtils::DiffusionLFR::m_tmp2
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_tmp2
Definition: DiffusionLFR.h:89
Nektar::SolverUtils::DiffusionLFR::m_IF1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_IF1
Definition: DiffusionLFR.h:77
Nektar::SolverUtils::DiffusionLFR::m_DD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DD1
Definition: DiffusionLFR.h:82
Nektar::SolverUtils::DiffusionLFR::m_divFD
Array< OneD, Array< OneD, NekDouble > > m_divFD
Definition: DiffusionLFR.h:85
Nektar::SolverUtils::DiffusionLFR::m_D1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_D1
Definition: DiffusionLFR.h:81
Nektar::SolverUtils::DiffusionLFR::m_BD1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_BD1
Definition: DiffusionLFR.h:80
Nektar::SolverUtils::DiffusionLFR::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_IF2
Array< OneD, Array< OneD, NekDouble > > m_IF2
Definition: DiffusionLFR.h:83

◆ v_NumFluxforScalar()

void Nektar::SolverUtils::DiffusionLFR::v_NumFluxforScalar ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const Array< OneD, Array< OneD, NekDouble > > &  ufield,
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &  uflux 
)
protectedvirtual

Builds the numerical flux for the 1st order derivatives.

Definition at line 1166 of file DiffusionLFR.cpp.

References m_traceNormals, Vmath::Svtvp(), v_WeakPenaltyforScalar(), and Vmath::Vcopy().

Referenced by v_Diffuse().

1170  {
1171  int i, j;
1172  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1173  int nvariables = fields.num_elements();
1174  int nDim = fields[0]->GetCoordim(0);
1175 
1176  Array<OneD, NekDouble > Fwd (nTracePts);
1177  Array<OneD, NekDouble > Bwd (nTracePts);
1178  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1179  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
1180 
1181  // Get the normal velocity Vn
1182  for (i = 0; i < nDim; ++i)
1183  {
1184  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1185  Vn, 1, Vn, 1);
1186  }
1187 
1188  // Get the sign of (v \cdot n), v = an arbitrary vector
1189  // Evaluate upwind flux:
1190  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
1191  for (j = 0; j < nDim; ++j)
1192  {
1193  for (i = 0; i < nvariables ; ++i)
1194  {
1195  // Compute Fwd and Bwd value of ufield of i direction
1196  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1197 
1198  // if Vn >= 0, flux = uFwd, i.e.,
1199  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
1200  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
1201 
1202  // else if Vn < 0, flux = uBwd, i.e.,
1203  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
1204  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
1205 
1206  fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, fluxtemp);
1207 
1208  // Imposing weak boundary condition with flux
1209  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
1210  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
1211  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
1212 
1213  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
1214  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
1215  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
1216 
1217  if(fields[0]->GetBndCondExpansions().num_elements())
1218  {
1219  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
1220  }
1221 
1222  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
1223  // i.e,
1224  // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1225  // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1226 
1227  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
1228  // i.e,
1229  // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
1230  // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
1231 
1232  Vmath::Vcopy(nTracePts, fluxtemp, 1, uflux[i][j], 1);
1233  }
1234  }
1235  }
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:488
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar
virtual void v_WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
Imposes appropriate bcs for the 1st order derivatives.
Definition: DiffusionLFR.cpp:1241
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064

◆ v_NumFluxforVector()

void Nektar::SolverUtils::DiffusionLFR::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 
)
protectedvirtual

Build the numerical flux for the 2nd order derivatives.

Definition at line 1307 of file DiffusionLFR.cpp.

References m_traceNormals, Vmath::Svtvp(), v_WeakPenaltyforVector(), Vmath::Vadd(), and Vmath::Vmul().

Referenced by v_Diffuse().

1312  {
1313  boost::ignore_unused(ufield);
1314 
1315  int i, j;
1316  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1317  int nvariables = fields.num_elements();
1318  int nDim = fields[0]->GetCoordim(0);
1319 
1320  NekDouble C11 = 0.0;
1321  Array<OneD, NekDouble > Fwd(nTracePts);
1322  Array<OneD, NekDouble > Bwd(nTracePts);
1323  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
1324 
1325  Array<OneD, NekDouble > qFwd (nTracePts);
1326  Array<OneD, NekDouble > qBwd (nTracePts);
1327  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
1328  Array<OneD, NekDouble > uterm(nTracePts);
1329 
1330  // Get the normal velocity Vn
1331  for(i = 0; i < nDim; ++i)
1332  {
1333  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
1334  Vn, 1, Vn, 1);
1335  }
1336 
1337  // Evaulate upwind flux:
1338  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
1339  for (i = 0; i < nvariables; ++i)
1340  {
1341  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
1342  for (j = 0; j < nDim; ++j)
1343  {
1344  // Compute Fwd and Bwd value of ufield of jth direction
1345  fields[i]->GetFwdBwdTracePhys(qfield[i][j], qFwd, qBwd);
1346 
1347  // if Vn >= 0, flux = uFwd, i.e.,
1348  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
1349  // qflux = qBwd = q+
1350  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
1351  // qflux = qBwd = q-
1352 
1353  // else if Vn < 0, flux = uBwd, i.e.,
1354  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
1355  // qflux = qFwd = q-
1356  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
1357  // qflux = qFwd = q+
1358 
1359  fields[i]->GetTrace()->Upwind(Vn, qBwd, qFwd, qfluxtemp);
1360 
1361  Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
1362  qfluxtemp, 1);
1363 
1364  /*
1365  // Generate Stability term = - C11 ( u- - u+ )
1366  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
1367 
1368  Vmath::Vsub(nTracePts,
1369  Fwd, 1, Bwd, 1,
1370  uterm, 1);
1371 
1372  Vmath::Smul(nTracePts,
1373  -1.0 * C11, uterm, 1,
1374  uterm, 1);
1375 
1376  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
1377  Vmath::Vadd(nTracePts,
1378  uterm, 1,
1379  qfluxtemp, 1,
1380  qfluxtemp, 1);
1381  */
1382 
1383  // Imposing weak boundary condition with flux
1384  if (fields[0]->GetBndCondExpansions().num_elements())
1385  {
1386  v_WeakPenaltyforVector(fields, i, j, qfield[i][j],
1387  qfluxtemp, C11);
1388  }
1389 
1390  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
1391  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
1392  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
1393  Vmath::Vadd(nTracePts, qfluxtemp, 1, qflux[i], 1,
1394  qflux[i], 1);
1395  }
1396  }
1397  }
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:488
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforVector
virtual void v_WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
Imposes appropriate bcs for the 2nd order derivatives.
Definition: DiffusionLFR.cpp:1405
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:302
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

◆ v_SetupCFunctions()

void Nektar::SolverUtils::DiffusionLFR::v_SetupCFunctions ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Setup the derivatives of the correction functions. For more details see J Sci Comput (2011) 47: 50–72.

This routine calls 3 different bases: #eDG_DG_Left - #eDG_DG_Left which recovers a nodal DG scheme, #eDG_SD_Left - #eDG_SD_Left which recovers the SD scheme, #eDG_HU_Left - #eDG_HU_Left which recovers the Huynh scheme. The values of the derivatives of the correction function are then stored into global variables and reused into the virtual functions #v_DivCFlux_1D, v_DivCFlux_2D, #v_DivCFlux_3D to compute the the divergence of the correction flux for 1D, 2D or 3D problems.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.

Definition at line 326 of file DiffusionLFR.cpp.

References ASSERTL0, Polylib::jacobd(), m_dGL_xi1, m_dGL_xi2, m_dGL_xi3, m_dGR_xi1, m_dGR_xi2, m_dGR_xi3, and m_diffType.

Referenced by v_InitObject().

329  {
330  boost::ignore_unused(pSession);
331 
332  int i, n;
333  NekDouble c0 = 0.0;
334  NekDouble c1 = 0.0;
335  NekDouble c2 = 0.0;
336  int nquad0, nquad1, nquad2;
337  int nmodes0, nmodes1, nmodes2;
338  Array<OneD, LibUtilities::BasisSharedPtr> base;
339 
340  int nElements = pFields[0]->GetExpSize();
341  int nDim = pFields[0]->GetCoordim(0);
342 
343  switch (nDim)
344  {
345  case 1:
346  {
347  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
348  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
349 
350  for (n = 0; n < nElements; ++n)
351  {
352  base = pFields[0]->GetExp(n)->GetBase();
353  nquad0 = base[0]->GetNumPoints();
354  nmodes0 = base[0]->GetNumModes();
355  Array<OneD, const NekDouble> z0;
356  Array<OneD, const NekDouble> w0;
357 
358  base[0]->GetZW(z0, w0);
359 
360  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
361  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
362 
363  // Auxiliary vectors to build up the auxiliary
364  // derivatives of the Legendre polynomials
365  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
366  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
367  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
368 
369  // Degree of the correction functions
370  int p0 = nmodes0 - 1;
371 
372  // Function sign to compute left correction function
373  NekDouble sign0 = pow(-1.0, p0);
374 
375  // Factorial factor to build the scheme
376  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
377  / (pow(2.0, p0)
378  * boost::math::tgamma(p0 + 1)
379  * boost::math::tgamma(p0 + 1));
380 
381  // Scalar parameter which recovers the FR schemes
382  if (m_diffType == "LFRDG")
383  {
384  c0 = 0.0;
385  }
386  else if (m_diffType == "LFRSD")
387  {
388  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
389  * (ap0 * boost::math::tgamma(p0 + 1))
390  * (ap0 * boost::math::tgamma(p0 + 1)));
391  }
392  else if (m_diffType == "LFRHU")
393  {
394  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
395  * (ap0 * boost::math::tgamma(p0 + 1))
396  * (ap0 * boost::math::tgamma(p0 + 1)));
397  }
398  else if (m_diffType == "LFRcmin")
399  {
400  c0 = -2.0 / ((2.0 * p0 + 1.0)
401  * (ap0 * boost::math::tgamma(p0 + 1))
402  * (ap0 * boost::math::tgamma(p0 + 1)));
403  }
404  else if (m_diffType == "LFRinf")
405  {
406  c0 = 10000000000000000.0;
407  }
408 
409  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
410  * (ap0 * boost::math::tgamma(p0 + 1))
411  * (ap0 * boost::math::tgamma(p0 + 1));
412 
413  NekDouble overeta0 = 1.0 / (1.0 + etap0);
414 
415  // Derivative of the Legendre polynomials
416  // dLp = derivative of the Legendre polynomial of order p
417  // dLpp = derivative of the Legendre polynomial of order p+1
418  // dLpm = derivative of the Legendre polynomial of order p-1
419  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
420  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
421  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
422 
423  // Building the DG_c_Left
424  for(i = 0; i < nquad0; ++i)
425  {
426  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
427  m_dGL_xi1[n][i] += dLpp0[i];
428  m_dGL_xi1[n][i] *= overeta0;
429  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
430  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
431  }
432 
433  // Building the DG_c_Right
434  for(i = 0; i < nquad0; ++i)
435  {
436  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
437  m_dGR_xi1[n][i] += dLpp0[i];
438  m_dGR_xi1[n][i] *= overeta0;
439  m_dGR_xi1[n][i] += dLp0[i];
440  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
441  }
442  }
443  break;
444  }
445  case 2:
446  {
447  LibUtilities::BasisSharedPtr BasisFR_Left0;
448  LibUtilities::BasisSharedPtr BasisFR_Right0;
449  LibUtilities::BasisSharedPtr BasisFR_Left1;
450  LibUtilities::BasisSharedPtr BasisFR_Right1;
451 
452  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
453  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
454  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
455  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
456 
457  for (n = 0; n < nElements; ++n)
458  {
459  base = pFields[0]->GetExp(n)->GetBase();
460  nquad0 = base[0]->GetNumPoints();
461  nquad1 = base[1]->GetNumPoints();
462  nmodes0 = base[0]->GetNumModes();
463  nmodes1 = base[1]->GetNumModes();
464 
465  Array<OneD, const NekDouble> z0;
466  Array<OneD, const NekDouble> w0;
467  Array<OneD, const NekDouble> z1;
468  Array<OneD, const NekDouble> w1;
469 
470  base[0]->GetZW(z0, w0);
471  base[1]->GetZW(z1, w1);
472 
473  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
474  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
475  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
476  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
477 
478  // Auxiliary vectors to build up the auxiliary
479  // derivatives of the Legendre polynomials
480  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
481  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
482  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
483  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
484  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
485  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
486 
487  // Degree of the correction functions
488  int p0 = nmodes0 - 1;
489  int p1 = nmodes1 - 1;
490 
491  // Function sign to compute left correction function
492  NekDouble sign0 = pow(-1.0, p0);
493  NekDouble sign1 = pow(-1.0, p1);
494 
495  // Factorial factor to build the scheme
496  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
497  / (pow(2.0, p0)
498  * boost::math::tgamma(p0 + 1)
499  * boost::math::tgamma(p0 + 1));
500 
501  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
502  / (pow(2.0, p1)
503  * boost::math::tgamma(p1 + 1)
504  * boost::math::tgamma(p1 + 1));
505 
506  // Scalar parameter which recovers the FR schemes
507  if (m_diffType == "LFRDG")
508  {
509  c0 = 0.0;
510  c1 = 0.0;
511  }
512  else if (m_diffType == "LFRSD")
513  {
514  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
515  * (ap0 * boost::math::tgamma(p0 + 1))
516  * (ap0 * boost::math::tgamma(p0 + 1)));
517 
518  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
519  * (ap1 * boost::math::tgamma(p1 + 1))
520  * (ap1 * boost::math::tgamma(p1 + 1)));
521  }
522  else if (m_diffType == "LFRHU")
523  {
524  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
525  * (ap0 * boost::math::tgamma(p0 + 1))
526  * (ap0 * boost::math::tgamma(p0 + 1)));
527 
528  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
529  * (ap1 * boost::math::tgamma(p1 + 1))
530  * (ap1 * boost::math::tgamma(p1 + 1)));
531  }
532  else if (m_diffType == "LFRcmin")
533  {
534  c0 = -2.0 / ((2.0 * p0 + 1.0)
535  * (ap0 * boost::math::tgamma(p0 + 1))
536  * (ap0 * boost::math::tgamma(p0 + 1)));
537 
538  c1 = -2.0 / ((2.0 * p1 + 1.0)
539  * (ap1 * boost::math::tgamma(p1 + 1))
540  * (ap1 * boost::math::tgamma(p1 + 1)));
541  }
542  else if (m_diffType == "LFRinf")
543  {
544  c0 = 10000000000000000.0;
545  c1 = 10000000000000000.0;
546  }
547 
548  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
549  * (ap0 * boost::math::tgamma(p0 + 1))
550  * (ap0 * boost::math::tgamma(p0 + 1));
551 
552  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
553  * (ap1 * boost::math::tgamma(p1 + 1))
554  * (ap1 * boost::math::tgamma(p1 + 1));
555 
556  NekDouble overeta0 = 1.0 / (1.0 + etap0);
557  NekDouble overeta1 = 1.0 / (1.0 + etap1);
558 
559  // Derivative of the Legendre polynomials
560  // dLp = derivative of the Legendre polynomial of order p
561  // dLpp = derivative of the Legendre polynomial of order p+1
562  // dLpm = derivative of the Legendre polynomial of order p-1
563  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
564  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
565  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
566 
567  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
568  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
569  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
570 
571  // Building the DG_c_Left
572  for(i = 0; i < nquad0; ++i)
573  {
574  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
575  m_dGL_xi1[n][i] += dLpp0[i];
576  m_dGL_xi1[n][i] *= overeta0;
577  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
578  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
579  }
580 
581  // Building the DG_c_Left
582  for(i = 0; i < nquad1; ++i)
583  {
584  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
585  m_dGL_xi2[n][i] += dLpp1[i];
586  m_dGL_xi2[n][i] *= overeta1;
587  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
588  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
589  }
590 
591  // Building the DG_c_Right
592  for(i = 0; i < nquad0; ++i)
593  {
594  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
595  m_dGR_xi1[n][i] += dLpp0[i];
596  m_dGR_xi1[n][i] *= overeta0;
597  m_dGR_xi1[n][i] += dLp0[i];
598  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
599  }
600 
601  // Building the DG_c_Right
602  for(i = 0; i < nquad1; ++i)
603  {
604  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
605  m_dGR_xi2[n][i] += dLpp1[i];
606  m_dGR_xi2[n][i] *= overeta1;
607  m_dGR_xi2[n][i] += dLp1[i];
608  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
609  }
610  }
611  break;
612  }
613  case 3:
614  {
615  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
616  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
617  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
618  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
619  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
620  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
621 
622  for (n = 0; n < nElements; ++n)
623  {
624  base = pFields[0]->GetExp(n)->GetBase();
625  nquad0 = base[0]->GetNumPoints();
626  nquad1 = base[1]->GetNumPoints();
627  nquad2 = base[2]->GetNumPoints();
628  nmodes0 = base[0]->GetNumModes();
629  nmodes1 = base[1]->GetNumModes();
630  nmodes2 = base[2]->GetNumModes();
631 
632  Array<OneD, const NekDouble> z0;
633  Array<OneD, const NekDouble> w0;
634  Array<OneD, const NekDouble> z1;
635  Array<OneD, const NekDouble> w1;
636  Array<OneD, const NekDouble> z2;
637  Array<OneD, const NekDouble> w2;
638 
639  base[0]->GetZW(z0, w0);
640  base[1]->GetZW(z1, w1);
641  base[1]->GetZW(z2, w2);
642 
643  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
644  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
645  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
646  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
647  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
648  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
649 
650  // Auxiliary vectors to build up the auxiliary
651  // derivatives of the Legendre polynomials
652  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
653  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
654  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
655  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
656  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
657  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
658  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
659  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
660  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
661 
662  // Degree of the correction functions
663  int p0 = nmodes0 - 1;
664  int p1 = nmodes1 - 1;
665  int p2 = nmodes2 - 1;
666 
667  // Function sign to compute left correction function
668  NekDouble sign0 = pow(-1.0, p0);
669  NekDouble sign1 = pow(-1.0, p1);
670  NekDouble sign2 = pow(-1.0, p2);
671 
672  // Factorial factor to build the scheme
673  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
674  / (pow(2.0, p0)
675  * boost::math::tgamma(p0 + 1)
676  * boost::math::tgamma(p0 + 1));
677 
678  // Factorial factor to build the scheme
679  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
680  / (pow(2.0, p1)
681  * boost::math::tgamma(p1 + 1)
682  * boost::math::tgamma(p1 + 1));
683 
684  // Factorial factor to build the scheme
685  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
686  / (pow(2.0, p2)
687  * boost::math::tgamma(p2 + 1)
688  * boost::math::tgamma(p2 + 1));
689 
690  // Scalar parameter which recovers the FR schemes
691  if (m_diffType == "LFRDG")
692  {
693  c0 = 0.0;
694  c1 = 0.0;
695  c2 = 0.0;
696  }
697  else if (m_diffType == "LFRSD")
698  {
699  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
700  * (ap0 * boost::math::tgamma(p0 + 1))
701  * (ap0 * boost::math::tgamma(p0 + 1)));
702 
703  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
704  * (ap1 * boost::math::tgamma(p1 + 1))
705  * (ap1 * boost::math::tgamma(p1 + 1)));
706 
707  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
708  * (ap2 * boost::math::tgamma(p2 + 1))
709  * (ap2 * boost::math::tgamma(p2 + 1)));
710  }
711  else if (m_diffType == "LFRHU")
712  {
713  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
714  * (ap0 * boost::math::tgamma(p0 + 1))
715  * (ap0 * boost::math::tgamma(p0 + 1)));
716 
717  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
718  * (ap1 * boost::math::tgamma(p1 + 1))
719  * (ap1 * boost::math::tgamma(p1 + 1)));
720 
721  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
722  * (ap2 * boost::math::tgamma(p2 + 1))
723  * (ap2 * boost::math::tgamma(p2 + 1)));
724  }
725  else if (m_diffType == "LFRcmin")
726  {
727  c0 = -2.0 / ((2.0 * p0 + 1.0)
728  * (ap0 * boost::math::tgamma(p0 + 1))
729  * (ap0 * boost::math::tgamma(p0 + 1)));
730 
731  c1 = -2.0 / ((2.0 * p1 + 1.0)
732  * (ap1 * boost::math::tgamma(p1 + 1))
733  * (ap1 * boost::math::tgamma(p1 + 1)));
734 
735  c2 = -2.0 / ((2.0 * p2 + 1.0)
736  * (ap2 * boost::math::tgamma(p2 + 1))
737  * (ap2 * boost::math::tgamma(p2 + 1)));
738  }
739  else if (m_diffType == "LFRinf")
740  {
741  c0 = 10000000000000000.0;
742  c1 = 10000000000000000.0;
743  c2 = 10000000000000000.0;
744  }
745 
746  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
747  * (ap0 * boost::math::tgamma(p0 + 1))
748  * (ap0 * boost::math::tgamma(p0 + 1));
749 
750  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
751  * (ap1 * boost::math::tgamma(p1 + 1))
752  * (ap1 * boost::math::tgamma(p1 + 1));
753 
754  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
755  * (ap2 * boost::math::tgamma(p2 + 1))
756  * (ap2 * boost::math::tgamma(p2 + 1));
757 
758  NekDouble overeta0 = 1.0 / (1.0 + etap0);
759  NekDouble overeta1 = 1.0 / (1.0 + etap1);
760  NekDouble overeta2 = 1.0 / (1.0 + etap2);
761 
762  // Derivative of the Legendre polynomials
763  // dLp = derivative of the Legendre polynomial of order p
764  // dLpp = derivative of the Legendre polynomial of order p+1
765  // dLpm = derivative of the Legendre polynomial of order p-1
766  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
767  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
768  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
769 
770  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
771  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
772  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
773 
774  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
775  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
776  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
777 
778  // Building the DG_c_Left
779  for(i = 0; i < nquad0; ++i)
780  {
781  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
782  m_dGL_xi1[n][i] += dLpp0[i];
783  m_dGL_xi1[n][i] *= overeta0;
784  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
785  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
786  }
787 
788  // Building the DG_c_Left
789  for(i = 0; i < nquad1; ++i)
790  {
791  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
792  m_dGL_xi2[n][i] += dLpp1[i];
793  m_dGL_xi2[n][i] *= overeta1;
794  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
795  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
796  }
797 
798  // Building the DG_c_Left
799  for(i = 0; i < nquad2; ++i)
800  {
801  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
802  m_dGL_xi3[n][i] += dLpp2[i];
803  m_dGL_xi3[n][i] *= overeta2;
804  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
805  m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
806  }
807 
808  // Building the DG_c_Right
809  for(i = 0; i < nquad0; ++i)
810  {
811  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
812  m_dGR_xi1[n][i] += dLpp0[i];
813  m_dGR_xi1[n][i] *= overeta0;
814  m_dGR_xi1[n][i] += dLp0[i];
815  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
816  }
817 
818  // Building the DG_c_Right
819  for(i = 0; i < nquad1; ++i)
820  {
821  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
822  m_dGR_xi2[n][i] += dLpp1[i];
823  m_dGR_xi2[n][i] *= overeta1;
824  m_dGR_xi2[n][i] += dLp1[i];
825  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
826  }
827 
828  // Building the DG_c_Right
829  for(i = 0; i < nquad2; ++i)
830  {
831  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
832  m_dGR_xi3[n][i] += dLpp2[i];
833  m_dGR_xi3[n][i] *= overeta2;
834  m_dGR_xi3[n][i] += dLp2[i];
835  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
836  }
837  }
838  break;
839  }
840  default:
841  {
842  ASSERTL0(false,"Expansion dimension not recognised");
843  break;
844  }
845  }
846  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:62
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFR.h:66
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:64
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:74
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFR.h:67
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:65
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1131
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

◆ v_SetupMetrics()

void Nektar::SolverUtils::DiffusionLFR::v_SetupMetrics ( LibUtilities::SessionReaderSharedPtr  pSession,
Array< OneD, MultiRegions::ExpListSharedPtr >  pFields 
)
protectedvirtual

Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform the fluxes at the interfaces of each element from the physical space to the standard space).

This routine calls the function #GetEdgeQFactors to compute and store the metric factors following an anticlockwise conventions along the edges/faces of each element. Note: for 1D problem the transformation is not needed.

Parameters
pSessionPointer to session reader.
pFieldsPointer to fields.
Todo:
Add the metric terms for 3D Hexahedra.

Definition at line 174 of file DiffusionLFR.cpp.

References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::LocalRegions::Expansion::GetMetricInfo(), m_gmat, m_jac, m_Q2D_e0, m_Q2D_e1, m_Q2D_e2, m_Q2D_e3, and m_traceNormals.

Referenced by v_InitObject().

177  {
178  boost::ignore_unused(pSession);
179 
180  int i, n;
181  int nquad0, nquad1;
182  int phys_offset;
183  int nLocalSolutionPts;
184  int nElements = pFields[0]->GetExpSize();
185  int nDimensions = pFields[0]->GetCoordim(0);
186  int nSolutionPts = pFields[0]->GetTotPoints();
187  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
188 
189  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
190  for (i = 0; i < nDimensions; ++i)
191  {
192  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
193  }
194  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
195 
196  Array<TwoD, const NekDouble> gmat;
197  Array<OneD, const NekDouble> jac;
198 
199  m_jac = Array<OneD, NekDouble>(nSolutionPts);
200 
201  Array<OneD, NekDouble> auxArray1;
202  Array<OneD, LibUtilities::BasisSharedPtr> base;
203  LibUtilities::PointsKeyVector ptsKeys;
204 
205  switch (nDimensions)
206  {
207  case 1:
208  {
209  for (n = 0; n < nElements; ++n)
210  {
211  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
212  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
213  phys_offset = pFields[0]->GetPhys_Offset(n);
214  jac = pFields[0]->GetExp(n)
215  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
216  ->GetMetricInfo()->GetJac(ptsKeys);
217  for (i = 0; i < nLocalSolutionPts; ++i)
218  {
219  m_jac[i+phys_offset] = jac[0];
220  }
221  }
222  break;
223  }
224  case 2:
225  {
226  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
227  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
228  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
229  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
230  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
231 
232  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
233  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
234  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
235  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
236 
237  for (n = 0; n < nElements; ++n)
238  {
239  base = pFields[0]->GetExp(n)->GetBase();
240  nquad0 = base[0]->GetNumPoints();
241  nquad1 = base[1]->GetNumPoints();
242 
243  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
244  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
245  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
246  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
247 
248  // Extract the Q factors at each edge point
249  pFields[0]->GetExp(n)->GetEdgeQFactors(
250  0, auxArray1 = m_Q2D_e0[n]);
251  pFields[0]->GetExp(n)->GetEdgeQFactors(
252  1, auxArray1 = m_Q2D_e1[n]);
253  pFields[0]->GetExp(n)->GetEdgeQFactors(
254  2, auxArray1 = m_Q2D_e2[n]);
255  pFields[0]->GetExp(n)->GetEdgeQFactors(
256  3, auxArray1 = m_Q2D_e3[n]);
257 
258  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
259  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
260  phys_offset = pFields[0]->GetPhys_Offset(n);
261 
262  jac = pFields[0]->GetExp(n)
263  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
264  ->GetMetricInfo()->GetJac(ptsKeys);
265  gmat = pFields[0]->GetExp(n)
266  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
267  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
268 
269  if (pFields[0]->GetExp(n)
270  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
271  ->GetMetricInfo()->GetGtype()
272  == SpatialDomains::eDeformed)
273  {
274  for (i = 0; i < nLocalSolutionPts; ++i)
275  {
276  m_jac[i+phys_offset] = jac[i];
277  m_gmat[0][i+phys_offset] = gmat[0][i];
278  m_gmat[1][i+phys_offset] = gmat[1][i];
279  m_gmat[2][i+phys_offset] = gmat[2][i];
280  m_gmat[3][i+phys_offset] = gmat[3][i];
281  }
282  }
283  else
284  {
285  for (i = 0; i < nLocalSolutionPts; ++i)
286  {
287  m_jac[i+phys_offset] = jac[0];
288  m_gmat[0][i+phys_offset] = gmat[0][0];
289  m_gmat[1][i+phys_offset] = gmat[1][0];
290  m_gmat[2][i+phys_offset] = gmat[2][0];
291  m_gmat[3][i+phys_offset] = gmat[3][0];
292  }
293  }
294  }
295  break;
296  }
297  case 3:
298  {
299  ASSERTL0(false,"3DFR Metric terms not implemented yet");
300  break;
301  }
302  default:
303  {
304  ASSERTL0(false, "Expansion dimension not recognised");
305  break;
306  }
307  }
308  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:59
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:60
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:54
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:58
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:57
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56

◆ v_WeakPenaltyforScalar()

void Nektar::SolverUtils::DiffusionLFR::v_WeakPenaltyforScalar ( const Array< OneD, MultiRegions::ExpListSharedPtr > &  fields,
const int  var,
const Array< OneD, const NekDouble > &  ufield,
Array< OneD, NekDouble > &  penaltyflux 
)
protectedvirtual

Imposes appropriate bcs for the 1st order derivatives.

Definition at line 1241 of file DiffusionLFR.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, and Vmath::Vcopy().

Referenced by v_NumFluxforScalar().

1246  {
1247  // Variable initialisation
1248  int i, e, id1, id2;
1249  int nBndEdgePts, nBndEdges;
1250  int cnt = 0;
1251  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1252  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1253 
1254  // Extract physical values of the fields at the trace space
1255  Array<OneD, NekDouble > uplus(nTracePts);
1256  fields[var]->ExtractTracePhys(ufield, uplus);
1257 
1258  // Impose boundary conditions
1259  for (i = 0; i < nBndRegions; ++i)
1260  {
1261  // Number of boundary edges related to bcRegion 'i'
1262  nBndEdges = fields[var]->
1263  GetBndCondExpansions()[i]->GetExpSize();
1264 
1265  // Weakly impose boundary conditions by modifying flux values
1266  for (e = 0; e < nBndEdges ; ++e)
1267  {
1268  // Number of points on the expansion
1269  nBndEdgePts = fields[var]->
1270  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1271 
1272  // Offset of the boundary expansion
1273  id1 = fields[var]->
1274  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1275 
1276  // Offset of the trace space related to boundary expansion
1277  id2 = fields[0]->GetTrace()->
1278  GetPhys_Offset(fields[0]->GetTraceMap()->
1279  GetBndCondTraceToGlobalTraceMap(cnt++));
1280 
1281  // Dirichlet bcs ==> uflux = gD
1282  if (fields[var]->GetBndConditions()[i]->
1283  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1284  {
1285  Vmath::Vcopy(nBndEdgePts,
1286  &(fields[var]->
1287  GetBndCondExpansions()[i]->
1288  GetPhys())[id1], 1,
1289  &penaltyflux[id2], 1);
1290  }
1291  // Neumann bcs ==> uflux = u+
1292  else if ((fields[var]->GetBndConditions()[i])->
1293  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1294  {
1295  Vmath::Vcopy(nBndEdgePts,
1296  &uplus[id2], 1,
1297  &penaltyflux[id2], 1);
1298  }
1299  }
1300  }
1301  }
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064

◆ v_WeakPenaltyforVector()

void Nektar::SolverUtils::DiffusionLFR::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 
)
protectedvirtual

Imposes appropriate bcs for the 2nd order derivatives.

Definition at line 1405 of file DiffusionLFR.cpp.

References Nektar::SpatialDomains::eDirichlet, Nektar::SpatialDomains::eNeumann, m_traceNormals, and Vmath::Vmul().

Referenced by v_NumFluxforVector().

1412  {
1413  boost::ignore_unused(C11);
1414 
1415  int i, e, id1, id2;
1416  int nBndEdges, nBndEdgePts;
1417  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1418  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1419 
1420  Array<OneD, NekDouble > uterm(nTracePts);
1421  Array<OneD, NekDouble > qtemp(nTracePts);
1422  int cnt = 0;
1423 
1424  fields[var]->ExtractTracePhys(qfield, qtemp);
1425 
1426  for (i = 0; i < nBndRegions; ++i)
1427  {
1428  nBndEdges = fields[var]->
1429  GetBndCondExpansions()[i]->GetExpSize();
1430 
1431  // Weakly impose boundary conditions by modifying flux values
1432  for (e = 0; e < nBndEdges ; ++e)
1433  {
1434  nBndEdgePts = fields[var]->
1435  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1436 
1437  id1 = fields[var]->
1438  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1439 
1440  id2 = fields[0]->GetTrace()->
1441  GetPhys_Offset(fields[0]->GetTraceMap()->
1442  GetBndCondTraceToGlobalTraceMap(cnt++));
1443 
1444  // For Dirichlet boundary condition:
1445  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1446  if (fields[var]->GetBndConditions()[i]->
1447  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1448  {
1449  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1450  &qtemp[id2], 1, &penaltyflux[id2], 1);
1451  }
1452  // For Neumann boundary condition: qflux = g_N
1453  else if ((fields[var]->GetBndConditions()[i])->
1454  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1455  {
1456  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1457  &(fields[var]->GetBndCondExpansions()[i]->
1458  GetPhys())[id1], 1, &penaltyflux[id2], 1);
1459  }
1460  }
1461  }
1462  }
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:74
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:186

Member Data Documentation

◆ m_BD1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_BD1
protected

Definition at line 80 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_D1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_D1
protected

Definition at line 81 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DD1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DD1
protected

Definition at line 82 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DFC1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DFC1
protected

Definition at line 79 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DFC2

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DFC2
protected

Definition at line 84 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_dGL_xi1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1

Definition at line 62 of file DiffusionLFR.h.

Referenced by v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGL_xi2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2

Definition at line 64 of file DiffusionLFR.h.

Referenced by v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGL_xi3

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3

Definition at line 66 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

◆ m_dGR_xi1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1

Definition at line 63 of file DiffusionLFR.h.

Referenced by v_DerCFlux_1D(), v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGR_xi2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2

Definition at line 65 of file DiffusionLFR.h.

Referenced by v_DerCFlux_2D(), v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupCFunctions().

◆ m_dGR_xi3

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3

Definition at line 67 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

◆ m_diffType

std::string Nektar::SolverUtils::DiffusionLFR::m_diffType
protected

Definition at line 91 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

◆ m_divFC

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_divFC
protected

Definition at line 86 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_divFD

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_divFD
protected

Definition at line 85 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_DU1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_DU1
protected

Definition at line 78 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_gmat

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_gmat

Definition at line 55 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

◆ m_IF1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_IF1
protected

Definition at line 77 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_IF2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_IF2
protected

Definition at line 83 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_Ixm

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixm

Definition at line 68 of file DiffusionLFR.h.

◆ m_Ixp

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixp

Definition at line 69 of file DiffusionLFR.h.

◆ m_jac

Array<OneD, NekDouble> Nektar::SolverUtils::DiffusionLFR::m_jac

Definition at line 54 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

◆ m_Q2D_e0

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0

Definition at line 57 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_Q2D_e1

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1

Definition at line 58 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_Q2D_e2

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2

Definition at line 59 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_Q2D_e3

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3

Definition at line 60 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), and v_SetupMetrics().

◆ m_session

LibUtilities::SessionReaderSharedPtr Nektar::SolverUtils::DiffusionLFR::m_session
protected

Definition at line 75 of file DiffusionLFR.h.

Referenced by v_InitObject().

◆ m_tmp1

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_tmp1
protected

Definition at line 88 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_tmp2

Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Nektar::SolverUtils::DiffusionLFR::m_tmp2
protected

Definition at line 89 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

◆ m_traceNormals

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_traceNormals
protected

Definition at line 74 of file DiffusionLFR.h.

Referenced by v_DivCFlux_2D(), v_DivCFlux_2D_Gauss(), v_NumFluxforScalar(), v_NumFluxforVector(), v_SetupMetrics(), and v_WeakPenaltyforVector().

◆ type

std::string Nektar::SolverUtils::DiffusionLFR::type
static
Initial value:
= {
GetDiffusionFactory().RegisterCreatorFunction(
"LFRDG", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRSD", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRHU", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcmin", DiffusionLFR::create),
GetDiffusionFactory().RegisterCreatorFunction(
"LFRcinf", DiffusionLFR::create)}

Definition at line 52 of file DiffusionLFR.h.

Generated on Wed Jan 20 2021 03:34:11 for Nektar++ by   doxygen 1.8.13