Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Types | Private Member Functions | List of all members
Nektar::MultiRegions::PreconditionerLowEnergy Class Reference

#include <PreconditionerLowEnergy.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerLowEnergy:
[legend]

Public Member Functions

 PreconditionerLowEnergy (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~PreconditionerLowEnergy ()
 
- Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
 Preconditioner (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~Preconditioner ()
 
void DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)
 
void DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)
 
void DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 
void DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void BuildPreconditioner ()
 
void InitObject ()
 
Array< OneD, NekDouble > AssembleStaticCondGlobalDiagonals ()
 Performs global assembly of diagonal entries to global Schur complement matrix. More...
 
const DNekScalBlkMatSharedPtr & GetBlockTransformedSchurCompl () const
 
const DNekScalBlkMatSharedPtr & GetBlockCMatrix () const
 
const DNekScalBlkMatSharedPtr & GetBlockInvDMatrix () const
 
const DNekScalBlkMatSharedPtr & GetBlockSchurCompl () const
 
const DNekScalBlkMatSharedPtr & GetBlockTransformationMatrix () const
 
const DNekScalBlkMatSharedPtr & GetBlockTransposedTransformationMatrix () const
 
DNekScalMatSharedPtr TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat)
 

Static Public Member Functions

static PreconditionerSharedPtr create (const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
 Creates an instance of this class. More...
 

Static Public Attributes

static std::string className
 Name of class. More...
 

Protected Member Functions

virtual void v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)
 Transform the basis vector to low energy space. More...
 
virtual void v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 transform the solution coeffiicents from low energy back to the original coefficient space. More...
 
virtual void v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis. More...
 
virtual void v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space. More...
 

Protected Attributes

DNekBlkMatSharedPtr m_BlkMat
 
DNekBlkMatSharedPtr m_RBlk
 
DNekBlkMatSharedPtr m_InvRBlk
 
Array< OneD, NekDouble > m_variablePmask
 
bool m_signChange
 
std::vector< std::pair< int, int > > m_sameBlock
 
- Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const std::weak_ptr< GlobalLinSys > m_linsys
 
PreconditionerType m_preconType
 
DNekMatSharedPtr m_preconditioner
 
std::weak_ptr< AssemblyMap > m_locToGloMap
 
LibUtilities::CommSharedPtr m_comm
 

Private Types

typedef std::map< LibUtilities::ShapeType, DNekScalMatSharedPtr > ShapeToDNekMap
 
typedef std::map< LibUtilities::ShapeType, LocalRegions::Expansion3DSharedPtr > ShapeToExpMap
 
typedef std::map< LibUtilities::ShapeType, Array< OneD, unsigned int > > ShapeToIntArrayMap
 
typedef std::map< LibUtilities::ShapeType, Array< OneD, Array< OneD, unsigned int > > > ShapeToIntArrayArrayMap
 

Private Member Functions

void SetupBlockTransformationMatrix (void)
 
void SetUpReferenceElements (ShapeToDNekMap &maxRmat, ShapeToExpMap &maxElmt, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR)
 Loop expansion and determine different variants of the transformation matrix. More...
 
