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

#include <DiffusionLFR.h>

Inheritance diagram for Nektar::SolverUtils::DiffusionLFR:
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)
 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
 
RiemannSolverSharedPtr m_riemann
 
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)
 
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 SetRiemannSolver (RiemannSolverSharedPtr riemann)
 
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 69 of file DiffusionLFR.cpp.

Referenced by create().

69  :m_diffType(diffType)
70  {
71  }

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:164
Nektar::SolverUtils::DiffusionLFR::DiffusionLFR
DiffusionLFR(std::string diffType)
DiffusionLFR uses the Flux Reconstruction (FR) approach to compute the diffusion term. The implementation is only for segments, quadrilaterals and hexahedra at the moment.
Definition: DiffusionLFR.cpp:69
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 1460 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().

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

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

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

Implements Nektar::SolverUtils::Diffusion.

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

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

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

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

86  {
87  m_session = pSession;
88  v_SetupMetrics(pSession, pFields);
89  v_SetupCFunctions(pSession, pFields);
90 
91  // Initialising arrays
92  int i, j;
93  int nConvectiveFields = pFields.num_elements();
94  int nDim = pFields[0]->GetCoordim(0);
95  int nSolutionPts = pFields[0]->GetTotPoints();
96  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
97 
98  m_IF1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
99  nConvectiveFields);
100  m_DU1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
101  nConvectiveFields);
102  m_DFC1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
103  nConvectiveFields);
104  m_BD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
105  nConvectiveFields);
106  m_D1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
107  nConvectiveFields);
108  m_DD1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
109  nConvectiveFields);
110  m_IF2 = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
111  m_DFC2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
112  nConvectiveFields);
113  m_divFD = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
114  m_divFC = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
115 
116  m_tmp1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
117  nConvectiveFields);
118  m_tmp2 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (
119  nConvectiveFields);
120 
121 
122 
123  for (i = 0; i < nConvectiveFields; ++i)
124  {
125  m_IF1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
126  m_DU1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
127  m_DFC1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
128  m_BD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
129  m_D1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
130  m_DD1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
131  m_IF2[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
132  m_DFC2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
133  m_divFD[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
134  m_divFC[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
135 
136  m_tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
137  m_tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
138 
139  for (j = 0; j < nDim; ++j)
140  {
141  m_IF1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
142  m_DU1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
143  m_DFC1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
144  m_BD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
145  m_D1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
146  m_DD1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
147  m_DFC2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
148 
149  m_tmp1[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
150  m_tmp2[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
151  }
152  }
153 
154  }
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:322
Nektar::SolverUtils::DiffusionLFR::m_DFC1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_DFC1
Definition: DiffusionLFR.h:80
Nektar::SolverUtils::DiffusionLFR::v_SetupMetrics
virtual void v_SetupMetrics(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Setup the metric terms to compute the contravariant fluxes. (i.e. this special metric terms transform...
Definition: DiffusionLFR.cpp:172
Nektar::SolverUtils::DiffusionLFR::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 1156 of file DiffusionLFR.cpp.

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

Referenced by v_Diffuse().

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

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

Referenced by v_Diffuse().

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

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

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

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

Referenced by v_NumFluxforScalar().

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

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

Referenced by v_NumFluxforVector().

1400  {
1401  int i, e, id1, id2;
1402  int nBndEdges, nBndEdgePts;
1403  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
1404  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
1405 
1406  Array<OneD, NekDouble > uterm(nTracePts);
1407  Array<OneD, NekDouble > qtemp(nTracePts);
1408  int cnt = 0;
1409 
1410  fields[var]->ExtractTracePhys(qfield, qtemp);
1411 
1412  for (i = 0; i < nBndRegions; ++i)
1413  {
1414  nBndEdges = fields[var]->
1415  GetBndCondExpansions()[i]->GetExpSize();
1416 
1417  // Weakly impose boundary conditions by modifying flux values
1418  for (e = 0; e < nBndEdges ; ++e)
1419  {
1420  nBndEdgePts = fields[var]->
1421  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
1422 
1423  id1 = fields[var]->
1424  GetBndCondExpansions()[i]->GetPhys_Offset(e);
1425 
1426  id2 = fields[0]->GetTrace()->
1427  GetPhys_Offset(fields[0]->GetTraceMap()->
1428  GetBndCondTraceToGlobalTraceMap(cnt++));
1429 
1430  // For Dirichlet boundary condition:
1431  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
1432  if (fields[var]->GetBndConditions()[i]->
1433  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
1434  {
1435  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1436  &qtemp[id2], 1, &penaltyflux[id2], 1);
1437  }
1438  // For Neumann boundary condition: qflux = g_N
1439  else if ((fields[var]->GetBndConditions()[i])->
1440  GetBoundaryConditionType() == SpatialDomains::eNeumann)
1441  {
1442  Vmath::Vmul(nBndEdgePts, &m_traceNormals[dir][id2], 1,
1443  &(fields[var]->GetBndCondExpansions()[i]->
1444  GetPhys())[id1], 1, &penaltyflux[id2], 1);
1445  }
1446  }
1447  }
1448  }
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:169

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 Tue Dec 8 2015 23:45:39 for Nektar++ by   doxygen 1.8.9.1