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)
 
 ~PreconditionerLowEnergy () override
 
- 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, const bool &IsLocal=false)
 
void DoAssembleLoc (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &ZeroDir)
 Apply an assembly and scatter back to lcoal array. More...
 
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

void v_InitObject () override
 
void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false) override
 
void v_BuildPreconditioner () override
 Construct the low energy preconditioner from \(\mathbf{S}_{2}\). More...
 
DNekScalMatSharedPtr v_TransformedSchurCompl (int n, int offset, const std::shared_ptr< DNekScalMat > &loc_mat) override
 Set up the transformed block matrix system. More...
 
void v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut) override
 Transform the basis vector to low energy space. More...
 
void v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut) override
 transform the solution coeffiicents from low energy back to the original coefficient space. More...
 
void v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
 Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis. More...
 
void v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
 Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space. More...
 
- Protected Member Functions inherited from Nektar::MultiRegions::Preconditioner
virtual DNekScalMatSharedPtr v_TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat)
 Get block elemental transposed transformation matrix \(\mathbf{R}^{T}\). More...
 
virtual void v_InitObject ()
 
virtual void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false)
 Apply a preconditioner to the conjugate gradient method. More...
 
virtual void v_DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce)
 Apply a preconditioner to the conjugate gradient method with an output for non-vertex degrees of freedom. More...
 
virtual void v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)
 Transform from original basis to low energy basis. More...
 
virtual void v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 Transform from low energy coeffs to orignal basis. More...
 
virtual void v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block inverse transformation matrix. More...
 
virtual void v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block transposed inverse transformation matrix. More...
 
virtual void v_BuildPreconditioner ()
 

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
 
std::string 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...
 

Detailed Description

This class implements low energy preconditioning for the conjugate gradient matrix solver.

Definition at line 51 of file PreconditionerLowEnergy.h.

Member Typedef Documentation

◆ ShapeToDNekMap

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

Definition at line 120 of file PreconditionerLowEnergy.h.

◆ ShapeToExpMap

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

Definition at line 123 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayArrayMap

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

Definition at line 128 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayMap

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

Definition at line 125 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

◆ PreconditionerLowEnergy()

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

Definition at line 65 of file PreconditionerLowEnergy.cpp.

