Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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:
Inheritance graph
[legend]
Collaboration diagram for Nektar::SolverUtils::DiffusionLFR:
Collaboration graph
[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 int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 Calculate FR Diffusion for the linear problems using an LDG interface flux. 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
 
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::Diffusion
SOLVER_UTILS_EXPORT void InitObject (LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
 
SOLVER_UTILS_EXPORT void Diffuse (const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
 
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 SetFluxVectorVec (DiffusionFluxVecCB fluxVector)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetFluxVectorNS (FuncPointerT func, ObjectPointerT obj)
 
template<typename FuncPointerT , typename ObjectPointerT >
void SetArtificialDiffusionVector (FuncPointerT func, ObjectPointerT obj)
 
void SetFluxVectorNS (DiffusionFluxVecCBNS fluxVector)
 
void SetHomoDerivs (Array< OneD, Array< OneD, NekDouble > > &deriv)
 
virtual Array< OneD, Array
< OneD, Array< OneD, NekDouble > > > & 		GetFluxTensor ()
 

Detailed Description

Definition at line 45 of file DiffusionLFR.h.

Constructor & Destructor Documentation

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

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

Definition at line 48 of file DiffusionLFR.h.

References DiffusionLFR().

49  {
50  return DiffusionSharedPtr(new DiffusionLFR(diffType));
51  }
Nektar::SolverUtils::DiffusionSharedPtr
boost::shared_ptr< Diffusion > DiffusionSharedPtr
A shared pointer to an EquationSystem object.
Definition: Diffusion.h:162
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
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 1464 of file DiffusionLFR.cpp.

References Nektar::StdRegions::StdExpansion::GetMetricInfo(), m_dGL_xi1, m_dGR_xi1, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vcopy().

Referenced by v_Diffuse().

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

References ASSERTL0, Nektar::StdRegions::eBackwards, Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::StdExpansion::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().

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

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

Implements Nektar::SolverUtils::Diffusion.

Definition at line 849 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().

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

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

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

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

Referenced by v_Diffuse().

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

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

Referenced by v_Diffuse().

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

327  {
328  int i, n;
329  NekDouble c0 = 0.0;
330  NekDouble c1 = 0.0;
331  NekDouble c2 = 0.0;
332  int nquad0, nquad1, nquad2;
333  int nmodes0, nmodes1, nmodes2;
334  Array<OneD, LibUtilities::BasisSharedPtr> base;
335 
336  int nElements = pFields[0]->GetExpSize();
337  int nDim = pFields[0]->GetCoordim(0);
338 
339  switch (nDim)
340  {
341  case 1:
342  {
343  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
344  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
345 
346  for (n = 0; n < nElements; ++n)
347  {
348  base = pFields[0]->GetExp(n)->GetBase();
349  nquad0 = base[0]->GetNumPoints();
350  nmodes0 = base[0]->GetNumModes();
351  Array<OneD, const NekDouble> z0;
352  Array<OneD, const NekDouble> w0;
353 
354  base[0]->GetZW(z0, w0);
355 
356  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
357  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
358 
359  // Auxiliary vectors to build up the auxiliary
360  // derivatives of the Legendre polynomials
361  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
362  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
363  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
364 
365  // Degree of the correction functions
366  int p0 = nmodes0 - 1;
367 
368  // Function sign to compute left correction function
369  NekDouble sign0 = pow(-1.0, p0);
370 
371  // Factorial factor to build the scheme
372  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
373  / (pow(2.0, p0)
374  * boost::math::tgamma(p0 + 1)
375  * boost::math::tgamma(p0 + 1));
376 
377  // Scalar parameter which recovers the FR schemes
378  if (m_diffType == "LFRDG")
379  {
380  c0 = 0.0;
381  }
382  else if (m_diffType == "LFRSD")
383  {
384  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
385  * (ap0 * boost::math::tgamma(p0 + 1))
386  * (ap0 * boost::math::tgamma(p0 + 1)));
387  }
388  else if (m_diffType == "LFRHU")
389  {
390  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
391  * (ap0 * boost::math::tgamma(p0 + 1))
392  * (ap0 * boost::math::tgamma(p0 + 1)));
393  }
394  else if (m_diffType == "LFRcmin")
395  {
396  c0 = -2.0 / ((2.0 * p0 + 1.0)
397  * (ap0 * boost::math::tgamma(p0 + 1))
398  * (ap0 * boost::math::tgamma(p0 + 1)));
399  }
400  else if (m_diffType == "LFRinf")
401  {
402  c0 = 10000000000000000.0;
403  }
404 
405  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
406  * (ap0 * boost::math::tgamma(p0 + 1))
407  * (ap0 * boost::math::tgamma(p0 + 1));
408 
409  NekDouble overeta0 = 1.0 / (1.0 + etap0);
410 
411  // Derivative of the Legendre polynomials
412  // dLp = derivative of the Legendre polynomial of order p
413  // dLpp = derivative of the Legendre polynomial of order p+1
414  // dLpm = derivative of the Legendre polynomial of order p-1
415  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
416  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
417  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
418 
419  // Building the DG_c_Left
420  for(i = 0; i < nquad0; ++i)
421  {
422  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
423  m_dGL_xi1[n][i] += dLpp0[i];
424  m_dGL_xi1[n][i] *= overeta0;
425  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
426  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
427  }
428 
429  // Building the DG_c_Right
430  for(i = 0; i < nquad0; ++i)
431  {
432  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
433  m_dGR_xi1[n][i] += dLpp0[i];
434  m_dGR_xi1[n][i] *= overeta0;
435  m_dGR_xi1[n][i] += dLp0[i];
436  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
437  }
438  }
439  break;
440  }
441  case 2:
442  {
443  LibUtilities::BasisSharedPtr BasisFR_Left0;
444  LibUtilities::BasisSharedPtr BasisFR_Right0;
445  LibUtilities::BasisSharedPtr BasisFR_Left1;
446  LibUtilities::BasisSharedPtr BasisFR_Right1;
447 
448  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
449  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
450  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
451  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
452 
453  for (n = 0; n < nElements; ++n)
454  {
455  base = pFields[0]->GetExp(n)->GetBase();
456  nquad0 = base[0]->GetNumPoints();
457  nquad1 = base[1]->GetNumPoints();
458  nmodes0 = base[0]->GetNumModes();
459  nmodes1 = base[1]->GetNumModes();
460 
461  Array<OneD, const NekDouble> z0;
462  Array<OneD, const NekDouble> w0;
463  Array<OneD, const NekDouble> z1;
464  Array<OneD, const NekDouble> w1;
465 
466  base[0]->GetZW(z0, w0);
467  base[1]->GetZW(z1, w1);
468 
469  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
470  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
471  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
472  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
473 
474  // Auxiliary vectors to build up the auxiliary
475  // derivatives of the Legendre polynomials
476  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
477  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
478  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
479  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
480  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
481  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
482 
483  // Degree of the correction functions
484  int p0 = nmodes0 - 1;
485  int p1 = nmodes1 - 1;
486 
487  // Function sign to compute left correction function
488  NekDouble sign0 = pow(-1.0, p0);
489  NekDouble sign1 = pow(-1.0, p1);
490 
491  // Factorial factor to build the scheme
492  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
493  / (pow(2.0, p0)
494  * boost::math::tgamma(p0 + 1)
495  * boost::math::tgamma(p0 + 1));
496 
497  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
498  / (pow(2.0, p1)
499  * boost::math::tgamma(p1 + 1)
500  * boost::math::tgamma(p1 + 1));
501 
502  // Scalar parameter which recovers the FR schemes
503  if (m_diffType == "LFRDG")
504  {
505  c0 = 0.0;
506  c1 = 0.0;
507  }
508  else if (m_diffType == "LFRSD")
509  {
510  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
511  * (ap0 * boost::math::tgamma(p0 + 1))
512  * (ap0 * boost::math::tgamma(p0 + 1)));
513 
514  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
515  * (ap1 * boost::math::tgamma(p1 + 1))
516  * (ap1 * boost::math::tgamma(p1 + 1)));
517  }
518  else if (m_diffType == "LFRHU")
519  {
520  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
521  * (ap0 * boost::math::tgamma(p0 + 1))
522  * (ap0 * boost::math::tgamma(p0 + 1)));
523 
524  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
525  * (ap1 * boost::math::tgamma(p1 + 1))
526  * (ap1 * boost::math::tgamma(p1 + 1)));
527  }
528  else if (m_diffType == "LFRcmin")
529  {
530  c0 = -2.0 / ((2.0 * p0 + 1.0)
531  * (ap0 * boost::math::tgamma(p0 + 1))
532  * (ap0 * boost::math::tgamma(p0 + 1)));
533 
534  c1 = -2.0 / ((2.0 * p1 + 1.0)
535  * (ap1 * boost::math::tgamma(p1 + 1))
536  * (ap1 * boost::math::tgamma(p1 + 1)));
537  }
538  else if (m_diffType == "LFRinf")
539  {
540  c0 = 10000000000000000.0;
541  c1 = 10000000000000000.0;
542  }
543 
544  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
545  * (ap0 * boost::math::tgamma(p0 + 1))
546  * (ap0 * boost::math::tgamma(p0 + 1));
547 
548  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
549  * (ap1 * boost::math::tgamma(p1 + 1))
550  * (ap1 * boost::math::tgamma(p1 + 1));
551 
552  NekDouble overeta0 = 1.0 / (1.0 + etap0);
553  NekDouble overeta1 = 1.0 / (1.0 + etap1);
554 
555  // Derivative of the Legendre polynomials
556  // dLp = derivative of the Legendre polynomial of order p
557  // dLpp = derivative of the Legendre polynomial of order p+1
558  // dLpm = derivative of the Legendre polynomial of order p-1
559  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
560  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
561  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
562 
563  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
564  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
565  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
566 
567  // Building the DG_c_Left
568  for(i = 0; i < nquad0; ++i)
569  {
570  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
571  m_dGL_xi1[n][i] += dLpp0[i];
572  m_dGL_xi1[n][i] *= overeta0;
573  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
574  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
575  }
576 
577  // Building the DG_c_Left
578  for(i = 0; i < nquad1; ++i)
579  {
580  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
581  m_dGL_xi2[n][i] += dLpp1[i];
582  m_dGL_xi2[n][i] *= overeta1;
583  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
584  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
585  }
586 
587  // Building the DG_c_Right
588  for(i = 0; i < nquad0; ++i)
589  {
590  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
591  m_dGR_xi1[n][i] += dLpp0[i];
592  m_dGR_xi1[n][i] *= overeta0;
593  m_dGR_xi1[n][i] += dLp0[i];
594  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
595  }
596 
597  // Building the DG_c_Right
598  for(i = 0; i < nquad1; ++i)
599  {
600  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
601  m_dGR_xi2[n][i] += dLpp1[i];
602  m_dGR_xi2[n][i] *= overeta1;
603  m_dGR_xi2[n][i] += dLp1[i];
604  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
605  }
606  }
607  break;
608  }
609  case 3:
610  {
611  m_dGL_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
612  m_dGR_xi1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
613  m_dGL_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
614  m_dGR_xi2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
615  m_dGL_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
616  m_dGR_xi3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
617 
618  for (n = 0; n < nElements; ++n)
619  {
620  base = pFields[0]->GetExp(n)->GetBase();
621  nquad0 = base[0]->GetNumPoints();
622  nquad1 = base[1]->GetNumPoints();
623  nquad2 = base[2]->GetNumPoints();
624  nmodes0 = base[0]->GetNumModes();
625  nmodes1 = base[1]->GetNumModes();
626  nmodes2 = base[2]->GetNumModes();
627 
628  Array<OneD, const NekDouble> z0;
629  Array<OneD, const NekDouble> w0;
630  Array<OneD, const NekDouble> z1;
631  Array<OneD, const NekDouble> w1;
632  Array<OneD, const NekDouble> z2;
633  Array<OneD, const NekDouble> w2;
634 
635  base[0]->GetZW(z0, w0);
636  base[1]->GetZW(z1, w1);
637  base[1]->GetZW(z2, w2);
638 
639  m_dGL_xi1[n] = Array<OneD, NekDouble>(nquad0);
640  m_dGR_xi1[n] = Array<OneD, NekDouble>(nquad0);
641  m_dGL_xi2[n] = Array<OneD, NekDouble>(nquad1);
642  m_dGR_xi2[n] = Array<OneD, NekDouble>(nquad1);
643  m_dGL_xi3[n] = Array<OneD, NekDouble>(nquad2);
644  m_dGR_xi3[n] = Array<OneD, NekDouble>(nquad2);
645 
646  // Auxiliary vectors to build up the auxiliary
647  // derivatives of the Legendre polynomials
648  Array<OneD,NekDouble> dLp0 (nquad0, 0.0);
649  Array<OneD,NekDouble> dLpp0 (nquad0, 0.0);
650  Array<OneD,NekDouble> dLpm0 (nquad0, 0.0);
651  Array<OneD,NekDouble> dLp1 (nquad1, 0.0);
652  Array<OneD,NekDouble> dLpp1 (nquad1, 0.0);
653  Array<OneD,NekDouble> dLpm1 (nquad1, 0.0);
654  Array<OneD,NekDouble> dLp2 (nquad2, 0.0);
655  Array<OneD,NekDouble> dLpp2 (nquad2, 0.0);
656  Array<OneD,NekDouble> dLpm2 (nquad2, 0.0);
657 
658  // Degree of the correction functions
659  int p0 = nmodes0 - 1;
660  int p1 = nmodes1 - 1;
661  int p2 = nmodes2 - 1;
662 
663  // Function sign to compute left correction function
664  NekDouble sign0 = pow(-1.0, p0);
665  NekDouble sign1 = pow(-1.0, p1);
666  NekDouble sign2 = pow(-1.0, p2);
667 
668  // Factorial factor to build the scheme
669  NekDouble ap0 = boost::math::tgamma(2 * p0 + 1)
670  / (pow(2.0, p0)
671  * boost::math::tgamma(p0 + 1)
672  * boost::math::tgamma(p0 + 1));
673 
674  // Factorial factor to build the scheme
675  NekDouble ap1 = boost::math::tgamma(2 * p1 + 1)
676  / (pow(2.0, p1)
677  * boost::math::tgamma(p1 + 1)
678  * boost::math::tgamma(p1 + 1));
679 
680  // Factorial factor to build the scheme
681  NekDouble ap2 = boost::math::tgamma(2 * p2 + 1)
682  / (pow(2.0, p2)
683  * boost::math::tgamma(p2 + 1)
684  * boost::math::tgamma(p2 + 1));
685 
686  // Scalar parameter which recovers the FR schemes
687  if (m_diffType == "LFRDG")
688  {
689  c0 = 0.0;
690  c1 = 0.0;
691  c2 = 0.0;
692  }
693  else if (m_diffType == "LFRSD")
694  {
695  c0 = 2.0 * p0 / ((2.0 * p0 + 1.0) * (p0 + 1.0)
696  * (ap0 * boost::math::tgamma(p0 + 1))
697  * (ap0 * boost::math::tgamma(p0 + 1)));
698 
699  c1 = 2.0 * p1 / ((2.0 * p1 + 1.0) * (p1 + 1.0)
700  * (ap1 * boost::math::tgamma(p1 + 1))
701  * (ap1 * boost::math::tgamma(p1 + 1)));
702 
703  c2 = 2.0 * p2 / ((2.0 * p2 + 1.0) * (p2 + 1.0)
704  * (ap2 * boost::math::tgamma(p2 + 1))
705  * (ap2 * boost::math::tgamma(p2 + 1)));
706  }
707  else if (m_diffType == "LFRHU")
708  {
709  c0 = 2.0 * (p0 + 1.0) / ((2.0 * p0 + 1.0) * p0
710  * (ap0 * boost::math::tgamma(p0 + 1))
711  * (ap0 * boost::math::tgamma(p0 + 1)));
712 
713  c1 = 2.0 * (p1 + 1.0) / ((2.0 * p1 + 1.0) * p1
714  * (ap1 * boost::math::tgamma(p1 + 1))
715  * (ap1 * boost::math::tgamma(p1 + 1)));
716 
717  c2 = 2.0 * (p2 + 1.0) / ((2.0 * p2 + 1.0) * p2
718  * (ap2 * boost::math::tgamma(p2 + 1))
719  * (ap2 * boost::math::tgamma(p2 + 1)));
720  }
721  else if (m_diffType == "LFRcmin")
722  {
723  c0 = -2.0 / ((2.0 * p0 + 1.0)
724  * (ap0 * boost::math::tgamma(p0 + 1))
725  * (ap0 * boost::math::tgamma(p0 + 1)));
726 
727  c1 = -2.0 / ((2.0 * p1 + 1.0)
728  * (ap1 * boost::math::tgamma(p1 + 1))
729  * (ap1 * boost::math::tgamma(p1 + 1)));
730 
731  c2 = -2.0 / ((2.0 * p2 + 1.0)
732  * (ap2 * boost::math::tgamma(p2 + 1))
733  * (ap2 * boost::math::tgamma(p2 + 1)));
734  }
735  else if (m_diffType == "LFRinf")
736  {
737  c0 = 10000000000000000.0;
738  c1 = 10000000000000000.0;
739  c2 = 10000000000000000.0;
740  }
741 
742  NekDouble etap0 = 0.5 * c0 * (2.0 * p0 + 1.0)
743  * (ap0 * boost::math::tgamma(p0 + 1))
744  * (ap0 * boost::math::tgamma(p0 + 1));
745 
746  NekDouble etap1 = 0.5 * c1 * (2.0 * p1 + 1.0)
747  * (ap1 * boost::math::tgamma(p1 + 1))
748  * (ap1 * boost::math::tgamma(p1 + 1));
749 
750  NekDouble etap2 = 0.5 * c2 * (2.0 * p2 + 1.0)
751  * (ap2 * boost::math::tgamma(p2 + 1))
752  * (ap2 * boost::math::tgamma(p2 + 1));
753 
754  NekDouble overeta0 = 1.0 / (1.0 + etap0);
755  NekDouble overeta1 = 1.0 / (1.0 + etap1);
756  NekDouble overeta2 = 1.0 / (1.0 + etap2);
757 
758  // Derivative of the Legendre polynomials
759  // dLp = derivative of the Legendre polynomial of order p
760  // dLpp = derivative of the Legendre polynomial of order p+1
761  // dLpm = derivative of the Legendre polynomial of order p-1
762  Polylib::jacobd(nquad0, z0.data(), &(dLp0[0]), p0, 0.0, 0.0);
763  Polylib::jacobd(nquad0, z0.data(), &(dLpp0[0]), p0+1, 0.0, 0.0);
764  Polylib::jacobd(nquad0, z0.data(), &(dLpm0[0]), p0-1, 0.0, 0.0);
765 
766  Polylib::jacobd(nquad1, z1.data(), &(dLp1[0]), p1, 0.0, 0.0);
767  Polylib::jacobd(nquad1, z1.data(), &(dLpp1[0]), p1+1, 0.0, 0.0);
768  Polylib::jacobd(nquad1, z1.data(), &(dLpm1[0]), p1-1, 0.0, 0.0);
769 
770  Polylib::jacobd(nquad2, z2.data(), &(dLp2[0]), p2, 0.0, 0.0);
771  Polylib::jacobd(nquad2, z2.data(), &(dLpp2[0]), p2+1, 0.0, 0.0);
772  Polylib::jacobd(nquad2, z2.data(), &(dLpm2[0]), p2-1, 0.0, 0.0);
773 
774  // Building the DG_c_Left
775  for(i = 0; i < nquad0; ++i)
776  {
777  m_dGL_xi1[n][i] = etap0 * dLpm0[i];
778  m_dGL_xi1[n][i] += dLpp0[i];
779  m_dGL_xi1[n][i] *= overeta0;
780  m_dGL_xi1[n][i] = dLp0[i] - m_dGL_xi1[n][i];
781  m_dGL_xi1[n][i] = 0.5 * sign0 * m_dGL_xi1[n][i];
782  }
783 
784  // Building the DG_c_Left
785  for(i = 0; i < nquad1; ++i)
786  {
787  m_dGL_xi2[n][i] = etap1 * dLpm1[i];
788  m_dGL_xi2[n][i] += dLpp1[i];
789  m_dGL_xi2[n][i] *= overeta1;
790  m_dGL_xi2[n][i] = dLp1[i] - m_dGL_xi2[n][i];
791  m_dGL_xi2[n][i] = 0.5 * sign1 * m_dGL_xi2[n][i];
792  }
793 
794  // Building the DG_c_Left
795  for(i = 0; i < nquad2; ++i)
796  {
797  m_dGL_xi3[n][i] = etap2 * dLpm2[i];
798  m_dGL_xi3[n][i] += dLpp2[i];
799  m_dGL_xi3[n][i] *= overeta2;
800  m_dGL_xi3[n][i] = dLp2[i] - m_dGL_xi3[n][i];
801  m_dGL_xi3[n][i] = 0.5 * sign2 * m_dGL_xi3[n][i];
802  }
803 
804  // Building the DG_c_Right
805  for(i = 0; i < nquad0; ++i)
806  {
807  m_dGR_xi1[n][i] = etap0 * dLpm0[i];
808  m_dGR_xi1[n][i] += dLpp0[i];
809  m_dGR_xi1[n][i] *= overeta0;
810  m_dGR_xi1[n][i] += dLp0[i];
811  m_dGR_xi1[n][i] = 0.5 * m_dGR_xi1[n][i];
812  }
813 
814  // Building the DG_c_Right
815  for(i = 0; i < nquad1; ++i)
816  {
817  m_dGR_xi2[n][i] = etap1 * dLpm1[i];
818  m_dGR_xi2[n][i] += dLpp1[i];
819  m_dGR_xi2[n][i] *= overeta1;
820  m_dGR_xi2[n][i] += dLp1[i];
821  m_dGR_xi2[n][i] = 0.5 * m_dGR_xi2[n][i];
822  }
823 
824  // Building the DG_c_Right
825  for(i = 0; i < nquad2; ++i)
826  {
827  m_dGR_xi3[n][i] = etap2 * dLpm2[i];
828  m_dGR_xi3[n][i] += dLpp2[i];
829  m_dGR_xi3[n][i] *= overeta2;
830  m_dGR_xi3[n][i] += dLp2[i];
831  m_dGR_xi3[n][i] = 0.5 * m_dGR_xi3[n][i];
832  }
833  }
834  break;
835  }
836  default:
837  {
838  ASSERTL0(false,"Expansion dimension not recognised");
839  break;
840  }
841  }
842  }
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi1
Definition: DiffusionLFR.h:63
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi3
Definition: DiffusionLFR.h:67
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi1
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi1
Definition: DiffusionLFR.h:64
Nektar::SolverUtils::DiffusionLFR::m_dGL_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGL_xi2
Definition: DiffusionLFR.h:65
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi3
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi3
Definition: DiffusionLFR.h:68
Nektar::SolverUtils::DiffusionLFR::m_dGR_xi2
Array< OneD, Array< OneD, NekDouble > > m_dGR_xi2
Definition: DiffusionLFR.h:66
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:2151
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LibUtilities::BasisSharedPtr
boost::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:80
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::StdRegions::StdExpansion::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  int i, n;
179  int nquad0, nquad1;
180  int phys_offset;
181  int nLocalSolutionPts;
182  int nElements = pFields[0]->GetExpSize();
183  int nDimensions = pFields[0]->GetCoordim(0);
184  int nSolutionPts = pFields[0]->GetTotPoints();
185  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
186 
187  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
188  for (i = 0; i < nDimensions; ++i)
189  {
190  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
191  }
192  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
193 
194  Array<TwoD, const NekDouble> gmat;
195  Array<OneD, const NekDouble> jac;
196 
197  m_jac = Array<OneD, NekDouble>(nSolutionPts);
198 
199  Array<OneD, NekDouble> auxArray1;
200  Array<OneD, LibUtilities::BasisSharedPtr> base;
201  LibUtilities::PointsKeyVector ptsKeys;
202 
203  switch (nDimensions)
204  {
205  case 1:
206  {
207  for (n = 0; n < nElements; ++n)
208  {
209  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
210  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
211  phys_offset = pFields[0]->GetPhys_Offset(n);
212  jac = pFields[0]->GetExp(n)
213  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
214  ->GetMetricInfo()->GetJac(ptsKeys);
215  for (i = 0; i < nLocalSolutionPts; ++i)
216  {
217  m_jac[i+phys_offset] = jac[0];
218  }
219  }
220  break;
221  }
222  case 2:
223  {
224  m_gmat = Array<OneD, Array<OneD, NekDouble> >(4);
225  m_gmat[0] = Array<OneD, NekDouble>(nSolutionPts);
226  m_gmat[1] = Array<OneD, NekDouble>(nSolutionPts);
227  m_gmat[2] = Array<OneD, NekDouble>(nSolutionPts);
228  m_gmat[3] = Array<OneD, NekDouble>(nSolutionPts);
229 
230  m_Q2D_e0 = Array<OneD, Array<OneD, NekDouble> >(nElements);
231  m_Q2D_e1 = Array<OneD, Array<OneD, NekDouble> >(nElements);
232  m_Q2D_e2 = Array<OneD, Array<OneD, NekDouble> >(nElements);
233  m_Q2D_e3 = Array<OneD, Array<OneD, NekDouble> >(nElements);
234 
235  for (n = 0; n < nElements; ++n)
236  {
237  base = pFields[0]->GetExp(n)->GetBase();
238  nquad0 = base[0]->GetNumPoints();
239  nquad1 = base[1]->GetNumPoints();
240 
241  m_Q2D_e0[n] = Array<OneD, NekDouble>(nquad0);
242  m_Q2D_e1[n] = Array<OneD, NekDouble>(nquad1);
243  m_Q2D_e2[n] = Array<OneD, NekDouble>(nquad0);
244  m_Q2D_e3[n] = Array<OneD, NekDouble>(nquad1);
245 
246  // Extract the Q factors at each edge point
247  pFields[0]->GetExp(n)->GetEdgeQFactors(
248  0, auxArray1 = m_Q2D_e0[n]);
249  pFields[0]->GetExp(n)->GetEdgeQFactors(
250  1, auxArray1 = m_Q2D_e1[n]);
251  pFields[0]->GetExp(n)->GetEdgeQFactors(
252  2, auxArray1 = m_Q2D_e2[n]);
253  pFields[0]->GetExp(n)->GetEdgeQFactors(
254  3, auxArray1 = m_Q2D_e3[n]);
255 
256  ptsKeys = pFields[0]->GetExp(n)->GetPointsKeys();
257  nLocalSolutionPts = pFields[0]->GetExp(n)->GetTotPoints();
258  phys_offset = pFields[0]->GetPhys_Offset(n);
259 
260  jac = pFields[0]->GetExp(n)
261  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
262  ->GetMetricInfo()->GetJac(ptsKeys);
263  gmat = pFields[0]->GetExp(n)
264  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
265  ->GetMetricInfo()->GetDerivFactors(ptsKeys);
266 
267  if (pFields[0]->GetExp(n)
268  ->as<LocalRegions::Expansion2D>()->GetGeom2D()
269  ->GetMetricInfo()->GetGtype()
270  == SpatialDomains::eDeformed)
271  {
272  for (i = 0; i < nLocalSolutionPts; ++i)
273  {
274  m_jac[i+phys_offset] = jac[i];
275  m_gmat[0][i+phys_offset] = gmat[0][i];
276  m_gmat[1][i+phys_offset] = gmat[1][i];
277  m_gmat[2][i+phys_offset] = gmat[2][i];
278  m_gmat[3][i+phys_offset] = gmat[3][i];
279  }
280  }
281  else
282  {
283  for (i = 0; i < nLocalSolutionPts; ++i)
284  {
285  m_jac[i+phys_offset] = jac[0];
286  m_gmat[0][i+phys_offset] = gmat[0][0];
287  m_gmat[1][i+phys_offset] = gmat[1][0];
288  m_gmat[2][i+phys_offset] = gmat[2][0];
289  m_gmat[3][i+phys_offset] = gmat[3][0];
290  }
291  }
292  }
293  break;
294  }
295  case 3:
296  {
297  ASSERTL0(false,"3DFR Metric terms not implemented yet");
298  break;
299  }
300  default:
301  {
302  ASSERTL0(false, "Expansion dimension not recognised");
303  break;
304  }
305  }
306  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e2
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e2
Definition: DiffusionLFR.h:60
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e3
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e3
Definition: DiffusionLFR.h:61
Nektar::SolverUtils::DiffusionLFR::m_jac
Array< OneD, NekDouble > m_jac
Definition: DiffusionLFR.h:55
Nektar::SolverUtils::DiffusionLFR::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLFR.h:75
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e1
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e1
Definition: DiffusionLFR.h:59
Nektar::SolverUtils::DiffusionLFR::m_gmat
Array< OneD, Array< OneD, NekDouble > > m_gmat
Definition: DiffusionLFR.h:56
Nektar::SolverUtils::DiffusionLFR::m_Q2D_e0
Array< OneD, Array< OneD, NekDouble > > m_Q2D_e0
Definition: DiffusionLFR.h:58
Nektar::SpatialDomains::eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:57
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 1235 of file DiffusionLFR.cpp.

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

