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
 

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

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

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

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

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:2266
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 1827 of file PreconditionerLowEnergy.cpp.

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

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

116{
117 std::shared_ptr<MultiRegions::ExpList> expList =
118 ((m_linsys.lock())->GetLocMat()).lock();
119 LocalRegions::Expansion3DSharedPtr locExpansion;
120 GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
121
122 LibUtilities::CommSharedPtr vComm = expList->GetComm()->GetSpaceComm();
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 (vComm->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 vComm->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, vComm, 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 vComm->AllReduce(maxEdgeDof, ReduceMax);
490 vComm->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 vComm->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, vComm, 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, vComm, 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:278
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_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 LocalRegions::ExpansionSharedPtr locExpansion;
84
85 locExpansion = expList->GetExp(0);
86
87 int nDim = locExpansion->GetShapeDimension();
88
89 ASSERTL0(nDim == 3, "Preconditioner type only valid in 3D");
90
91 // Set up block transformation matrix
92 SetupBlockTransformationMatrix();
93
94 // Sets up variable p mask
95 CreateVariablePMask();
96}

References ASSERTL0, CreateVariablePMask(), Nektar::MultiRegions::eIterativeStaticCond, Nektar::MultiRegions::ePETScStaticCond, 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: BasicUtils/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 May 15 2024 07:37:08 for Nektar++ by doxygen 1.9.4