68 : Preconditioner(plinsys, pLocToGloMap)
69{
70}
Nektar::MultiRegions::Preconditioner::Preconditioner
Preconditioner(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
Definition: Preconditioner.cpp:54

◆ ~PreconditionerLowEnergy()

Nektar::MultiRegions::PreconditionerLowEnergy::~PreconditionerLowEnergy ( )
inlineoverride

Definition at line 74 of file PreconditionerLowEnergy.h.

75 {
76 }

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 55 of file PreconditionerLowEnergy.h.

58 {
59 PreconditionerSharedPtr p =
60 MemoryManager<PreconditionerLowEnergy>::AllocateSharedPtr(
61 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:166
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::MultiRegions::PreconditionerSharedPtr
std::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:58

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 1542 of file PreconditionerLowEnergy.cpp.

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

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 1279 of file PreconditionerLowEnergy.cpp.

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

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 1376 of file PreconditionerLowEnergy.cpp.

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

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 1466 of file PreconditionerLowEnergy.cpp.

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

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

Referenced by SetUpReferenceElements().

◆ CreateVariablePMask()

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

Create the mask array

Definition at line 1252 of file PreconditionerLowEnergy.cpp.

1253{
1254 auto asmMap = m_locToGloMap.lock();
1255 unsigned int nLocBnd = asmMap->GetNumLocalBndCoeffs();
1256
1257 m_signChange = asmMap->GetSignChange();
1258
1259 // Construct a mask array
1260 m_variablePmask = Array<OneD, NekDouble>(nLocBnd, 1.0);
1261
1262 if (m_signChange)
1263 {
1264 unsigned int i;
1265 const Array<OneD, const NekDouble> &sign =
1266 asmMap->GetLocalToGlobalBndSign();
1267
1268 for (i = 0; i < nLocBnd; ++i)
1269 {
1270 m_variablePmask[i] = fabs(sign[i]);
1271 }
1272 }
1273}
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
Nektar::MultiRegions::Preconditioner::m_locToGloMap
std::weak_ptr< AssemblyMap > m_locToGloMap
Definition: Preconditioner.h:123

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 2264 of file PreconditionerLowEnergy.cpp.

2268{
2269 NekDouble val;
2270 NekDouble zero = 0.0;
2271
2272 int nRows = locExp->NumBndryCoeffs();
2273 DNekMatSharedPtr newmat =
2274 MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
2275
2276 Array<OneD, unsigned int> vlocmap;
2277 Array<OneD, Array<OneD, unsigned int>> elocmap;
2278 Array<OneD, Array<OneD, unsigned int>> flocmap;
2279
2280 locExp->GetInverseBoundaryMaps(vlocmap, elocmap, flocmap);
2281
2282 // fill diagonal
2283 for (int i = 0; i < nRows; ++i)
2284 {
2285 val = 1.0;
2286 newmat->SetValue(i, i, val);
2287 }
2288
2289 int nverts = locExp->GetNverts();
2290 int nedges = locExp->GetNedges();
2291 int nfaces = locExp->GetNtraces();
2292
2293 // fill vertex to edge coupling
2294 for (int e = 0; e < nedges; ++e)
2295 {
2296 int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
2297
2298 for (int v = 0; v < nverts; ++v)
2299 {
2300 for (int i = 0; i < nEdgeInteriorCoeffs; ++i)
2301 {
2302 val = (*maxRmat)(vmap[v], emap[e][i]);
2303 newmat->SetValue(vlocmap[v], elocmap[e][i], val);
2304 }
2305 }
2306 }
2307
2308 for (int f = 0; f < nfaces; ++f)
2309 {
2310 // Get details to extrac this face from max reference matrix
2311 StdRegions::Orientation FwdOrient =
2312 StdRegions::eDir1FwdDir1_Dir2FwdDir2;
2313 int m0, m1; // Local face expansion orders.
2314
2315 int nFaceInteriorCoeffs = locExp->GetTraceIntNcoeffs(f);
2316
2317 locExp->GetTraceNumModes(f, m0, m1, FwdOrient);
2318
2319 Array<OneD, unsigned int> fmapRmat =
2320 maxExp->GetTraceInverseBoundaryMap(f, FwdOrient, m0, m1);
2321
2322 // fill in vertex to face coupling
2323 for (int v = 0; v < nverts; ++v)
2324 {
2325 for (int i = 0; i < nFaceInteriorCoeffs; ++i)
2326 {
2327 val = (*maxRmat)(vmap[v], fmapRmat[i]);
2328 newmat->SetValue(vlocmap[v], flocmap[f][i], val);
2329 }
2330 }
2331
2332 // fill in edges to face coupling
2333 for (int e = 0; e < nedges; ++e)
2334 {
2335 int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
2336
2337 for (int j = 0; j < nEdgeInteriorCoeffs; ++j)
2338 {
2339 for (int i = 0; i < nFaceInteriorCoeffs; ++i)
2340 {
2341 val = (*maxRmat)(emap[e][j], fmapRmat[i]);
2342 newmat->SetValue(elocmap[e][j], flocmap[f][i], val);
2343 }
2344 }
2345 }
2346 }
2347
2348 return newmat;
2349}
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
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 2048 of file PreconditionerLowEnergy.cpp.

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

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 1990 of file PreconditionerLowEnergy.cpp.

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

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 936 of file PreconditionerLowEnergy.cpp.

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

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 1825 of file PreconditionerLowEnergy.cpp.

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

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 1619 of file PreconditionerLowEnergy.cpp.

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), CreateRefHexGeom(), CreateRefPrismGeom(), CreateRefPyrGeom(), CreateRefTetGeom(), Nektar::LibUtilities::eGaussLobattoLegendre, 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 ( )
overrideprotectedvirtual

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 117 of file PreconditionerLowEnergy.cpp.

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

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,
const bool &  isLocal = false 
)
overrideprotectedvirtual

Apply the low energy preconditioner during the conjugate gradient routine

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 912 of file PreconditionerLowEnergy.cpp.

916{
917 ASSERTL0(isLocal == false, "PreconditionerLowEnergy is only currently "
918 "set up for Global iteratives sovles");
919 int nDir = m_locToGloMap.lock()->GetNumGlobalDirBndCoeffs();
920 int nGlobal = m_locToGloMap.lock()->GetNumGlobalBndCoeffs();
921 int nNonDir = nGlobal - nDir;
922 DNekBlkMat &M = (*m_BlkMat);
923
924 NekVector<NekDouble> r(nNonDir, pInput, eWrapper);
925
926 NekVector<NekDouble> z(nNonDir, pOutput, eWrapper);
927
928 z = M * r;
929}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar::DNekBlkMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:58

References ASSERTL0, Nektar::eWrapper, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and Nektar::UnitTests::z().

◆ v_DoTransformBasisFromLowEnergy()

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

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 1127 of file PreconditionerLowEnergy.cpp.

1129{
1130 int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1131
1132 ASSERTL1(pInput.size() >= nLocBndDofs,
1133 "Input array is smaller than nLocBndDofs");
1134 ASSERTL1(pOutput.size() >= nLocBndDofs,
1135 "Output array is smaller than nLocBndDofs");
1136
1137 // Block inverse transformation matrix
1138 DNekBlkMat &invR = *m_InvRBlk;
1139
1140 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
1141
1142 // Multiply by the inverse transformation matrix
1143 int cnt = 0;
1144 int cnt1 = 0;
1145 for (int i = 0; i < m_sameBlock.size(); ++i)
1146 {
1147 int nexp = m_sameBlock[i].first;
1148 int nbndcoeffs = m_sameBlock[i].second;
1149 Blas::Dgemm('N', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
1150 &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1151 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pOutput.get() + cnt,
1152 nbndcoeffs);
1153 cnt += nbndcoeffs * nexp;
1154 cnt1 += nexp;
1155 }
1156}
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:383

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)
overrideprotectedvirtual

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 1051 of file PreconditionerLowEnergy.cpp.