void SetUpPyrMaxRMat (int nummodesmax, LocalRegions::PyrExpSharedPtr &PyrExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
 
void ReSetTetMaxRMat (int nummodesmax, LocalRegions::TetExpSharedPtr &TetExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
 
void ReSetPrismMaxRMat (int nummodesmax, LocalRegions::PrismExpSharedPtr &PirsmExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR, bool UseTetOnly)
 
DNekMatSharedPtr ExtractLocMat (LocalRegions::Expansion3DSharedPtr &locExp, DNekScalMatSharedPtr &maxRmat, LocalRegions::Expansion3DSharedPtr &expMax, Array< OneD, unsigned int > &vertMapMaxR, Array< OneD, Array< OneD, unsigned int > > &edgeMapMaxR)
 
void CreateVariablePMask (void)
 
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom (void)
 Sets up the reference tretrahedral element needed to construct a low energy basis. More...
 
SpatialDomains::PyrGeomSharedPtr CreateRefPyrGeom (void)
 Sets up the reference prismatic element needed to construct a low energy basis mapping arrays. More...
 
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom (void)
 Sets up the reference prismatic element needed to construct a low energy basis. More...
 
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom (void)
 Sets up the reference hexahedral element needed to construct a low energy basis. More...
 
virtual void v_InitObject ()
 
virtual void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
virtual void v_BuildPreconditioner ()
 Construct the low energy preconditioner from \(\mathbf{S}_{2}\). More...
 
virtual DNekScalMatSharedPtr v_TransformedSchurCompl (int n, int offset, const std::shared_ptr< DNekScalMat > &loc_mat)
 Set up the transformed block matrix system. More...
 

Detailed Description

Definition at line 53 of file PreconditionerLowEnergy.h.

Member Typedef Documentation

◆ ShapeToDNekMap

typedef std::map<LibUtilities::ShapeType, DNekScalMatSharedPtr> Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToDNekMap
private

Definition at line 110 of file PreconditionerLowEnergy.h.

◆ ShapeToExpMap

typedef std::map<LibUtilities::ShapeType, LocalRegions::Expansion3DSharedPtr > Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToExpMap
private

Definition at line 112 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayArrayMap

typedef std::map<LibUtilities::ShapeType, Array<OneD, Array<OneD, unsigned int> > > Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToIntArrayArrayMap
private

Definition at line 117 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayMap

typedef std::map<LibUtilities::ShapeType, Array<OneD, unsigned int> > Nektar::MultiRegions::PreconditionerLowEnergy::ShapeToIntArrayMap
private

Definition at line 114 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

◆ PreconditionerLowEnergy()

Nektar::MultiRegions::PreconditionerLowEnergy::PreconditionerLowEnergy ( const std::shared_ptr< GlobalLinSys > &  plinsys,
const AssemblyMapSharedPtr &  pLocToGloMap 
) 
@class PreconditionerLowEnergy

This class implements low energy preconditioning for the conjugate

gradient matrix solver.

Definition at line 67 of file PreconditionerLowEnergy.cpp.

70  : Preconditioner(plinsys, pLocToGloMap)
71  {
72  }
Nektar::MultiRegions::Preconditioner::Preconditioner
Preconditioner(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
Definition: Preconditioner.cpp:77

◆ ~PreconditionerLowEnergy()

virtual Nektar::MultiRegions::PreconditionerLowEnergy::~PreconditionerLowEnergy ( )
inlinevirtual

Definition at line 74 of file PreconditionerLowEnergy.h.

74 {}

Member Function Documentation

◆ create()

static PreconditionerSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::create ( const std::shared_ptr< GlobalLinSys > &  plinsys,
const std::shared_ptr< AssemblyMap > &  pLocToGloMap 
)
inlinestatic

Creates an instance of this class.

Definition at line 57 of file PreconditionerLowEnergy.h.

60  {
61  PreconditionerSharedPtr p = MemoryManager<PreconditionerLowEnergy>::AllocateSharedPtr(plinsys,pLocToGloMap);
62  p->InitObject();
63  return p;
64  }
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::MultiRegions::PreconditionerSharedPtr
std::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:60

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ CreateRefHexGeom()

SpatialDomains::HexGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefHexGeom ( void  )
private

Sets up the reference hexahedral element needed to construct a low energy basis.

Definition at line 1554 of file PreconditionerLowEnergy.cpp.

1555  {
1556  ////////////////////////////////
1557  // Set up Hexahedron vertices //
1558  ////////////////////////////////
1559 
1560  const int three=3;
1561 
1562  const int nVerts = 8;
1563  const double point[][3] = {
1564  {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1565  {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1566  };
1567 
1568  // Populate the list of verts
1569  SpatialDomains::PointGeomSharedPtr verts[8];
1570  for( int i = 0; i < nVerts; ++i ) {
1571  verts[i] = MemoryManager<SpatialDomains::PointGeom>
1572  ::AllocateSharedPtr(three, i, point[i][0],
1573  point[i][1], point[i][2]);
1574  }
1575 
1576  /////////////////////////////
1577  // Set up Hexahedron Edges //
1578  /////////////////////////////
1579 
1580  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1581  const int nEdges = 12;
1582  const int vertexConnectivity[][2] = {
1583  {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1584  {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1585  };
1586 
1587  // Populate the list of edges
1588  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1589  for( int i = 0; i < nEdges; ++i ) {
1590  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1591  for( int j = 0; j < 2; ++j ) {
1592  vertsArray[j] = verts[vertexConnectivity[i][j]];
1593  }
1594  edges[i] = MemoryManager<SpatialDomains::SegGeom>::
1595  AllocateSharedPtr( i, three, vertsArray);
1596  }
1597 
1598  /////////////////////////////
1599  // Set up Hexahedron faces //
1600  /////////////////////////////
1601 
1602  const int nFaces = 6;
1603  const int edgeConnectivity[][4] = {
1604  {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1605  {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1606  };
1607 
1608  // Populate the list of faces
1609  SpatialDomains::QuadGeomSharedPtr faces[nFaces];
1610  for( int i = 0; i < nFaces; ++i ) {
1611  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1612  for( int j = 0; j < 4; ++j ) {
1613  edgeArray[j] = edges[edgeConnectivity[i][j]];
1614  }
1615  faces[i] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(i, edgeArray);
1616  }
1617 
1618  SpatialDomains::HexGeomSharedPtr geom =
1619  MemoryManager<SpatialDomains::HexGeom>::AllocateSharedPtr
1620  (0, faces);
1621 
1622  return geom;
1623  }
Nektar::SpatialDomains::QuadGeomSharedPtr
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
Nektar::SpatialDomains::SegGeomSharedPtr
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
Nektar::SpatialDomains::HexGeomSharedPtr
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
Nektar::SpatialDomains::PointGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

Referenced by SetUpReferenceElements().

◆ CreateRefPrismGeom()

SpatialDomains::PrismGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPrismGeom ( void  )
private

Sets up the reference prismatic element needed to construct a low energy basis.

Definition at line 1304 of file PreconditionerLowEnergy.cpp.

1305  {
1306  //////////////////////////
1307  // Set up Prism element //
1308  //////////////////////////
1309 
1310  const int three=3;
1311  const int nVerts = 6;
1312  const double point[][3] = {
1313  {-1,-1,0}, {1,-1,0}, {1,1,0},
1314  {-1,1,0}, {0,-1,sqrt(double(3))}, {0,1,sqrt(double(3))},
1315  };
1316 
1317  //std::shared_ptr<SpatialDomains::PointGeom> verts[6];
1318  SpatialDomains::PointGeomSharedPtr verts[6];
1319  for(int i=0; i < nVerts; ++i)
1320  {
1321  verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr
1322  ( three, i, point[i][0], point[i][1], point[i][2] );
1323  }
1324  const int nEdges = 9;
1325  const int vertexConnectivity[][2] = {
1326  {0,1}, {1,2}, {3,2}, {0,3}, {0,4},
1327  {1,4}, {2,5}, {3,5}, {4,5}
1328  };
1329 
1330  // Populate the list of edges
1331  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1332  for(int i=0; i < nEdges; ++i){
1333  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1334  for(int j=0; j<2; ++j)
1335  {
1336  vertsArray[j] = verts[vertexConnectivity[i][j]];
1337  }
1338  edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(i, three, vertsArray);
1339  }
1340 
1341  ////////////////////////
1342  // Set up Prism faces //
1343  ////////////////////////
1344 
1345  const int nFaces = 5;
1346  //quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
1347  const int quadEdgeConnectivity[][4] = { {0,1,2,3}, {1,6,8,5}, {3,7,8,4} };
1348  // QuadId ordered as 0, 1, 2, otherwise return false
1349  const int quadId[] = { 0,-1,1,-1,2 };
1350 
1351  //triangle-edge connectivity side-triface-1, side triface-3
1352  const int triEdgeConnectivity[][3] = { {0,5,4}, {2,6,7} };
1353  // TriId ordered as 0, 1, otherwise return false
1354  const int triId[] = { -1,0,-1,1,-1 };
1355 
1356  // Populate the list of faces
1357  SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1358  for(int f = 0; f < nFaces; ++f){
1359  if(f == 1 || f == 3) {
1360  int i = triId[f];
1361  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1362  for(int j = 0; j < 3; ++j){
1363  edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1364  }
1365  faces[f] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(f, edgeArray);
1366  }
1367  else {
1368  int i = quadId[f];
1369  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1370  for(int j=0; j < 4; ++j){
1371  edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1372  }
1373  faces[f] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(f, edgeArray);
1374  }
1375  }
1376 
1377  SpatialDomains::PrismGeomSharedPtr geom = MemoryManager<SpatialDomains::PrismGeom>::AllocateSharedPtr(0, faces);
1378 
1379  return geom;
1380  }
Nektar::SpatialDomains::PrismGeomSharedPtr
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:82
Nektar::SpatialDomains::Geometry2DSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and tinysimd::sqrt().

Referenced by SetUpReferenceElements().

◆ CreateRefPyrGeom()

SpatialDomains::PyrGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPyrGeom ( void  )
private

Sets up the reference prismatic element needed to construct a low energy basis mapping arrays.

Definition at line 1386 of file PreconditionerLowEnergy.cpp.

1387  {
1388  //////////////////////////
1389  // Set up Pyramid element //
1390  //////////////////////////
1391 
1392  const int nVerts = 5;
1393  const double point[][3] = {
1394  {-1,-1,0}, {1,-1,0}, {1,1,0},
1395  {-1,1,0}, {0,0,sqrt(double(2))}
1396  };
1397 
1398  //boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
1399  const int three=3;
1400  SpatialDomains::PointGeomSharedPtr verts[5];
1401  for(int i=0; i < nVerts; ++i)
1402  {
1403  verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr
1404  ( three, i, point[i][0], point[i][1], point[i][2] );
1405  }
1406  const int nEdges = 8;
1407  const int vertexConnectivity[][2] = {
1408  {0,1}, {1,2}, {2,3}, {3,0},
1409  {0,4}, {1,4}, {2,4}, {3,4}
1410  };
1411 
1412  // Populate the list of edges
1413  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1414  for(int i=0; i < nEdges; ++i)
1415  {
1416  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1417  for(int j=0; j<2; ++j)
1418  {
1419  vertsArray[j] = verts[vertexConnectivity[i][j]];
1420  }
1421  edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(i, three, vertsArray);
1422  }
1423 
1424  ////////////////////////
1425  // Set up Pyramid faces //
1426  ////////////////////////
1427 
1428  const int nFaces = 5;
1429  //quad-edge connectivity base-face0,
1430  const int quadEdgeConnectivity[][4] = {{0,1,2,3}};
1431 
1432  //triangle-edge connectivity side-triface-1, side triface-2
1433  const int triEdgeConnectivity[][3] = { {0,5,4}, {1,6,5}, {2,7,6}, {3,4,7}};
1434 
1435  // Populate the list of faces
1436  SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1437  for(int f = 0; f < nFaces; ++f)
1438  {
1439  if(f == 0)
1440  {
1441  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1442  for(int j=0; j < 4; ++j)
1443  {
1444  edgeArray[j] = edges[quadEdgeConnectivity[f][j]];
1445  }
1446 
1447  faces[f] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(f,edgeArray);
1448  }
1449  else {
1450  int i = f-1;
1451  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1452  for(int j = 0; j < 3; ++j)
1453  {
1454  edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1455  }
1456  faces[f] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(f, edgeArray);
1457  }
1458  }
1459 
1460  SpatialDomains::PyrGeomSharedPtr geom =
1461  MemoryManager<SpatialDomains::PyrGeom>::AllocateSharedPtr(0,faces);
1462 
1463  return geom;
1464  }
Nektar::SpatialDomains::PyrGeomSharedPtr
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:74

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and tinysimd::sqrt().

Referenced by SetUpReferenceElements().

◆ CreateRefTetGeom()

SpatialDomains::TetGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefTetGeom ( void  )
private

Sets up the reference tretrahedral element needed to construct a low energy basis.

Definition at line 1470 of file PreconditionerLowEnergy.cpp.

1471  {
1472  /////////////////////////////////
1473  // Set up Tetrahedron vertices //
1474  /////////////////////////////////
1475 
1476  int i,j;
1477  const int three=3;
1478  const int nVerts = 4;
1479  const double point[][3] = {
1480  {-1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1481  {1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1482  {0,2/sqrt(double(3)),-1/sqrt(double(6))},
1483  {0,0,3/sqrt(double(6))}};
1484 
1485  std::shared_ptr<SpatialDomains::PointGeom> verts[4];
1486  for(i=0; i < nVerts; ++i)
1487  {
1488  verts[i] =
1489  MemoryManager<SpatialDomains::PointGeom>::
1490  AllocateSharedPtr
1491  ( three, i, point[i][0], point[i][1], point[i][2] );
1492  }
1493 
1494  //////////////////////////////
1495  // Set up Tetrahedron Edges //
1496  //////////////////////////////
1497 
1498  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1499  const int nEdges = 6;
1500  const int vertexConnectivity[][2] = {
1501  {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1502  };
1503 
1504  // Populate the list of edges
1505  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1506  for(i=0; i < nEdges; ++i)
1507  {
1508  std::shared_ptr<SpatialDomains::PointGeom>
1509  vertsArray[2];
1510  for(j=0; j<2; ++j)
1511  {
1512  vertsArray[j] = verts[vertexConnectivity[i][j]];
1513  }
1514 
1515  edges[i] = MemoryManager<SpatialDomains::SegGeom>
1516  ::AllocateSharedPtr(i, three, vertsArray);
1517  }
1518 
1519  //////////////////////////////
1520  // Set up Tetrahedron faces //
1521  //////////////////////////////
1522 
1523  const int nFaces = 4;
1524  const int edgeConnectivity[][3] = {
1525  {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1526  };
1527 
1528  // Populate the list of faces
1529  SpatialDomains::TriGeomSharedPtr faces[nFaces];
1530  for(i=0; i < nFaces; ++i)
1531  {
1532  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1533  for(j=0; j < 3; ++j)
1534  {
1535  edgeArray[j] = edges[edgeConnectivity[i][j]];
1536  }
1537 
1538 
1539  faces[i] = MemoryManager<SpatialDomains::TriGeom>
1540  ::AllocateSharedPtr(i, edgeArray);
1541  }
1542 
1543  SpatialDomains::TetGeomSharedPtr geom =
1544  MemoryManager<SpatialDomains::TetGeom>::AllocateSharedPtr
1545  (0, faces);
1546 
1547  return geom;
1548  }
Nektar::SpatialDomains::TetGeomSharedPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
Nektar::SpatialDomains::TriGeomSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and tinysimd::sqrt().

Referenced by SetUpReferenceElements().

◆ CreateVariablePMask()

void Nektar::MultiRegions::PreconditionerLowEnergy::CreateVariablePMask ( void  )
private

Create the inverse multiplicity map.

Definition at line 1274 of file PreconditionerLowEnergy.cpp.

1275  {
1276  unsigned int nLocBnd = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1277  unsigned int i;
1278  auto asmMap = m_locToGloMap.lock();
1279 
1280  const Array< OneD, const NekDouble > &sign
1281  = asmMap->GetLocalToGlobalBndSign();
1282 
1283  m_signChange=asmMap->GetSignChange();
1284 
1285  // Construct a map of 1/multiplicity
1286  m_variablePmask = Array<OneD, NekDouble>(nLocBnd);
1287  for (i = 0; i < nLocBnd; ++i)
1288  {
1289  if(m_signChange)
1290  {
1291  m_variablePmask[i] = fabs(sign[i]);
1292  }
1293  else
1294  {
1295  m_variablePmask[i] = 1.0;
1296  }
1297  }
1298  }
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Nektar::MultiRegions::Preconditioner::m_locToGloMap
std::weak_ptr< AssemblyMap > m_locToGloMap
Definition: Preconditioner.h:131

References Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_signChange, m_variablePmask, and sign.

Referenced by v_InitObject().

◆ ExtractLocMat()

DNekMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::ExtractLocMat ( LocalRegions::Expansion3DSharedPtr &  locExp,
DNekScalMatSharedPtr &  maxRmat,
LocalRegions::Expansion3DSharedPtr &  expMax,
Array< OneD, unsigned int > &  vertMapMaxR,
Array< OneD, Array< OneD, unsigned int > > &  edgeMapMaxR 
)
private

Definition at line 2244 of file PreconditionerLowEnergy.cpp.

2250  {
2251  NekDouble val;
2252  NekDouble zero = 0.0;
2253 
2254  int nRows = locExp->NumBndryCoeffs();
2255  DNekMatSharedPtr newmat = MemoryManager<DNekMat>::
2256  AllocateSharedPtr(nRows,nRows,zero,eFULL);
2257 
2258  Array<OneD, unsigned int> vlocmap;
2259  Array<OneD, Array<OneD, unsigned int> > elocmap;
2260  Array<OneD, Array<OneD, unsigned int> > flocmap;
2261 
2262  locExp->GetInverseBoundaryMaps(vlocmap,elocmap,flocmap);
2263 
2264  // fill diagonal
2265  for(int i = 0; i < nRows; ++i)
2266  {
2267  val = 1.0;
2268  newmat->SetValue(i,i,val);
2269  }
2270 
2271  int nverts = locExp->GetNverts();
2272  int nedges = locExp->GetNedges();
2273  int nfaces = locExp->GetNtraces();
2274 
2275  // fill vertex to edge coupling
2276  for(int e = 0; e < nedges; ++e)
2277  {
2278  int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) -2;
2279 
2280  for(int v = 0; v < nverts; ++v)
2281  {
2282  for(int i = 0; i < nEdgeInteriorCoeffs; ++i)
2283  {
2284  val = (*maxRmat)(vmap[v],emap[e][i]);
2285  newmat->SetValue(vlocmap[v],elocmap[e][i],val);
2286  }
2287  }
2288  }
2289 
2290  for(int f = 0; f < nfaces; ++f)
2291  {
2292  // Get details to extrac this face from max reference matrix
2293  StdRegions::Orientation FwdOrient = StdRegions::eDir1FwdDir1_Dir2FwdDir2;
2294  int m0,m1; //Local face expansion orders.
2295 
2296  int nFaceInteriorCoeffs = locExp->GetTraceIntNcoeffs(f);
2297 
2298  locExp->GetTraceNumModes(f,m0,m1,FwdOrient);
2299 
2300  Array<OneD, unsigned int> fmapRmat = maxExp->
2301  GetTraceInverseBoundaryMap(f,FwdOrient, m0,m1);
2302 
2303  // fill in vertex to face coupling
2304  for(int v = 0; v < nverts; ++v)
2305  {
2306  for(int i = 0; i < nFaceInteriorCoeffs; ++i)
2307  {
2308  val = (*maxRmat)(vmap[v],fmapRmat[i]);
2309  newmat->SetValue(vlocmap[v],flocmap[f][i],val);
2310  }
2311  }
2312 
2313  // fill in edges to face coupling
2314  for(int e = 0; e < nedges; ++e)
2315  {
2316  int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) -2;
2317 
2318  for(int j = 0; j < nEdgeInteriorCoeffs; ++j)
2319  {
2320 
2321  for(int i = 0; i < nFaceInteriorCoeffs; ++i)
2322  {
2323  val = (*maxRmat)(emap[e][j],fmapRmat[i]);
2324  newmat->SetValue(elocmap[e][j],flocmap[f][i],val);
2325  }
2326  }
2327  }
2328  }
2329 
2330  return newmat;
2331  }
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, and Nektar::eFULL.

Referenced by SetupBlockTransformationMatrix().

◆ ReSetPrismMaxRMat()

void Nektar::MultiRegions::PreconditionerLowEnergy::ReSetPrismMaxRMat ( int  nummodesmax,
LocalRegions::PrismExpSharedPtr &  PirsmExp,
ShapeToDNekMap &  maxRmat,
ShapeToIntArrayMap &  vertMapMaxR,
ShapeToIntArrayArrayMap &  edgeMapMaxR,
ShapeToIntArrayArrayMap &  faceMapMaxR,
bool  UseTetOnly 
)
private

Definition at line 2040 of file PreconditionerLowEnergy.cpp.

2047  {
2048  int nRows = PrismExp->NumBndryCoeffs();
2049  NekDouble val;
2050  NekDouble zero = 0.0;
2051  DNekMatSharedPtr newmat = MemoryManager<DNekMat>::
2052  AllocateSharedPtr(nRows,nRows,zero,eFULL);
2053 
2054 
2055  int nfacemodes;
2056 
2057  if(UseTetOnly)
2058  {
2059  // copy existing system
2060  for(int i = 0; i < nRows; ++i)
2061  {
2062  for(int j = 0; j < nRows; ++j)
2063  {
2064  val = (*maxRmat[ePrism])(i,j);
2065  newmat->SetValue(i,j,val);
2066  }
2067  }
2068 
2069  // Reset vertex to edge mapping from tet.
2070  const int VETetVert[][2] = {{0,0},{1,1},{1,1},{0,0},{3,3},{3,3}};
2071  const int VETetEdge[][2] = {{0,3},{0,4},{0,4},{0,3},{3,4},{4,3}};
2072  const int VEPrismEdge[][2] = {{0,4},{0,5},{2,6},{2,7},{4,5},{6,7}};
2073 
2074  // fill vertex to edge coupling
2075  for(int v = 0; v < 6; ++v)
2076  {
2077  for(int e = 0; e < 2; ++e)
2078  {
2079  for(int i = 0; i < nummodesmax-2; ++i)
2080  {
2081  // Note this is using wrong shape but gives
2082  // answer that seems to have correct error!
2083  val = (*maxRmat[eTetrahedron])(
2084  vertMapMaxR[eTetrahedron][VETetVert[v][e]],
2085  edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i]);
2086  newmat->
2087  SetValue(vertMapMaxR[ePrism][v],
2088  edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
2089  val);
2090  }
2091  }
2092  }
2093  }
2094  else
2095  {
2096 
2097  // set diagonal to 1
2098  for(int i = 0; i < nRows; ++i)
2099  {
2100  newmat->SetValue(i,i,1.0);
2101  }
2102 
2103 
2104  // Set vertex to edge mapping from Hex.
2105 
2106  // The following lists specify the number of adjacent
2107  // edges to each vertex (nadj) then the Hex vert to
2108  // use for each prism ver in the vert-edge map (VEVert)
2109  // followed by the hex edge to use for each prism edge
2110  // in the vert-edge map (VEEdge)
2111  const int VEHexVert[][3] = {{0,0,0},{1,1,1},{2,2,2},{3,3,3},
2112  {4,5,5},{6,7,7}};
2113  const int VEHexEdge[][3] = {{0,3,4},{0,1,5},{1,2,6},{2,3,7},
2114  {4,5,9},{6,7,11}};
2115  const int VEPrismEdge[][3] = {{0,3,4},{0,1,5},{1,2,6},{2,3,7},
2116  {4,5,8},{6,7,8}};
2117 
2118  // fill vertex to edge coupling
2119  for(int v = 0; v < 6; ++v)
2120  {
2121  for(int e = 0; e < 3; ++e)
2122  {
2123  for(int i = 0; i < nummodesmax-2; ++i)
2124  {
2125  // Note this is using wrong shape but gives
2126  // answer that seems to have correct error!
2127  val = (*maxRmat[eHexahedron])(
2128  vertMapMaxR[eHexahedron][VEHexVert[v][e]],
2129  edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
2130  newmat->SetValue(vertMapMaxR[ePrism][v],
2131  edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
2132  val);
2133  }
2134  }
2135  }
2136 
2137 
2138  // Setup vertex to face mapping from Hex
2139  const int VFHexVert[][2] = {{0,0},{1,1},{4,5},{2,2},{3,3},{6,7}};
2140  const int VFHexFace[][2] = {{0,4},{0,2},{4,2},{0,2},{0,4},{2,4}};
2141 
2142  const int VQFPrismVert[][2] = {{0,0},{1,1},{4,4},{2,2},{3,3},{5,5}};
2143  const int VQFPrismFace[][2] = {{0,4},{0,2},{4,2},{0,2},{0,4},{2,4}};
2144 
2145  nfacemodes = (nummodesmax-2)*(nummodesmax-2);
2146  // Replace two Quad faces on every vertex
2147  for(int v = 0; v < 6; ++v)
2148  {
2149  for(int f = 0; f < 2; ++f)
2150  {
2151  for(int i = 0; i < nfacemodes; ++i)
2152  {
2153  val = (*maxRmat[eHexahedron])(
2154  vertMapMaxR[eHexahedron][VFHexVert[v][f]],
2155  faceMapMaxR[eHexahedron][VFHexFace[v][f]][i]);
2156  newmat->SetValue(vertMapMaxR[ePrism][VQFPrismVert[v][f]],
2157  faceMapMaxR[ePrism][VQFPrismFace[v][f]][i],val);
2158  }
2159  }
2160  }
2161 
2162  // Mapping of Hex Edge-Face mappings to Prism Edge-Face Mappings
2163  const int nadjface[] = {1,2,1,2,1,1,1,1,2};
2164  const int EFHexEdge[][2] = {{0,-1},{1,1},{2,-1},{3,3},{4,-1},{5,-1},{6,-1},{7,-1},{9,11}};
2165  const int EFHexFace[][2] = {{0,-1},{0,2},{0,-1},{0,4},{4,-1},{2,-1},{2,-1},{4,-1},{2,4}};
2166  const int EQFPrismEdge[][2] = {{0,-1},{1,1},{2,-1},{3,3},{4,-1},{5,-1},{6,-1},{7,-1},{8,8}};
2167  const int EQFPrismFace[][2] = {{0,-1},{0,2},{0,-1},{0,4},{4,-1},{2,-1},{2,-1},{4,-1},{2,4}};
2168 
2169  // all base edges are coupled to face 0
2170  nfacemodes = (nummodesmax-2)*(nummodesmax-2);
2171  for(int e = 0; e < 9; ++e)
2172  {
2173  for(int f = 0; f < nadjface[e]; ++f)
2174  {
2175  for(int i = 0; i < nummodesmax-2; ++i)
2176  {
2177  for(int j = 0; j < nfacemodes; ++j)
2178  {
2179  int edgemapid = edgeMapMaxR[eHexahedron][EFHexEdge[e][f]][i];
2180  int facemapid = faceMapMaxR[eHexahedron][EFHexFace[e][f]][j];
2181 
2182  val = (*maxRmat[eHexahedron])(edgemapid,facemapid);
2183 
2184  int edgemapid1 = edgeMapMaxR[ePrism][EQFPrismEdge[e][f]][i];
2185  int facemapid1 = faceMapMaxR[ePrism][EQFPrismFace[e][f]][j];
2186  newmat->SetValue(edgemapid1, facemapid1, val);
2187  }
2188  }
2189  }
2190  }
2191  }
2192 
2193  const int VFTetVert[] = {0,1,3,1,0,3};
2194  const int VFTetFace[] = {1,1,1,1,1,1};
2195  const int VTFPrismVert[] = {0,1,4,2,3,5};
2196  const int VTFPrismFace[] = {1,1,1,3,3,3};
2197 
2198  // Handle all triangular faces from tetrahedron
2199  nfacemodes = (nummodesmax-3)*(nummodesmax-2)/2;
2200  for(int v = 0; v < 6; ++v)
2201  {
2202  for(int i = 0; i < nfacemodes; ++i)
2203  {
2204  val = (*maxRmat[eTetrahedron])
2205  (vertMapMaxR[eTetrahedron][VFTetVert[v]],
2206  faceMapMaxR[eTetrahedron][VFTetFace[v]][i]);
2207 
2208  newmat->SetValue(vertMapMaxR[ePrism][VTFPrismVert[v]],
2209  faceMapMaxR[ePrism][VTFPrismFace[v]][i],val);
2210  }
2211  }
2212 
2213  // Mapping of Tet Edge-Face mappings to Prism Edge-Face Mappings
2214  const int EFTetEdge[] = {0,3,4,0,4,3};
2215  const int EFTetFace[] = {1,1,1,1,1,1};
2216  const int ETFPrismEdge[] = {0,4,5,2,6,7};
2217  const int ETFPrismFace[] = {1,1,1,3,3,3};
2218 
2219  // handle all edge to triangular faces from tetrahedron
2220  // (only 6 this time)
2221  nfacemodes = (nummodesmax-3)*(nummodesmax-2)/2;
2222  for(int e = 0; e < 6; ++e)
2223  {
2224  for(int i = 0; i < nummodesmax-2; ++i)
2225  {
2226  for(int j = 0; j < nfacemodes; ++j)
2227  {
2228  int edgemapid = edgeMapMaxR[eTetrahedron][EFTetEdge[e]][i];
2229  int facemapid = faceMapMaxR[eTetrahedron][EFTetFace[e]][j];
2230  val = (*maxRmat[eTetrahedron])(edgemapid,facemapid);
2231 
2232  newmat->SetValue(edgeMapMaxR[ePrism][ETFPrismEdge[e]][i],
2233  faceMapMaxR[ePrism][ETFPrismFace[e]][j],val);
2234  }
2235  }
2236  }
2237 
2238 
2239  DNekScalMatSharedPtr PrismR
2240  = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
2241  maxRmat[ePrism] = PrismR;
2242  }
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, and Nektar::LibUtilities::eTetrahedron.

Referenced by SetUpReferenceElements().

◆ ReSetTetMaxRMat()

void Nektar::MultiRegions::PreconditionerLowEnergy::ReSetTetMaxRMat ( int  nummodesmax,
LocalRegions::TetExpSharedPtr &  TetExp,
ShapeToDNekMap &  maxRmat,
ShapeToIntArrayMap &  vertMapMaxR,
ShapeToIntArrayArrayMap &  edgeMapMaxR,
ShapeToIntArrayArrayMap &  faceMapMaxR 
)
private

Definition at line 1981 of file PreconditionerLowEnergy.cpp.

1987  {
1988  boost::ignore_unused(faceMapMaxR);
1989 
1990  int nRows = TetExp->NumBndryCoeffs();
1991  NekDouble val;
1992  NekDouble zero = 0.0;
1993  DNekMatSharedPtr newmat = MemoryManager<DNekMat>::
1994  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1995 
1996  // copy existing system
1997  for(int i = 0; i < nRows; ++i)
1998  {
1999  for(int j = 0; j < nRows; ++j)
2000  {
2001  val = (*maxRmat[eTetrahedron])(i,j);
2002  newmat->SetValue(i,j,val);
2003  }
2004  }
2005 
2006  // The following lists specify the number of adjacent
2007  // edges to each vertex (nadj) then the Hex vert to
2008  // use for each pyramid ver in the vert-edge map (VEVert)
2009  // followed by the hex edge to use for each Tet edge
2010  // in the vert-edge map (VEEdge)
2011  const int VEHexVert[][4] = {{0,0,0},{1,1,1},{2,2,2},{4,5,6}};
2012  const int VEHexEdge[][4] = {{0,3,4},{0,1,5},{1,2,6},{4,5,6}};
2013  const int VETetEdge[][4] = {{0,2,3},{0,1,4},{1,2,5},{3,4,5}};
2014 
2015  // fill vertex to edge coupling
2016  for(int v = 0; v < 4; ++v)
2017  {
2018  for(int e = 0; e < 3; ++e)
2019  {
2020  for(int i = 0; i < nummodesmax-2; ++i)
2021  {
2022  // Note this is using wrong shape but gives
2023  // answer that seems to have correct error!
2024  val = (*maxRmat[eHexahedron])(
2025  vertMapMaxR[eHexahedron][VEHexVert[v][e]],
2026  edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
2027  newmat->SetValue(vertMapMaxR[eTetrahedron][v],
2028  edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i],
2029  val);
2030  }
2031  }
2032  }
2033 
2034  DNekScalMatSharedPtr TetR =
2035  MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
2036 
2037  maxRmat[eTetrahedron] = TetR;
2038  }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, and Nektar::LibUtilities::eTetrahedron.

Referenced by SetUpReferenceElements().

◆ SetupBlockTransformationMatrix()

void Nektar::MultiRegions::PreconditionerLowEnergy::SetupBlockTransformationMatrix ( void  )
private

Set a block transformation matrices for each element type. These are needed in routines that transform the schur complement matrix to and from the low energy basis.

Definition at line 946 of file PreconditionerLowEnergy.cpp.

947  {
948  std::shared_ptr<MultiRegions::ExpList>
949  expList=((m_linsys.lock())->GetLocMat()).lock();
950  LocalRegions::Expansion3DSharedPtr locExp;
951  LocalRegions::Expansion3DSharedPtr locExpSav;
952  map<int,int> EdgeSize;
953 
954  int n;
955 
956  std::map<ShapeType, DNekScalMatSharedPtr> maxRmat;
957  map<ShapeType, LocalRegions::Expansion3DSharedPtr > maxElmt;
958  map<ShapeType, Array<OneD, unsigned int> > vertMapMaxR;
959  map<ShapeType, Array<OneD, Array<OneD, unsigned int> > > edgeMapMaxR;
960 
961 
962  //Sets up reference element and builds transformation matrix for
963  // maximum polynomial order meshes
964  SetUpReferenceElements(maxRmat,maxElmt,vertMapMaxR,edgeMapMaxR);
965 
966  const Array<OneD,const unsigned int>& nbdry_size
967  = m_locToGloMap.lock()->GetNumLocalBndCoeffsPerPatch();
968 
969  int n_exp=expList->GetNumElmts();
970 
971  MatrixStorage blkmatStorage = eDIAGONAL;
972 
973  //Variants of R matrices required for low energy preconditioning
974  m_RBlk = MemoryManager<DNekBlkMat>
975  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
976  m_InvRBlk = MemoryManager<DNekBlkMat>
977  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
978 
979  DNekMatSharedPtr rmat, invrmat;
980 
981  int offset = 0;
982 
983  // Set up transformation matrices whilst checking to see if
984  // consecutive matrices are the same and if so reuse the
985  // matrices and store how many consecutive offsets there
986  // are
987  for(n=0; n < n_exp; ++n)
988  {
989  locExp = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
990 
991  ShapeType eltype = locExp->DetShapeType();
992 
993  int nbndcoeffs = locExp->NumBndryCoeffs();
994 
995  if(m_sameBlock.size() == 0)
996  {
997  rmat = ExtractLocMat(locExp,maxRmat[eltype],
998  maxElmt[eltype],
999  vertMapMaxR[eltype],
1000  edgeMapMaxR[eltype]);
1001  //Block R matrix
1002  m_RBlk->SetBlock(n, n, rmat);
1003 
1004  invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
1005  invrmat->Invert();
1006 
1007  //Block inverse R matrix
1008  m_InvRBlk->SetBlock(n, n, invrmat);
1009 
1010  m_sameBlock.push_back(pair<int,int>(1,nbndcoeffs));
1011  locExpSav = locExp;
1012  }
1013  else
1014  {
1015  bool reuse = true;
1016 
1017  // check to see if same as previous matrix and
1018  // reuse if we can
1019  for(int i = 0; i < 3; ++i)
1020  {
1021  if(locExpSav->GetBasis(i) != locExp->GetBasis(i))
1022  {
1023  reuse = false;
1024  break;
1025  }
1026  }
1027 
1028  if(reuse)
1029  {
1030  m_RBlk->SetBlock(n, n, rmat);
1031  m_InvRBlk->SetBlock(n, n, invrmat);
1032 
1033  m_sameBlock[offset] =
1034  (pair<int,int>(m_sameBlock[offset].first+1,nbndcoeffs));
1035  }
1036  else
1037  {
1038  rmat = ExtractLocMat(locExp,maxRmat[eltype],
1039  maxElmt[eltype],
1040  vertMapMaxR[eltype],
1041  edgeMapMaxR[eltype]);
1042 
1043  //Block R matrix
1044  m_RBlk->SetBlock(n, n, rmat);
1045 
1046  invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
1047  invrmat->Invert();
1048  //Block inverse R matrix
1049  m_InvRBlk->SetBlock(n, n, invrmat);
1050 
1051  m_sameBlock.push_back(pair<int,int>(1,nbndcoeffs));
1052  offset++;
1053  locExpSav = locExp;
1054  }
1055  }
1056  }
1057  }
Nektar::MultiRegions::Preconditioner::m_linsys
const std::weak_ptr< GlobalLinSys > m_linsys
Definition: Preconditioner.h:128
Nektar::MultiRegions::PreconditionerLowEnergy::ExtractLocMat
DNekMatSharedPtr ExtractLocMat(LocalRegions::Expansion3DSharedPtr &locExp, DNekScalMatSharedPtr &maxRmat, LocalRegions::Expansion3DSharedPtr &expMax, Array< OneD, unsigned int > &vertMapMaxR, Array< OneD, Array< OneD, unsigned int > > &edgeMapMaxR)
Definition: PreconditionerLowEnergy.cpp:2245
Nektar::MultiRegions::PreconditionerLowEnergy::m_sameBlock
std::vector< std::pair< int, int > > m_sameBlock
Definition: PreconditionerLowEnergy.h:89
Nektar::MultiRegions::PreconditionerLowEnergy::SetUpReferenceElements
void SetUpReferenceElements(ShapeToDNekMap &maxRmat, ShapeToExpMap &maxElmt, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR)
Loop expansion and determine different variants of the transformation matrix.
Definition: PreconditionerLowEnergy.cpp:1633
Nektar::LocalRegions::Expansion3DSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:49

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::eDIAGONAL, ExtractLocMat(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, and SetUpReferenceElements().

Referenced by v_InitObject().

◆ SetUpPyrMaxRMat()

void Nektar::MultiRegions::PreconditionerLowEnergy::SetUpPyrMaxRMat ( int  nummodesmax,
LocalRegions::PyrExpSharedPtr &  PyrExp,
ShapeToDNekMap &  maxRmat,
ShapeToIntArrayMap &  vertMapMaxR,
ShapeToIntArrayArrayMap &  edgeMapMaxR,
ShapeToIntArrayArrayMap &  faceMapMaxR 
)
private

Definition at line 1845 of file PreconditionerLowEnergy.cpp.

1851  {
1852  int nRows = PyrExp->NumBndryCoeffs();
1853  NekDouble val;
1854  NekDouble zero = 0.0;
1855  DNekMatSharedPtr newmat = MemoryManager<DNekMat>::
1856  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1857 
1858  // set diagonal to 1
1859  for(int i = 0; i < nRows; ++i)
1860  {
1861  newmat->SetValue(i,i,1.0);
1862  }
1863 
1864  // The following lists specify the number of adjacent
1865  // edges to each vertex (nadj) then the Hex vert to
1866  // use for each pyramid ver in the vert-edge map (VEVert)
1867  // followed by the hex edge to use for each pyramid edge
1868  // in the vert-edge map (VEEdge)
1869  const int nadjedge[] = {3,3,3,3,4};
1870  const int VEHexVert[][4] = {{0,0,0,-1},{1,1,1,-1},{2,2,2,-1},{3,3,3,-1},{4,5,6,7}};
1871  const int VEHexEdge[][4] = {{0,3,4,-1},{0,1,5,-1},{1,2,6,-1},{2,3,7,-1},{4,5,6,7}};
1872  const int VEPyrEdge[][4] = {{0,3,4,-1},{0,1,5,-1},{1,2,6,-1},{2,3,7,-1},{4,5,6,7}};
1873 
1874  // fill vertex to edge coupling
1875  for(int v = 0; v < 5; ++v)
1876  {
1877  for(int e = 0; e < nadjedge[v]; ++e)
1878  {
1879  for(int i = 0; i < nummodesmax-2; ++i)
1880  {
1881  // Note this is using wrong shape but gives
1882  // answer that seems to have correct error!
1883  val = (*maxRmat[eHexahedron])(
1884  vertMapMaxR[eHexahedron][VEHexVert[v][e]],
1885  edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
1886  newmat->SetValue(vertMapMaxR[ePyramid][v],
1887  edgeMapMaxR[ePyramid][VEPyrEdge[v][e]][i],val);
1888  }
1889  }
1890  }
1891 
1892  int nfacemodes;
1893  nfacemodes = (nummodesmax-2)*(nummodesmax-2);
1894  // First four verties are all adjacent to base face
1895  for(int v = 0; v < 4; ++v)
1896  {
1897  for(int i = 0; i < nfacemodes; ++i)
1898  {
1899  val = (*maxRmat[eHexahedron])(vertMapMaxR[eHexahedron][v],
1900  faceMapMaxR[eHexahedron][0][i]);
1901  newmat->SetValue(vertMapMaxR[ePyramid][v],
1902  faceMapMaxR[ePyramid][0][i],val);
1903  }
1904  }
1905 
1906 
1907  const int nadjface[] = {2,2,2,2,4};
1908  const int VFTetVert[][4] = {{0,0,-1,-1},{1,1,-1,-1},{2,2,-1,-1},{0,2,-1,-1},{3,3,3,3}};
1909  const int VFTetFace[][4] = {{1,3,-1,-1},{1,2,-1,-1},{2,3,-1,-1},{1,3,-1,-1},{1,2,1,2}};
1910  const int VFPyrFace[][4] = {{1,4,-1,-1},{1,2,-1,-1},{2,3,-1,-1},{3,4,-1,-1},{1,2,3,4}};
1911 
1912  // next handle all triangular faces from tetrahedron
1913  nfacemodes = (nummodesmax-3)*(nummodesmax-2)/2;
1914  for(int v = 0; v < 5; ++v)
1915  {
1916  for(int f = 0; f < nadjface[v]; ++f)
1917  {
1918  for(int i = 0; i < nfacemodes; ++i)
1919  {
1920  val = (*maxRmat[eTetrahedron])(vertMapMaxR[eTetrahedron][VFTetVert[v][f]],
1921  faceMapMaxR[eTetrahedron][VFTetFace[v][f]][i]);
1922  newmat->SetValue(vertMapMaxR[ePyramid][v],
1923  faceMapMaxR[ePyramid][VFPyrFace[v][f]][i],val);
1924  }
1925 
1926  }
1927  }
1928 
1929  // Edge to face coupling
1930  // all base edges are coupled to face 0
1931  nfacemodes = (nummodesmax-2)*(nummodesmax-2);
1932  for(int e = 0; e < 4; ++e)
1933  {
1934  for(int i = 0; i < nummodesmax-2; ++i)
1935  {
1936  for(int j = 0; j < nfacemodes; ++j)
1937  {
1938  int edgemapid = edgeMapMaxR[eHexahedron][e][i];
1939  int facemapid = faceMapMaxR[eHexahedron][0][j];
1940 
1941  val = (*maxRmat[eHexahedron])(edgemapid,facemapid);
1942  newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
1943  faceMapMaxR[ePyramid][0][j],val);
1944  }
1945 
1946  }
1947  }
1948 
1949  const int nadjface1[] = {1,1,1,1,2,2,2,2};
1950  const int EFTetEdge[][2] = {{0,-1},{1,-1},{0,-1},{2,-1},{3,3},{4,4},{5,5},{3,5}};
1951  const int EFTetFace[][2] = {{1,-1},{2,-1},{1,-1},{3,-1},{1,3},{1,2},{2,3},{1,3}};
1952  const int EFPyrFace[][2] = {{1,-1},{2,-1},{3,-1},{4,-1},{1,4},{1,2},{2,3},{3,4}};
1953 
1954  // next handle all triangular faces from tetrahedron
1955  nfacemodes = (nummodesmax-3)*(nummodesmax-2)/2;
1956  for(int e = 0; e < 8; ++e)
1957  {
1958  for(int f = 0; f < nadjface1[e]; ++f)
1959  {
1960  for(int i = 0; i < nummodesmax-2; ++i)
1961  {
1962  for(int j = 0; j < nfacemodes; ++j)
1963  {
1964  int edgemapid = edgeMapMaxR[eTetrahedron][EFTetEdge[e][f]][i];
1965  int facemapid = faceMapMaxR[eTetrahedron][EFTetFace[e][f]][j];
1966 
1967  val = (*maxRmat[eTetrahedron])(edgemapid,facemapid);
1968  newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
1969  faceMapMaxR[ePyramid][EFPyrFace[e][f]][j],val);
1970  }
1971  }
1972  }
1973  }
1974 
1975  DNekScalMatSharedPtr PyrR;
1976  PyrR = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
1977  maxRmat[ePyramid] =PyrR;
1978  }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePyramid, and Nektar::LibUtilities::eTetrahedron.

Referenced by SetUpReferenceElements().

◆ SetUpReferenceElements()

void Nektar::MultiRegions::PreconditionerLowEnergy::SetUpReferenceElements ( ShapeToDNekMap &  maxRmat,
ShapeToExpMap &  maxElmt,
ShapeToIntArrayMap &  vertMapMaxR,
ShapeToIntArrayArrayMap &  edgeMapMaxR 
)
private

Loop expansion and determine different variants of the transformation matrix. 

Sets up multiple reference elements based on the element expansion.

Definition at line 1632 of file PreconditionerLowEnergy.cpp.

1637  {
1638 
1639  std::shared_ptr<MultiRegions::ExpList>
1640  expList=((m_linsys.lock())->GetLocMat()).lock();
1641  GlobalLinSysKey linSysKey=(m_linsys.lock())->GetKey();
1642 
1643  LocalRegions::ExpansionSharedPtr locExp;
1644 
1645  // face maps for pyramid and hybrid meshes - not need to return.
1646  map<ShapeType, Array<OneD, Array<OneD, unsigned int> > > faceMapMaxR;
1647 
1648  /* Determine the maximum expansion order for all elements */
1649  int nummodesmax = 0;
1650  Array<OneD, int> Shapes(LibUtilities::SIZE_ShapeType,0);
1651 
1652  for(int n = 0; n < expList->GetNumElmts(); ++n)
1653  {
1654  locExp = expList->GetExp(n);
1655 
1656  nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(0));
1657  nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(1));
1658  nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(2));
1659 
1660  Shapes[locExp->DetShapeType()] = 1;
1661  }
1662 
1663 
1664  m_comm->AllReduce(nummodesmax, ReduceMax);
1665  m_comm->AllReduce(Shapes, ReduceMax);
1666 
1667  if(Shapes[ePyramid] || Shapes[ePrism]) // if Pyramids or Prisms used then need Tet and Hex expansion
1668  {
1669  Shapes[eTetrahedron] = 1;
1670  Shapes[eHexahedron] = 1;
1671  }
1672 
1673  StdRegions::MatrixType PreconR;
1674  if(linSysKey.GetMatrixType() == StdRegions::eMass)
1675  {
1676  PreconR = StdRegions::ePreconRMass;
1677  }
1678  else
1679  {
1680  PreconR = StdRegions::ePreconR;
1681  }
1682 
1683  Array<OneD, unsigned int> vmap;
1684  Array<OneD, Array<OneD, unsigned int> > emap;
1685  Array<OneD, Array<OneD, unsigned int> > fmap;
1686 
1687  /*
1688  * Set up a transformation matrices for equal max order
1689  * polynomial meshes
1690  */
1691 
1692  if(Shapes[eHexahedron])
1693  {
1694  SpatialDomains::HexGeomSharedPtr hexgeom = CreateRefHexGeom();
1695  //Bases for Hex element
1696  const BasisKey HexBa(eModified_A, nummodesmax,
1697  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1698  const BasisKey HexBb(eModified_A, nummodesmax,
1699  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1700  const BasisKey HexBc(eModified_A, nummodesmax,
1701  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1702 
1703  //Create reference Hexahdedral expansion
1704  LocalRegions::HexExpSharedPtr HexExp;
1705 
1706  HexExp = MemoryManager<LocalRegions::HexExp>
1707  ::AllocateSharedPtr(HexBa,HexBb,HexBc,
1708  hexgeom);
1709 
1710  maxElmt[eHexahedron] = HexExp;
1711 
1712  // Hex:
1713  HexExp->GetInverseBoundaryMaps(vmap,emap,fmap);
1714  vertMapMaxR[eHexahedron] = vmap;
1715  edgeMapMaxR[eHexahedron] = emap;
1716  faceMapMaxR[eHexahedron] = fmap;
1717 
1718  //Get hexahedral transformation matrix
1719  LocalRegions::MatrixKey HexR
1720  (PreconR, eHexahedron,
1721  *HexExp, linSysKey.GetConstFactors());
1722  maxRmat[eHexahedron] = HexExp->GetLocMatrix(HexR);
1723  }
1724 
1725  if(Shapes[eTetrahedron])
1726  {
1727  SpatialDomains::TetGeomSharedPtr tetgeom = CreateRefTetGeom();
1728  //Bases for Tetrahedral element
1729  const BasisKey TetBa(eModified_A, nummodesmax,
1730  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1731  const BasisKey TetBb(eModified_B, nummodesmax,
1732  PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
1733  const BasisKey TetBc(eModified_C, nummodesmax,
1734  PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
1735 
1736  //Create reference tetrahedral expansion
1737  LocalRegions::TetExpSharedPtr TetExp;
1738 
1739  TetExp = MemoryManager<LocalRegions::TetExp>
1740  ::AllocateSharedPtr(TetBa,TetBb,TetBc,tetgeom);
1741 
1742  maxElmt[eTetrahedron] = TetExp;
1743 
1744  TetExp->GetInverseBoundaryMaps(vmap,emap,fmap);
1745  vertMapMaxR[eTetrahedron] = vmap;
1746  edgeMapMaxR[eTetrahedron] = emap;
1747  faceMapMaxR[eTetrahedron] = fmap;
1748 
1749  //Get tetrahedral transformation matrix
1750  LocalRegions::MatrixKey TetR
1751  (PreconR, eTetrahedron,
1752  *TetExp, linSysKey.GetConstFactors());
1753  maxRmat[eTetrahedron] = TetExp->GetLocMatrix(TetR);
1754 
1755  if((Shapes[ePyramid])||(Shapes[eHexahedron]))
1756  {
1757  ReSetTetMaxRMat(nummodesmax, TetExp, maxRmat,
1758  vertMapMaxR, edgeMapMaxR, faceMapMaxR);
1759  }
1760  }
1761 
1762  if(Shapes[ePyramid])
1763  {
1764  SpatialDomains::PyrGeomSharedPtr pyrgeom = CreateRefPyrGeom();
1765 
1766  //Bases for Pyramid element
1767  const BasisKey PyrBa(eModified_A, nummodesmax,
1768  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1769  const BasisKey PyrBb(eModified_A, nummodesmax,
1770  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1771  const BasisKey PyrBc(eModifiedPyr_C, nummodesmax,
1772  PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
1773 
1774  //Create reference pyramid expansion
1775  LocalRegions::PyrExpSharedPtr PyrExp;
1776 
1777  PyrExp = MemoryManager<LocalRegions::PyrExp>
1778  ::AllocateSharedPtr(PyrBa,PyrBb,PyrBc,pyrgeom);
1779 
1780  maxElmt[ePyramid] = PyrExp;
1781 
1782  // Pyramid:
1783  PyrExp->GetInverseBoundaryMaps(vmap,emap,fmap);
1784  vertMapMaxR[ePyramid] = vmap;
1785  edgeMapMaxR[ePyramid] = emap;
1786  faceMapMaxR[ePyramid] = fmap;
1787 
1788  // Set up Pyramid Transformation Matrix based on Tet
1789  // and Hex expansion
1790  SetUpPyrMaxRMat(nummodesmax,PyrExp,maxRmat,vertMapMaxR,
1791  edgeMapMaxR,faceMapMaxR);
1792  }
1793 
1794  if(Shapes[ePrism])
1795  {
1796  SpatialDomains::PrismGeomSharedPtr prismgeom = CreateRefPrismGeom();
1797  //Bases for Prism element
1798  const BasisKey PrismBa(eModified_A, nummodesmax,
1799  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1800  const BasisKey PrismBb(eModified_A, nummodesmax,
1801  PointsKey(nummodesmax+1, eGaussLobattoLegendre));
1802  const BasisKey PrismBc(eModified_B, nummodesmax,
1803  PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
1804 
1805  //Create reference prismatic expansion
1806  LocalRegions::PrismExpSharedPtr PrismExp;
1807 
1808  PrismExp = MemoryManager<LocalRegions::PrismExp>
1809  ::AllocateSharedPtr(PrismBa,PrismBb,PrismBc,prismgeom);
1810  maxElmt[ePrism] = PrismExp;
1811 
1812  // Prism:
1813  PrismExp->GetInverseBoundaryMaps(vmap,emap,fmap);
1814  vertMapMaxR[ePrism] = vmap;
1815  edgeMapMaxR[ePrism] = emap;
1816 
1817  faceMapMaxR[ePrism] = fmap;
1818 
1819  if((Shapes[ePyramid])||(Shapes[eHexahedron]))
1820  {
1821  ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat,
1822  vertMapMaxR, edgeMapMaxR,
1823  faceMapMaxR, false);
1824  }
1825  else
1826  {
1827 
1828  //Get prismatic transformation matrix
1829  LocalRegions::MatrixKey PrismR
1830  (PreconR, ePrism,
1831  *PrismExp, linSysKey.GetConstFactors());
1832  maxRmat[ePrism] =
1833  PrismExp->GetLocMatrix(PrismR);
1834 
1835  if(Shapes[eTetrahedron]) // reset triangular faces from Tet
1836  {
1837  ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat,
1838  vertMapMaxR, edgeMapMaxR,
1839  faceMapMaxR, true);
1840  }
1841  }
1842  }
1843  }
Nektar::MultiRegions::Preconditioner::m_comm
LibUtilities::CommSharedPtr m_comm
Definition: Preconditioner.h:132
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefHexGeom
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1554
Nektar::MultiRegions::PreconditionerLowEnergy::ReSetPrismMaxRMat
void ReSetPrismMaxRMat(int nummodesmax, LocalRegions::PrismExpSharedPtr &PirsmExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR, bool UseTetOnly)
Definition: PreconditionerLowEnergy.cpp:2040
Nektar::MultiRegions::PreconditionerLowEnergy::ReSetTetMaxRMat
void ReSetTetMaxRMat(int nummodesmax, LocalRegions::TetExpSharedPtr &TetExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
Definition: PreconditionerLowEnergy.cpp:1981
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefTetGeom
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom(void)
Sets up the reference tretrahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1470
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPrismGeom
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1304
Nektar::MultiRegions::PreconditionerLowEnergy::SetUpPyrMaxRMat
void SetUpPyrMaxRMat(int nummodesmax, LocalRegions::PyrExpSharedPtr &PyrExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
Definition: PreconditionerLowEnergy.cpp:1845
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPyrGeom
SpatialDomains::PyrGeomSharedPtr CreateRefPyrGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis mapping arrays.
Definition: PreconditionerLowEnergy.cpp:1386
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
@ eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
@ eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LibUtilities::eModifiedPyr_C
@ eModifiedPyr_C
Principle Modified Functions .
Definition: BasisType.h:52
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::LocalRegions::PrismExpSharedPtr
std::shared_ptr< PrismExp > PrismExpSharedPtr
Definition: PrismExp.h:223
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::LocalRegions::HexExpSharedPtr
std::shared_ptr< HexExp > HexExpSharedPtr
Definition: HexExp.h:56
Nektar::LocalRegions::PyrExpSharedPtr
std::shared_ptr< PyrExp > PyrExpSharedPtr
Definition: PyrExp.h:195
Nektar::LocalRegions::TetExpSharedPtr
std::shared_ptr< TetExp > TetExpSharedPtr
Definition: TetExp.h:230

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), CreateRefHexGeom(), CreateRefPrismGeom(), CreateRefPyrGeom(), CreateRefTetGeom(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauMAlpha2Beta0, Nektar::LibUtilities::eHexahedron, Nektar::StdRegions::eMass, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eModifiedPyr_C, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eTetrahedron, Nektar::MultiRegions::GlobalMatrixKey::GetConstFactors(), Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), Nektar::MultiRegions::Preconditioner::m_comm, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::LibUtilities::ReduceMax, ReSetPrismMaxRMat(), ReSetTetMaxRMat(), SetUpPyrMaxRMat(), and Nektar::LibUtilities::SIZE_ShapeType.

Referenced by SetupBlockTransformationMatrix().

◆ v_BuildPreconditioner()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_BuildPreconditioner ( )
privatevirtual

Construct the low energy preconditioner from \(\mathbf{S}_{2}\).

\[\mathbf{M}^{-1}=\left[\begin{array}{ccc} Diag[(\mathbf{S_{2}})_{vv}] & & \\ & (\mathbf{S}_{2})_{eb} & \\ & & (\mathbf{S}_{2})_{fb} \end{array}\right] \]

where \(\mathbf{R}\) is the transformation matrix and \(\mathbf{S}_{2}\) the Schur complement of the modified basis, given by

\[\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\]

where \(\mathbf{S}_{1}\) is the local schur complement matrix for each element.

  • Count edges, face and add up min edges and min face sizes

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 121 of file PreconditionerLowEnergy.cpp.

122  {
123  std::shared_ptr<MultiRegions::ExpList>
124  expList=((m_linsys.lock())->GetLocMat()).lock();
125  LocalRegions::Expansion3DSharedPtr locExpansion;
126  GlobalLinSysKey linSysKey=(m_linsys.lock())->GetKey();
127 
128  int i, j, k;
129  int nVerts, nEdges,nFaces;
130  int eid, fid, n, cnt, nmodes, nedgemodes, nfacemodes;
131  int nedgemodesloc;
132  NekDouble zero = 0.0;
133 
134  int vMap1, vMap2, sign1, sign2;
135  int m, v, eMap1, eMap2, fMap1, fMap2;
136  int offset, globalrow, globalcol, nCoeffs;
137 
138  // Periodic information
139  PeriodicMap periodicVerts;
140  PeriodicMap periodicEdges;
141  PeriodicMap periodicFaces;
142  expList->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
143 
144  //matrix storage
145  MatrixStorage storage = eFULL;
146  MatrixStorage vertstorage = eDIAGONAL;
147  MatrixStorage blkmatStorage = eDIAGONAL;
148 
149  //local element static condensed matrices
150  DNekScalBlkMatSharedPtr loc_mat;
151  DNekScalMatSharedPtr bnd_mat;
152 
153  DNekMatSharedPtr pRSRT;
154 
155  DNekMat RS;
156  DNekMat RSRT;
157 
158  auto asmMap = m_locToGloMap.lock();
159 
160  int nDirBnd = asmMap->GetNumGlobalDirBndCoeffs();
161  int nNonDirVerts = asmMap->GetNumNonDirVertexModes();
162 
163  //Vertex, edge and face preconditioner matrices
164  DNekMatSharedPtr VertBlk = MemoryManager<DNekMat>::
165  AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
166 
167  Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
168  Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
169 
170  //maps for different element types
171  int n_exp = expList->GetNumElmts();
172  int nNonDirEdgeIDs=asmMap->GetNumNonDirEdges();
173  int nNonDirFaceIDs=asmMap->GetNumNonDirFaces();
174 
175  set<int> edgeDirMap;
176  set<int> faceDirMap;
177  map<int,int> uniqueEdgeMap;
178  map<int,int> uniqueFaceMap;
179 
180  //this should be of size total number of local edges + faces
181  Array<OneD, int> modeoffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs,0);
182  Array<OneD, int> globaloffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs,0);
183 
184  const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
185  LocalRegions::Expansion2DSharedPtr bndCondFaceExp;
186  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>&
187  bndConditions = expList->GetBndConditions();
188 
189  int meshVertId;
190  int meshEdgeId;
191  int meshFaceId;
192 
193  const Array<OneD, const int> &extradiredges
194  = asmMap->GetExtraDirEdges();
195  for(i=0; i<extradiredges.size(); ++i)
196  {
197  meshEdgeId=extradiredges[i];
198  edgeDirMap.insert(meshEdgeId);
199  }
200 
201  //Determine which boundary edges and faces have dirichlet values
202  for(i = 0; i < bndCondExp.size(); i++)
203  {
204  cnt = 0;
205  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
206  {
207  bndCondFaceExp = std::dynamic_pointer_cast<
208  LocalRegions::Expansion2D>(bndCondExp[i]->GetExp(j));
209  if (bndConditions[i]->GetBoundaryConditionType() ==
210  SpatialDomains::eDirichlet)
211  {
212  for(k = 0; k < bndCondFaceExp->GetNtraces(); k++)
213  {
214  meshEdgeId = bndCondFaceExp->GetGeom()->GetEid(k);
215  if(edgeDirMap.count(meshEdgeId) == 0)
216  {
217  edgeDirMap.insert(meshEdgeId);
218  }
219  }
220  meshFaceId = bndCondFaceExp->GetGeom()->GetGlobalID();
221  faceDirMap.insert(meshFaceId);
222  }
223  }
224  }
225 
226  int dof=0;
227  int maxFaceDof=0;
228  int maxEdgeDof=0;
229  int nlocalNonDirEdges=0;
230  int nlocalNonDirFaces=0;
231  int ntotalentries=0;
232 
233  map<int,int> EdgeSize;
234  map<int,int> FaceSize;
235  map<int,pair<int,int> >FaceModes;
236 
237  /// - Count edges, face and add up min edges and min face sizes
238  for(n = 0; n < n_exp; ++n)
239  {
240  locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
241 
242  nEdges = locExpansion->GetNedges();
243  for(j = 0; j < nEdges; ++j)
244  {
245  int nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
246  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
247  ->GetGeom3D()->GetEid(j);
248  if(EdgeSize.count(meshEdgeId) == 0)
249  {
250  EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
251  }
252  else
253  {
254  EdgeSize[meshEdgeId] = min(EdgeSize[meshEdgeId],
255  nEdgeInteriorCoeffs);
256  }
257  }
258 
259  nFaces = locExpansion->GetNtraces();
260  for(j = 0; j < nFaces; ++j)
261  {
262  int nFaceInteriorCoeffs = locExpansion->GetTraceIntNcoeffs(j);
263  meshFaceId = locExpansion->GetGeom3D()->GetFid(j);
264 
265  if(FaceSize.count(meshFaceId) == 0)
266  {
267  FaceSize[meshFaceId] = nFaceInteriorCoeffs;
268 
269  int m0,m1;
270  locExpansion->GetTraceNumModes(j,m0,m1,locExpansion->
271  GetTraceOrient(j));
272  FaceModes[meshFaceId] = pair<int,int>(m0,m1);
273  }
274  else
275  {
276  if(nFaceInteriorCoeffs < FaceSize[meshFaceId])
277  {
278  FaceSize[meshFaceId] = nFaceInteriorCoeffs;
279  int m0,m1;
280  locExpansion->GetTraceNumModes(j,m0,m1,locExpansion->
281  GetTraceOrient(j));
282  FaceModes[meshFaceId] = pair<int,int>(m0,m1);
283  }
284  }
285  }
286  }
287 
288  bool verbose =
289  expList->GetSession()->DefinesCmdLineArgument("verbose");
290 
291  // For parallel runs need to check have minimum of edges and faces over
292  // partition boundaries
293  if(m_comm->GetSize() > 1)
294  {
295  int EdgeSizeLen = EdgeSize.size();
296  int FaceSizeLen = FaceSize.size();
297  Array<OneD, long> FacetMap(EdgeSizeLen+FaceSizeLen,-1);
298  Array<OneD, NekDouble> FacetLen(EdgeSizeLen+FaceSizeLen,-1);
299 
300  map<int,int>::iterator it;
301 
302  cnt = 0;
303  int maxid = 0;
304  for(it = EdgeSize.begin(); it!=EdgeSize.end(); ++it,++cnt)
305  {
306  FacetMap[cnt] = it->first;
307  maxid = max(it->first,maxid);
308  FacetLen[cnt] = it->second;
309  }
310  maxid++;
311 
312  m_comm->AllReduce(maxid,ReduceMax);
313 
314  for(it = FaceSize.begin(); it!=FaceSize.end(); ++it,++cnt)
315  {
316  FacetMap[cnt] = it->first + maxid;
317  FacetLen[cnt] = it->second;
318  }
319 
320  //Exchange vertex data over different processes
321  Gs::gs_data *tmp = Gs::Init(FacetMap, m_comm, verbose);
322  Gs::Gather(FacetLen, Gs::gs_min, tmp);
323 
324  cnt = 0;
325  for(it = EdgeSize.begin(); it!=EdgeSize.end(); ++it,++cnt)
326  {
327  it->second = (int) FacetLen[cnt];
328  }
329 
330  for(it = FaceSize.begin(); it!=FaceSize.end(); ++it,++cnt)
331  {
332  it->second = (int)FacetLen[cnt];
333  }
334  }
335 
336  // Loop over all the elements in the domain and compute total edge
337  // DOF and set up unique ordering.
338  map<int,int> nblks;
339  int matrixlocation = 0;
340 
341  // First do periodic edges
342  for (auto &pIt : periodicEdges)
343  {
344  meshEdgeId = pIt.first;
345 
346  if(edgeDirMap.count(meshEdgeId)==0)
347  {
348  dof = EdgeSize[meshEdgeId];
349  if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
350  {
351  bool SetUpNewEdge = true;
352 
353 
354  for (i = 0; i < pIt.second.size(); ++i)
355  {
356  if (!pIt.second[i].isLocal)
357  {
358  continue;
359  }
360 
361  int meshEdgeId2 = pIt.second[i].id;
362  if(edgeDirMap.count(meshEdgeId2)==0)
363  {
364  if(uniqueEdgeMap.count(meshEdgeId2)!=0)
365  {
366  // set unique map to same location
367  uniqueEdgeMap[meshEdgeId] =
368  uniqueEdgeMap[meshEdgeId2];
369  SetUpNewEdge = false;
370  }
371  }
372  else
373  {
374  edgeDirMap.insert(meshEdgeId);
375  SetUpNewEdge = false;
376  }
377  }
378 
379  if(SetUpNewEdge)
380  {
381  uniqueEdgeMap[meshEdgeId]=matrixlocation;
382  globaloffset[matrixlocation]+=ntotalentries;
383  modeoffset[matrixlocation]=dof*dof;
384  ntotalentries+=dof*dof;
385  nblks [matrixlocation++] = dof;
386  }
387  }
388  }
389  }
390 
391  for(cnt=n=0; n < n_exp; ++n)
392  {
393  locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
394 
395  for (j = 0; j < locExpansion->GetNedges(); ++j)
396  {
397  meshEdgeId = locExpansion->GetGeom()->GetEid(j);
398  dof = EdgeSize[meshEdgeId];
399  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
400 
401  if(edgeDirMap.count(meshEdgeId)==0)
402  {
403  if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
404 
405  {
406  uniqueEdgeMap[meshEdgeId]=matrixlocation;
407 
408  globaloffset[matrixlocation]+=ntotalentries;
409  modeoffset[matrixlocation]=dof*dof;
410  ntotalentries+=dof*dof;
411  nblks[matrixlocation++] = dof;
412  }
413  nlocalNonDirEdges+=dof*dof;
414  }
415  }
416  }
417 
418  // Loop over all the elements in the domain and compute max face
419  // DOF. Reduce across all processes to get universal maximum.
420  // - Periodic faces
421  for (auto &pIt : periodicFaces)
422  {
423  meshFaceId = pIt.first;
424 
425  if(faceDirMap.count(meshFaceId)==0)
426  {
427  dof = FaceSize[meshFaceId];
428 
429  if(uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
430  {
431  bool SetUpNewFace = true;
432 
433  if(pIt.second[0].isLocal)
434  {
435  int meshFaceId2 = pIt.second[0].id;
436 
437  if(faceDirMap.count(meshFaceId2)==0)
438  {
439  if(uniqueFaceMap.count(meshFaceId2)!=0)
440  {
441  // set unique map to same location
442  uniqueFaceMap[meshFaceId] =
443  uniqueFaceMap[meshFaceId2];
444  SetUpNewFace = false;
445  }
446  }
447  else // set face to be a Dirichlet face
448  {
449  faceDirMap.insert(meshFaceId);
450  SetUpNewFace = false;
451  }
452  }
453 
454  if(SetUpNewFace)
455  {
456  uniqueFaceMap[meshFaceId]=matrixlocation;
457 
458  modeoffset[matrixlocation]=dof*dof;
459  globaloffset[matrixlocation]+=ntotalentries;
460  ntotalentries+=dof*dof;
461  nblks[matrixlocation++] = dof;
462  }
463  }
464  }
465  }
466 
467  for(cnt=n=0; n < n_exp; ++n)
468  {
469  locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
470 
471  for (j = 0; j < locExpansion->GetNtraces(); ++j)
472  {
473  meshFaceId = locExpansion->GetGeom()->GetFid(j);
474 
475  dof = FaceSize[meshFaceId];
476  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
477 
478  if(faceDirMap.count(meshFaceId)==0)
479  {
480  if(uniqueFaceMap.count(meshFaceId)==0 && dof > 0)
481  {
482  uniqueFaceMap[meshFaceId]=matrixlocation;
483  modeoffset[matrixlocation]=dof*dof;
484  globaloffset[matrixlocation]+=ntotalentries;
485  ntotalentries+=dof*dof;
486  nblks[matrixlocation++] = dof;
487  }
488  nlocalNonDirFaces+=dof*dof;
489  }
490  }
491  }
492 
493  m_comm->AllReduce(maxEdgeDof, ReduceMax);
494  m_comm->AllReduce(maxFaceDof, ReduceMax);
495 
496  //Allocate arrays for block to universal map (number of expansions * p^2)
497  Array<OneD, long> BlockToUniversalMap(ntotalentries,-1);
498  Array<OneD, int> localToGlobalMatrixMap(nlocalNonDirEdges +
499  nlocalNonDirFaces,-1);
500 
501  //Allocate arrays to store matrices (number of expansions * p^2)
502  Array<OneD, NekDouble> BlockArray(nlocalNonDirEdges +
503  nlocalNonDirFaces,0.0);
504 
505  int matrixoffset=0;
506  int vGlobal;
507  int uniEdgeOffset = 0;
508 
509  // Need to obtain a fixed offset for the universal number
510  // of the faces which come after the edge numbering
511  for(n=0; n < n_exp; ++n)
512  {
513  for(j = 0; j < locExpansion->GetNedges(); ++j)
514  {
515  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
516  ->GetGeom3D()->GetEid(j);
517 
518  uniEdgeOffset = max(meshEdgeId, uniEdgeOffset);
519  }
520  }
521  uniEdgeOffset++;
522 
523  m_comm->AllReduce(uniEdgeOffset,ReduceMax);
524  uniEdgeOffset = uniEdgeOffset*maxEdgeDof*maxEdgeDof;
525 
526  for(n=0; n < n_exp; ++n)
527  {
528  locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
529 
530  //loop over the edges of the expansion
531  for(j = 0; j < locExpansion->GetNedges(); ++j)
532  {
533  //get mesh edge id
534  meshEdgeId = locExpansion->GetGeom()->GetEid(j);
535 
536  nedgemodes = EdgeSize[meshEdgeId];
537 
538  if(edgeDirMap.count(meshEdgeId)==0 && nedgemodes > 0)
539  {
540  // Determine the Global edge offset
541  int edgeOffset = globaloffset[uniqueEdgeMap[meshEdgeId]];
542 
543  // Determine a universal map offset
544  int uniOffset = meshEdgeId;
545  auto pIt = periodicEdges.find(meshEdgeId);
546  if (pIt != periodicEdges.end())
547  {
548  for (int l = 0; l < pIt->second.size(); ++l)
549  {
550  uniOffset = min(uniOffset, pIt->second[l].id);
551  }
552  }
553  uniOffset = uniOffset*maxEdgeDof*maxEdgeDof;
554 
555  for(k=0; k<nedgemodes*nedgemodes; ++k)
556  {
557  vGlobal=edgeOffset+k;
558  localToGlobalMatrixMap[matrixoffset+k]=vGlobal;
559  BlockToUniversalMap[vGlobal] = uniOffset + k + 1;
560  }
561  matrixoffset+=nedgemodes*nedgemodes;
562  }
563  }
564 
565  Array<OneD, unsigned int> faceInteriorMap;
566  Array<OneD, int> faceInteriorSign;
567  //loop over the faces of the expansion
568  for(j = 0; j < locExpansion->GetNtraces(); ++j)
569  {
570  //get mesh face id
571  meshFaceId = locExpansion->GetGeom()->GetFid(j);
572 
573  nfacemodes = FaceSize[meshFaceId];
574 
575  //Check if face has dirichlet values
576  if(faceDirMap.count(meshFaceId)==0 && nfacemodes > 0)
577  {
578  // Determine the Global edge offset
579  int faceOffset = globaloffset[uniqueFaceMap[meshFaceId]];
580  // Determine a universal map offset
581  int uniOffset = meshFaceId;
582  // use minimum face edge when periodic
583  auto pIt = periodicFaces.find(meshFaceId);
584  if (pIt != periodicFaces.end())
585  {
586  uniOffset = min(uniOffset, pIt->second[0].id);
587  }
588  uniOffset = uniOffset * maxFaceDof * maxFaceDof;
589 
590  for(k=0; k<nfacemodes*nfacemodes; ++k)
591  {
592  vGlobal=faceOffset+k;
593 
594  localToGlobalMatrixMap[matrixoffset+k]
595  = vGlobal;
596 
597  BlockToUniversalMap[vGlobal] = uniOffset +
598  uniEdgeOffset + k + 1;
599  }
600  matrixoffset+=nfacemodes*nfacemodes;
601  }
602  }
603  }
604 
605  matrixoffset=0;
606 
607  map<int,int>::iterator it;
608  Array<OneD, unsigned int> n_blks(nblks.size()+1);
609  n_blks[0] = nNonDirVerts;
610  for(i =1, it = nblks.begin(); it != nblks.end(); ++it)
611  {
612  n_blks[i++] = it->second;
613  }
614 
615  m_BlkMat = MemoryManager<DNekBlkMat>
616  ::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
617 
618  //Here we loop over the expansion and build the block low energy
619  //preconditioner as well as the block versions of the transformation
620  //matrices.
621  for(cnt=n=0; n < n_exp; ++n)
622  {
623  locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
624  nCoeffs=locExpansion->NumBndryCoeffs();
625 
626  //Get correct transformation matrix for element type
627  DNekMat &R = (*m_RBlk->GetBlock(n,n));
628 
629  pRSRT = MemoryManager<DNekMat>::AllocateSharedPtr
630  (nCoeffs, nCoeffs, zero, storage);
631  RSRT = (*pRSRT);
632 
633  nVerts=locExpansion->GetGeom()->GetNumVerts();
634  nEdges=locExpansion->GetGeom()->GetNumEdges();
635  nFaces=locExpansion->GetGeom()->GetNumFaces();
636 
637  //Get statically condensed matrix
638  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
639 
640  //Extract boundary block (elemental S1)
641  bnd_mat=loc_mat->GetBlock(0,0);
642 
643  //offset by number of rows
644  offset = bnd_mat->GetRows();
645 
646  DNekScalMat &S=(*bnd_mat);
647 
648  DNekMat Sloc(nCoeffs,nCoeffs);
649 
650  // For variable p we need to just use the active modes
651  NekDouble val;
652 
653  for(int i = 0; i < nCoeffs; ++i)
654  {
655  for(int j = 0; j < nCoeffs; ++j)
656  {
657  val = S(i,j)*m_variablePmask[cnt+i]*m_variablePmask[cnt+j];
658  Sloc.SetValue(i,j,val);
659  }
660  }
661 
662  //Calculate R*S*trans(R)
663  RSRT = R*Sloc*Transpose(R);
664 
665  //loop over vertices of the element and return the vertex map
666  //for each vertex
667  for (v=0; v<nVerts; ++v)
668  {
669  vMap1=locExpansion->GetVertexMap(v);
670 
671  //Get vertex map
672  globalrow = asmMap->
673  GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
674 
675  if(globalrow >= 0)
676  {
677  for (m=0; m<nVerts; ++m)
678  {
679  vMap2=locExpansion->GetVertexMap(m);
680 
681  //global matrix location (without offset due to
682  //dirichlet values)
683  globalcol = asmMap->
684  GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
685 
686  //offset for dirichlet conditions
687  if (globalcol == globalrow)
688  {
689  //modal connectivity between elements
690  sign1 = asmMap->
691  GetLocalToGlobalBndSign(cnt + vMap1);
692  sign2 = asmMap->
693  GetLocalToGlobalBndSign(cnt + vMap2);
694 
695  vertArray[globalrow]
696  += sign1*sign2*RSRT(vMap1,vMap2);
697 
698 
699  meshVertId = locExpansion->GetGeom()->GetVid(v);
700 
701  auto pIt = periodicVerts.find(meshVertId);
702  if (pIt != periodicVerts.end())
703  {
704  for (k = 0; k < pIt->second.size(); ++k)
705  {
706  meshVertId = min(meshVertId, pIt->second[k].id);
707  }
708  }
709 
710  VertBlockToUniversalMap[globalrow]
711  = meshVertId + 1;
712  }
713  }
714  }
715  }
716 
717  //loop over edges of the element and return the edge map
718  for (eid=0; eid<nEdges; ++eid)
719  {
720 
721  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
722  ->GetGeom3D()->GetEid(eid);
723 
724 
725  nedgemodes = EdgeSize[meshEdgeId];
726  if(nedgemodes)
727  {
728  nedgemodesloc = locExpansion->GetEdgeNcoeffs(eid)-2;
729  DNekMatSharedPtr m_locMat =
730  MemoryManager<DNekMat>::AllocateSharedPtr
731  (nedgemodes,nedgemodes,zero,storage);
732 
733  Array<OneD, unsigned int> edgemodearray = locExpansion->GetEdgeInverseBoundaryMap(eid);
734 
735  if(edgeDirMap.count(meshEdgeId)==0)
736  {
737  for (v=0; v<nedgemodesloc; ++v)
738  {
739  eMap1=edgemodearray[v];
740  sign1 = asmMap->
741  GetLocalToGlobalBndSign(cnt + eMap1);
742 
743  if(sign1 == 0)
744  {
745  continue;
746  }
747 
748  for (m=0; m<nedgemodesloc; ++m)
749  {
750  eMap2=edgemodearray[m];
751 
752  //modal connectivity between elements
753  sign2 = asmMap->
754  GetLocalToGlobalBndSign(cnt + eMap2);
755 
756  NekDouble globalEdgeValue = sign1*sign2*RSRT(eMap1,eMap2);
757 
758  if(sign2 != 0)
759  {
760  //if(eMap1 == eMap2)
761  BlockArray[matrixoffset+v*nedgemodes+m]=globalEdgeValue;
762  }
763  }
764  }
765  matrixoffset+=nedgemodes*nedgemodes;
766  }
767  }
768  }
769 
770  //loop over faces of the element and return the face map
771  for (fid=0; fid<nFaces; ++fid)
772  {
773  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
774  ->GetGeom3D()->GetFid(fid);
775 
776  nfacemodes = FaceSize[meshFaceId];
777  if(nfacemodes > 0)
778  {
779  DNekMatSharedPtr m_locMat =
780  MemoryManager<DNekMat>::AllocateSharedPtr
781  (nfacemodes,nfacemodes,zero,storage);
782 
783  if(faceDirMap.count(meshFaceId) == 0)
784  {
785  Array<OneD, unsigned int> facemodearray;
786  StdRegions::Orientation faceOrient =
787  locExpansion->GetTraceOrient(fid);
788 
789  auto pIt = periodicFaces.find(meshFaceId);
790  if (pIt != periodicFaces.end())
791  {
792  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
793  {
794  faceOrient = DeterminePeriodicFaceOrient
795  (faceOrient,pIt->second[0].orient);
796  }
797  }
798 
799  facemodearray = locExpansion->GetTraceInverseBoundaryMap
800  (fid,faceOrient,FaceModes[meshFaceId].first,
801  FaceModes[meshFaceId].second);
802 
803  for (v=0; v<nfacemodes; ++v)
804  {
805  fMap1=facemodearray[v];
806 
807  sign1 = asmMap->
808  GetLocalToGlobalBndSign(cnt + fMap1);
809 
810  ASSERTL1(sign1 != 0,"Something is wrong since we "
811  "shoudl only be extracting modes for "
812  "lowest order expansion");
813 
814  for (m=0; m<nfacemodes; ++m)
815  {
816  fMap2=facemodearray[m];
817 
818  //modal connectivity between elements
819  sign2 = asmMap->
820  GetLocalToGlobalBndSign(cnt + fMap2);
821 
822  ASSERTL1(sign2 != 0,"Something is wrong since "
823  "we shoudl only be extracting modes for "
824  "lowest order expansion");
825 
826  // Get the face-face value from the
827  // low energy matrix (S2)
828  NekDouble globalFaceValue = sign1*sign2*
829  RSRT(fMap1,fMap2);
830 
831  //local face value to global face value
832  //if(fMap1 == fMap2)
833  BlockArray[matrixoffset+v*nfacemodes+m]=
834  globalFaceValue;
835  }
836  }
837  matrixoffset+=nfacemodes*nfacemodes;
838  }
839  }
840  }
841 
842  //offset for the expansion
843  cnt+=nCoeffs;
844  }
845 
846  if(nNonDirVerts!=0)
847  {
848  //Exchange vertex data over different processes
849  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm, verbose);
850  Gs::Gather(vertArray, Gs::gs_add, tmp);
851 
852  }
853 
854  Array<OneD, NekDouble> GlobalBlock(ntotalentries,0.0);
855  if(ntotalentries)
856  {
857  //Assemble edge matrices of each process
858  Vmath::Assmb(BlockArray.size(),
859  BlockArray,
860  localToGlobalMatrixMap,
861  GlobalBlock);
862  }
863 
864  //Exchange edge & face data over different processes
865  Gs::gs_data *tmp1 = Gs::Init(BlockToUniversalMap, m_comm, verbose);
866  Gs::Gather(GlobalBlock, Gs::gs_add, tmp1);
867 
868  // Populate vertex block
869  for (int i = 0; i < nNonDirVerts; ++i)
870  {
871  VertBlk->SetValue(i,i,1.0/vertArray[i]);
872  }
873 
874  //Set the first block to be the diagonal of the vertex space
875  m_BlkMat->SetBlock(0,0, VertBlk);
876 
877  //Build the edge and face matrices from the vector
878  DNekMatSharedPtr gmat;
879 
880  offset=0;
881  // -1 since we ignore vert block
882  for(int loc=0; loc<n_blks.size()-1; ++loc)
883  {
884  nmodes = n_blks[1+loc];
885  gmat = MemoryManager<DNekMat>::AllocateSharedPtr
886  (nmodes,nmodes,zero,storage);
887 
888  for (v=0; v<nmodes; ++v)
889  {
890  for (m=0; m<nmodes; ++m)
891  {
892  NekDouble Value = GlobalBlock[offset+v*nmodes+m];
893  gmat->SetValue(v,m,Value);
894 
895  }
896  }
897  m_BlkMat->SetBlock(1+loc,1+loc, gmat);
898  offset+=modeoffset[loc];
899  }
900 
901  // invert blocks.
902  int totblks=m_BlkMat->GetNumberOfBlockRows();
903  for (i=1; i< totblks; ++i)
904  {
905  unsigned int nmodes=m_BlkMat->GetNumberOfRowsInBlockRow(i);
906  if(nmodes)
907  {
908  DNekMatSharedPtr tmp_mat =
909  MemoryManager<DNekMat>::AllocateSharedPtr
910  (nmodes,nmodes,zero,storage);
911 
912  tmp_mat=m_BlkMat->GetBlock(i,i);
913  tmp_mat->Invert();
914 
915  m_BlkMat->SetBlock(i,i,tmp_mat);
916  }
917  }
918  }
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
CG_Iterations.loc
loc
Definition: CG_Iterations.py:323
Gs::Gather
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:245
Gs::gs_add
@ gs_add
Definition: GsLib.hpp:53
Gs::gs_min
@ gs_min
Definition: GsLib.hpp:53
Gs::Init
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:167
Nektar::LocalRegions::Expansion2DSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
Nektar::MultiRegions::PeriodicMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Definition: MultiRegions.hpp:207
Nektar::MultiRegions::DeterminePeriodicFaceOrient
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
Definition: AssemblyMapCG.cpp:2283
Nektar::DNekScalMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:68
Nektar::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:51
Nektar::Transpose
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:217
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:813
Gs::gs_data
Definition: GsLib.hpp:129

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Vmath::Assmb(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Gs::Gather(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Gs::gs_add, Gs::gs_min, Gs::Init(), CG_Iterations::loc, m_BlkMat, Nektar::MultiRegions::Preconditioner::m_comm, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_variablePmask, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), Nektar::LibUtilities::ReduceMax, and Nektar::Transpose().

◆ v_DoPreconditioner()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoPreconditioner ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
privatevirtual

Apply the low energy preconditioner during the conjugate gradient routine

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 925 of file PreconditionerLowEnergy.cpp.

928  {
929  int nDir = m_locToGloMap.lock()->GetNumGlobalDirBndCoeffs();
930  int nGlobal = m_locToGloMap.lock()->GetNumGlobalBndCoeffs();
931  int nNonDir = nGlobal-nDir;
932  DNekBlkMat &M = (*m_BlkMat);
933 
934  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
935  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
936 
937  z = M * r;
938  }
Nektar::DNekBlkMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:59

References Nektar::eWrapper, and Nektar::MultiRegions::Preconditioner::m_locToGloMap.

◆ v_DoTransformBasisFromLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformBasisFromLowEnergy ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
protectedvirtual

Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis.

Note; not typically required

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1145 of file PreconditionerLowEnergy.cpp.

1148  {
1149  int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1150 
1151  ASSERTL1(pInput.size() >= nLocBndDofs,
1152  "Input array is smaller than nLocBndDofs");
1153  ASSERTL1(pOutput.size() >= nLocBndDofs,
1154  "Output array is smaller than nLocBndDofs");
1155 
1156  //Block inverse transformation matrix
1157  DNekBlkMat &invR = *m_InvRBlk;
1158 
1159  Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
1160 
1161  //Multiply by the inverse transformation matrix
1162  int cnt = 0;
1163  int cnt1 = 0;
1164  for(int i = 0; i < m_sameBlock.size(); ++i)
1165  {
1166  int nexp = m_sameBlock[i].first;
1167  int nbndcoeffs = m_sameBlock[i].second;
1168  Blas::Dgemm('N','N', nbndcoeffs, nexp, nbndcoeffs,
1169  1.0, &(invR.GetBlock(cnt1,cnt1)->GetPtr()[0]),
1170  nbndcoeffs,pLocalIn.get() + cnt, nbndcoeffs,
1171  0.0, pOutput.get() + cnt, nbndcoeffs);
1172  cnt += nbndcoeffs*nexp;
1173  cnt1 += nexp;
1174  }
1175  }
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:394

References ASSERTL1, Blas::Dgemm(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and m_sameBlock.

◆ v_DoTransformBasisToLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformBasisToLowEnergy ( Array< OneD, NekDouble > &  pInOut)
protectedvirtual

Transform the basis vector to low energy space.

As the conjugate gradient system is solved for the low energy basis, the solution vector \(\mathbf{x}\) must be transformed to the low energy basis i.e. \(\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\).

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1066 of file PreconditionerLowEnergy.cpp.

1068  {
1069  int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1070 
1071  //Block transformation matrix
1072  DNekBlkMat &R = *m_RBlk;
1073 
1074  Array<OneD, NekDouble> pLocalIn(nLocBndDofs,pInOut.get());
1075 
1076  //Apply mask in case of variable P
1077  Vmath::Vmul(nLocBndDofs,pLocalIn, 1, m_variablePmask,1,
1078  pLocalIn,1);
1079 
1080  //Multiply by the block transformation matrix
1081  int cnt = 0;
1082  int cnt1 = 0;
1083  for(int i = 0; i < m_sameBlock.size(); ++i)
1084  {
1085  int nexp = m_sameBlock[i].first;
1086  int nbndcoeffs = m_sameBlock[i].second;
1087  Blas::Dgemm('N','N', nbndcoeffs, nexp, nbndcoeffs,
1088  1.0, &(R.GetBlock(cnt1,cnt1)->GetPtr()[0]),
1089  nbndcoeffs,pLocalIn.get() + cnt, nbndcoeffs,
1090  0.0, pInOut.get() + cnt, nbndcoeffs);
1091  cnt += nbndcoeffs*nexp;
1092  cnt1 += nexp;
1093  }
1094  }
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:192

References Blas::Dgemm(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, m_variablePmask, and Vmath::Vmul().

◆ v_DoTransformCoeffsFromLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformCoeffsFromLowEnergy ( Array< OneD, NekDouble > &  pInOut)
protectedvirtual

transform the solution coeffiicents from low energy back to the original coefficient space.

After the conjugate gradient routine the output vector is in the low energy basis and must be trasnformed back to the original basis in order to get the correct solution out. the solution vector i.e. \(\mathbf{x}=\mathbf{R^{T}}\mathbf{\overline{x}}\).

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1105 of file PreconditionerLowEnergy.cpp.

1107  {
1108  int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1109 
1110  ASSERTL1(pInOut.size() >= nLocBndDofs,
1111  "Output array is not greater than the nLocBndDofs");
1112 
1113  //Block transposed transformation matrix
1114  DNekBlkMat &R = *m_RBlk;
1115 
1116  Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
1117 
1118  //Multiply by the transpose of block transformation matrix
1119  int cnt = 0;
1120  int cnt1 = 0;
1121  for(int i = 0; i < m_sameBlock.size(); ++i)
1122  {
1123  int nexp = m_sameBlock[i].first;
1124  int nbndcoeffs = m_sameBlock[i].second;
1125  Blas::Dgemm('T','N', nbndcoeffs, nexp, nbndcoeffs,
1126  1.0, &(R.GetBlock(cnt1,cnt1)->GetPtr()[0]),
1127  nbndcoeffs,pLocalIn.get() + cnt, nbndcoeffs,
1128  0.0, pInOut.get() + cnt, nbndcoeffs);
1129  cnt += nbndcoeffs*nexp;
1130  cnt1 += nexp;
1131 
1132  }
1133 
1134  Vmath::Vmul(nLocBndDofs,pInOut, 1, m_variablePmask,1,
1135  pInOut,1);
1136  }

References ASSERTL1, Blas::Dgemm(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, m_variablePmask, and Vmath::Vmul().

◆ v_DoTransformCoeffsToLowEnergy()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformCoeffsToLowEnergy ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
protectedvirtual

Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space.

In JCP 2001 paper on low energy this is seen as (C^T)^{-1}

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1184 of file PreconditionerLowEnergy.cpp.

1187  {
1188  int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1189 
1190  ASSERTL1(pInput.size() >= nLocBndDofs,
1191  "Input array is less than nLocBndDofs");
1192  ASSERTL1(pOutput.size() >= nLocBndDofs,
1193  "Output array is less than nLocBndDofs");
1194 
1195  //Block inverse transformation matrix
1196  DNekBlkMat &invR = *m_InvRBlk;
1197 
1198  Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
1199 
1200  //Multiply by the transpose of block transformation matrix
1201  int cnt = 0;
1202  int cnt1 = 0;
1203  for(int i = 0; i < m_sameBlock.size(); ++i)
1204  {
1205  int nexp = m_sameBlock[i].first;
1206  int nbndcoeffs = m_sameBlock[i].second;
1207  Blas::Dgemm('T','N', nbndcoeffs, nexp, nbndcoeffs,
1208  1.0, &(invR.GetBlock(cnt1,cnt1)->GetPtr()[0]),
1209  nbndcoeffs,pLocalIn.get() + cnt, nbndcoeffs,
1210  0.0, pOutput.get() + cnt, nbndcoeffs);
1211  cnt += nbndcoeffs*nexp;
1212  cnt1 += nexp;
1213  }
1214  }

References ASSERTL1, Blas::Dgemm(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and m_sameBlock.

◆ v_InitObject()

void Nektar::MultiRegions::PreconditionerLowEnergy::v_InitObject ( )
privatevirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 74 of file PreconditionerLowEnergy.cpp.

75  {
76  GlobalSysSolnType solvertype =
77  m_locToGloMap.lock()->GetGlobalSysSolnType();
78 
79  ASSERTL0(solvertype == eIterativeStaticCond ||
80  solvertype == ePETScStaticCond, "Solver type not valid");
81 
82  std::shared_ptr<MultiRegions::ExpList>
83  expList=((m_linsys.lock())->GetLocMat()).lock();
84 
85  m_comm = expList->GetComm();
86 
87  LocalRegions::ExpansionSharedPtr locExpansion;
88 
89  locExpansion = expList->GetExp(0);
90 
91  int nDim = locExpansion->GetShapeDimension();
92 
93  ASSERTL0(nDim==3,
94  "Preconditioner type only valid in 3D");
95 
96  //Set up block transformation matrix
97  SetupBlockTransformationMatrix();
98 
99  //Sets up variable p mask
100  CreateVariablePMask();
101  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216

References ASSERTL0, CreateVariablePMask(), Nektar::MultiRegions::eIterativeStaticCond, Nektar::MultiRegions::ePETScStaticCond, Nektar::MultiRegions::Preconditioner::m_comm, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and SetupBlockTransformationMatrix().

◆ v_TransformedSchurCompl()

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::v_TransformedSchurCompl ( int  n,
int  bndoffset,
const std::shared_ptr< DNekScalMat > &  loc_mat 
)
privatevirtual

Set up the transformed block matrix system.

Sets up a block elemental matrix in which each of the block matrix is the low energy equivalent i.e. \(\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\)

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1223 of file PreconditionerLowEnergy.cpp.

1227  {
1228  std::shared_ptr<MultiRegions::ExpList>
1229  expList=((m_linsys.lock())->GetLocMat()).lock();
1230 
1231  LocalRegions::ExpansionSharedPtr locExpansion;
1232  locExpansion = expList->GetExp(n);
1233  unsigned int nbnd=locExpansion->NumBndryCoeffs();
1234 
1235  MatrixStorage storage = eFULL;
1236  DNekMatSharedPtr pS2 = MemoryManager<DNekMat>::
1237  AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1238 
1239  //transformation matrices
1240  DNekMat &R = (*m_RBlk->GetBlock(n,n));
1241 
1242  // Original Schur Complement
1243  DNekScalMat &S1 = (*loc_mat);
1244 
1245  DNekMat Sloc(nbnd,nbnd);
1246 
1247  // For variable p we need to just use the active modes
1248  NekDouble val;
1249 
1250  for(int i = 0; i < nbnd; ++i)
1251  {
1252  for(int j = 0; j < nbnd; ++j)
1253  {
1254  val = S1(i,j)*m_variablePmask[bndoffset+i]*
1255  m_variablePmask[bndoffset+j];
1256  Sloc.SetValue(i,j,val);
1257  }
1258  }
1259 
1260  //create low energy matrix
1261  DNekMat &S2 = (*pS2);
1262 
1263  S2= R*Sloc*Transpose(R);
1264 
1265  DNekScalMatSharedPtr return_val;
1266  return_val = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
1267 
1268  return return_val;
1269  }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::MultiRegions::Preconditioner::m_linsys, m_RBlk, m_variablePmask, and Nektar::Transpose().

Member Data Documentation

◆ className

string Nektar::MultiRegions::PreconditionerLowEnergy::className
static
Initial value:
= GetPreconFactory().RegisterCreatorFunction(
"LowEnergyBlock",
PreconditionerLowEnergy::create,
"LowEnergy Preconditioning")
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::MultiRegions::PreconditionerLowEnergy::create
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
Definition: PreconditionerLowEnergy.h:57
Nektar::MultiRegions::GetPreconFactory
PreconFactory & GetPreconFactory()
Definition: Preconditioner.cpp:89

Name of class.

Registers the class with the Factory.

Definition at line 67 of file PreconditionerLowEnergy.h.

◆ m_BlkMat

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_BlkMat
protected

Definition at line 78 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

◆ m_InvRBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_InvRBlk
protected

Definition at line 80 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_DoTransformBasisFromLowEnergy(), and v_DoTransformCoeffsToLowEnergy().

◆ m_RBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RBlk
protected

Definition at line 79 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_BuildPreconditioner(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_TransformedSchurCompl().

◆ m_sameBlock

std::vector<std::pair<int,int> > Nektar::MultiRegions::PreconditionerLowEnergy::m_sameBlock
protected

Definition at line 89 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_DoTransformBasisFromLowEnergy(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_DoTransformCoeffsToLowEnergy().

◆ m_signChange

bool Nektar::MultiRegions::PreconditionerLowEnergy::m_signChange
protected

Definition at line 85 of file PreconditionerLowEnergy.h.

Referenced by CreateVariablePMask().

◆ m_variablePmask

Array<OneD, NekDouble> Nektar::MultiRegions::PreconditionerLowEnergy::m_variablePmask
protected

Definition at line 83 of file PreconditionerLowEnergy.h.

Referenced by CreateVariablePMask(), v_BuildPreconditioner(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_TransformedSchurCompl().

Generated on Wed Nov 24 2021 10:33:36 for Nektar++ by doxygen 1.9.1