Referenced by v_NumFluxforScalar().

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

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

Referenced by v_NumFluxforVector().

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

Member Data Documentation

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

Definition at line 81 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 82 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 83 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 80 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 85 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::DiffusionLFR::m_dGL_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().

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

Definition at line 65 of file DiffusionLFR.h.

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

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

Definition at line 67 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

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

Definition at line 64 of file DiffusionLFR.h.

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

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

Definition at line 66 of file DiffusionLFR.h.

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

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

Definition at line 68 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

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

Definition at line 92 of file DiffusionLFR.h.

Referenced by v_SetupCFunctions().

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

Definition at line 87 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 86 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 79 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 56 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

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

Definition at line 78 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 84 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixm

Definition at line 69 of file DiffusionLFR.h.

DNekMatSharedPtr Nektar::SolverUtils::DiffusionLFR::m_Ixp

Definition at line 70 of file DiffusionLFR.h.

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

Definition at line 55 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_SetupMetrics().

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

Definition at line 58 of file DiffusionLFR.h.

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

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

Definition at line 59 of file DiffusionLFR.h.

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

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

Definition at line 60 of file DiffusionLFR.h.

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

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

Definition at line 61 of file DiffusionLFR.h.

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

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

Definition at line 76 of file DiffusionLFR.h.

Referenced by v_InitObject().

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

Definition at line 89 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 90 of file DiffusionLFR.h.

Referenced by v_Diffuse(), and v_InitObject().

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

Definition at line 75 of file DiffusionLFR.h.

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

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 53 of file DiffusionLFR.h.

Generated on Wed Mar 8 2017 12:30:02 for Nektar++ by   doxygen 1.8.8