1053{
1054 int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1055
1056 // Block transformation matrix
1057 DNekBlkMat &R = *m_RBlk;
1058
1059 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
1060
1061 // Apply mask in case of variable P
1062 Vmath::Vmul(nLocBndDofs, pLocalIn, 1, m_variablePmask, 1, pLocalIn, 1);
1063
1064 // Multiply by the block transformation matrix
1065 int cnt = 0;
1066 int cnt1 = 0;
1067 for (int i = 0; i < m_sameBlock.size(); ++i)
1068 {
1069 int nexp = m_sameBlock[i].first;
1070 int nbndcoeffs = m_sameBlock[i].second;
1071 Blas::Dgemm('N', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
1072 &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1073 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pInOut.get() + cnt,
1074 nbndcoeffs);
1075 cnt += nbndcoeffs * nexp;
1076 cnt1 += nexp;
1077 }
1078}
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.hpp:72

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)
overrideprotectedvirtual

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 1089 of file PreconditionerLowEnergy.cpp.

1091{
1092 int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1093
1094 ASSERTL1(pInOut.size() >= nLocBndDofs,
1095 "Output array is not greater than the nLocBndDofs");
1096
1097 // Block transposed transformation matrix
1098 DNekBlkMat &R = *m_RBlk;
1099
1100 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
1101
1102 // Multiply by the transpose of block transformation matrix
1103 int cnt = 0;
1104 int cnt1 = 0;
1105 for (int i = 0; i < m_sameBlock.size(); ++i)
1106 {
1107 int nexp = m_sameBlock[i].first;
1108 int nbndcoeffs = m_sameBlock[i].second;
1109 Blas::Dgemm('T', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
1110 &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1111 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pInOut.get() + cnt,
1112 nbndcoeffs);
1113 cnt += nbndcoeffs * nexp;
1114 cnt1 += nexp;
1115 }
1116
1117 Vmath::Vmul(nLocBndDofs, pInOut, 1, m_variablePmask, 1, pInOut, 1);
1118}

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 
)
overrideprotectedvirtual

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 1165 of file PreconditionerLowEnergy.cpp.

1167{
1168 int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
1169
1170 ASSERTL1(pInput.size() >= nLocBndDofs,
1171 "Input array is less than nLocBndDofs");
1172 ASSERTL1(pOutput.size() >= nLocBndDofs,
1173 "Output array is less than nLocBndDofs");
1174
1175 // Block inverse transformation matrix
1176 DNekBlkMat &invR = *m_InvRBlk;
1177
1178 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
1179
1180 // Multiply by the transpose of block transformation matrix
1181 int cnt = 0;
1182 int cnt1 = 0;
1183 for (int i = 0; i < m_sameBlock.size(); ++i)
1184 {
1185 int nexp = m_sameBlock[i].first;
1186 int nbndcoeffs = m_sameBlock[i].second;
1187 Blas::Dgemm('T', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
1188 &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1189 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pOutput.get() + cnt,
1190 nbndcoeffs);
1191 cnt += nbndcoeffs * nexp;
1192 cnt1 += nexp;
1193 }
1194}

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

◆ v_InitObject()

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 72 of file PreconditionerLowEnergy.cpp.

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

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 
)
overrideprotectedvirtual

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 1203 of file PreconditionerLowEnergy.cpp.

1205{
1206 std::shared_ptr<MultiRegions::ExpList> expList =
1207 ((m_linsys.lock())->GetLocMat()).lock();
1208
1209 LocalRegions::ExpansionSharedPtr locExpansion;
1210 locExpansion = expList->GetExp(n);
1211 unsigned int nbnd = locExpansion->NumBndryCoeffs();
1212
1213 MatrixStorage storage = eFULL;
1214 DNekMatSharedPtr pS2 =
1215 MemoryManager<DNekMat>::AllocateSharedPtr(nbnd, nbnd, 0.0, storage);
1216
1217 // transformation matrices
1218 DNekMat &R = (*m_RBlk->GetBlock(n, n));
1219
1220 // Original Schur Complement
1221 DNekScalMat &S1 = (*loc_mat);
1222
1223 DNekMat Sloc(nbnd, nbnd);
1224
1225 // For variable p we need to just use the active modes
1226 NekDouble val;
1227
1228 for (int i = 0; i < nbnd; ++i)
1229 {
1230 for (int j = 0; j < nbnd; ++j)
1231 {
1232 val = S1(i, j) * m_variablePmask[bndoffset + i] *
1233 m_variablePmask[bndoffset + j];
1234 Sloc.SetValue(i, j, val);
1235 }
1236 }
1237
1238 // create low energy matrix
1239 DNekMat &S2 = (*pS2);
1240
1241 S2 = R * Sloc * Transpose(R);
1242
1243 DNekScalMatSharedPtr return_val;
1244 return_val = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
1245
1246 return return_val;
1247}

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:197
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:55
Nektar::MultiRegions::GetPreconFactory
PreconFactory & GetPreconFactory()
Definition: Preconditioner.cpp:64

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 79 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

◆ m_InvRBlk

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

Definition at line 81 of file PreconditionerLowEnergy.h.

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

◆ m_RBlk

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

Definition at line 80 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 Jan 31 2024 09:36:05 for Nektar++ by doxygen 1.9.4