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

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

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

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

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

1256{
1257 auto asmMap = m_locToGloMap.lock();
1258 unsigned int nLocBnd = asmMap->GetNumLocalBndCoeffs();
1259
1260 m_signChange = asmMap->GetSignChange();
1261
1262 // Construct a mask array
1263 m_variablePmask = Array<OneD, NekDouble>(nLocBnd, 1.0);
1264
1265 if (m_signChange)
1266 {
1267 unsigned int i;
1268 const Array<OneD, const NekDouble> &sign =
1269 asmMap->GetLocalToGlobalBndSign();
1270
1271 for (i = 0; i < nLocBnd; ++i)
1272 {
1273 m_variablePmask[i] = fabs(sign[i]);
1274 }
1275 }
1276}
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 2273 of file PreconditionerLowEnergy.cpp.

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

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

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

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

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

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

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

1627{
1628 std::shared_ptr<MultiRegions::ExpList> expList =
1629 ((m_linsys.lock())->GetLocMat()).lock();
1630 GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
1631
1632 LibUtilities::CommSharedPtr vComm = expList->GetComm()->GetSpaceComm();
1633
1634 LocalRegions::ExpansionSharedPtr locExp;
1635
1636 // face maps for pyramid and hybrid meshes - not need to return.
1637 map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> faceMapMaxR;
1638
1639 /* Determine the maximum expansion order for all elements */
1640 int nummodesmax = 0;
1641 Array<OneD, int> Shapes(LibUtilities::SIZE_ShapeType, 0);
1642
1643 for (int n = 0; n < expList->GetNumElmts(); ++n)
1644 {
1645 locExp = expList->GetExp(n);
1646
1647 nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(0));
1648 nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(1));
1649 nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(2));
1650
1651 Shapes[locExp->DetShapeType()] = 1;
1652 }
1653
1654 vComm->AllReduce(nummodesmax, ReduceMax);
1655 vComm->AllReduce(Shapes, ReduceMax);
1656
1657 if (Shapes[ePyramid] || Shapes[ePrism]) // if Pyramids or Prisms used then
1658 // need Tet and Hex expansion
1659 {
1660 Shapes[eTetrahedron] = 1;
1661 Shapes[eHexahedron] = 1;
1662 }
1663
1664 StdRegions::MatrixType PreconR;
1665 if (linSysKey.GetMatrixType() == StdRegions::eMass)
1666 {
1667 PreconR = StdRegions::ePreconRMass;
1668 }
1669 else
1670 {
1671 PreconR = StdRegions::ePreconR;
1672 }
1673
1674 Array<OneD, unsigned int> vmap;
1675 Array<OneD, Array<OneD, unsigned int>> emap;
1676 Array<OneD, Array<OneD, unsigned int>> fmap;
1677
1678 /*
1679 * Set up a transformation matrices for equal max order
1680 * polynomial meshes
1681 */
1682
1683 if (Shapes[eHexahedron])
1684 {
1685 SpatialDomains::HexGeomSharedPtr hexgeom = CreateRefHexGeom();
1686 // Bases for Hex element
1687 const BasisKey HexBa(eModified_A, nummodesmax,
1688 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1689 const BasisKey HexBb(eModified_A, nummodesmax,
1690 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1691 const BasisKey HexBc(eModified_A, nummodesmax,
1692 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1693
1694 // Create reference Hexahdedral expansion
1695 LocalRegions::HexExpSharedPtr HexExp;
1696
1697 HexExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
1698 HexBa, HexBb, HexBc, hexgeom);
1699
1700 maxElmt[eHexahedron] = HexExp;
1701
1702 // Hex:
1703 HexExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1704 vertMapMaxR[eHexahedron] = vmap;
1705 edgeMapMaxR[eHexahedron] = emap;
1706 faceMapMaxR[eHexahedron] = fmap;
1707
1708 // Get hexahedral transformation matrix
1709 LocalRegions::MatrixKey HexR(PreconR, eHexahedron, *HexExp,
1710 linSysKey.GetConstFactors());
1711 maxRmat[eHexahedron] = HexExp->GetLocMatrix(HexR);
1712 HexExp->DropLocMatrix(HexR); // Need to delete reference from manager
1713 }
1714
1715 if (Shapes[eTetrahedron])
1716 {
1717 SpatialDomains::TetGeomSharedPtr tetgeom = CreateRefTetGeom();
1718 // Bases for Tetrahedral element
1719 const BasisKey TetBa(eModified_A, nummodesmax,
1720 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1721 const BasisKey TetBb(eModified_B, nummodesmax,
1722 PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
1723 const BasisKey TetBc(eModified_C, nummodesmax,
1724 PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
1725
1726 // Create reference tetrahedral expansion
1727 LocalRegions::TetExpSharedPtr TetExp;
1728
1729 TetExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
1730 TetBa, TetBb, TetBc, tetgeom);
1731
1732 maxElmt[eTetrahedron] = TetExp;
1733
1734 TetExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1735 vertMapMaxR[eTetrahedron] = vmap;
1736 edgeMapMaxR[eTetrahedron] = emap;
1737 faceMapMaxR[eTetrahedron] = fmap;
1738
1739 // Get tetrahedral transformation matrix
1740 LocalRegions::MatrixKey TetR(PreconR, eTetrahedron, *TetExp,
1741 linSysKey.GetConstFactors());
1742 maxRmat[eTetrahedron] = TetExp->GetLocMatrix(TetR);
1743 TetExp->DropLocMatrix(TetR); // Need to delete reference from manager
1744
1745 if ((Shapes[ePyramid]) || (Shapes[eHexahedron]))
1746 {
1747 ReSetTetMaxRMat(nummodesmax, TetExp, maxRmat, vertMapMaxR,
1748 edgeMapMaxR, faceMapMaxR);
1749 }
1750 }
1751
1752 if (Shapes[ePyramid])
1753 {
1754 SpatialDomains::PyrGeomSharedPtr pyrgeom = CreateRefPyrGeom();
1755
1756 // Bases for Pyramid element
1757 const BasisKey PyrBa(eModified_A, nummodesmax,
1758 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1759 const BasisKey PyrBb(eModified_A, nummodesmax,
1760 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1761 const BasisKey PyrBc(eModifiedPyr_C, nummodesmax,
1762 PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
1763
1764 // Create reference pyramid expansion
1765 LocalRegions::PyrExpSharedPtr PyrExp;
1766
1767 PyrExp = MemoryManager<LocalRegions::PyrExp>::AllocateSharedPtr(
1768 PyrBa, PyrBb, PyrBc, pyrgeom);
1769
1770 maxElmt[ePyramid] = PyrExp;
1771
1772 // Pyramid:
1773 PyrExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1774 vertMapMaxR[ePyramid] = vmap;
1775 edgeMapMaxR[ePyramid] = emap;
1776 faceMapMaxR[ePyramid] = fmap;
1777
1778 // Set up Pyramid Transformation Matrix based on Tet
1779 // and Hex expansion
1780 SetUpPyrMaxRMat(nummodesmax, PyrExp, maxRmat, vertMapMaxR, edgeMapMaxR,
1781 faceMapMaxR);
1782 }
1783
1784 if (Shapes[ePrism])
1785 {
1786 SpatialDomains::PrismGeomSharedPtr prismgeom = CreateRefPrismGeom();
1787 // Bases for Prism element
1788 const BasisKey PrismBa(
1789 eModified_A, nummodesmax,
1790 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1791 const BasisKey PrismBb(
1792 eModified_A, nummodesmax,
1793 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1794 const BasisKey PrismBc(eModified_B, nummodesmax,
1795 PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
1796
1797 // Create reference prismatic expansion
1798 LocalRegions::PrismExpSharedPtr PrismExp;
1799
1800 PrismExp = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
1801 PrismBa, PrismBb, PrismBc, prismgeom);
1802 maxElmt[ePrism] = PrismExp;
1803
1804 // Prism:
1805 PrismExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1806 vertMapMaxR[ePrism] = vmap;
1807 edgeMapMaxR[ePrism] = emap;
1808
1809 faceMapMaxR[ePrism] = fmap;
1810
1811 if ((Shapes[ePyramid]) || (Shapes[eHexahedron]))
1812 {
1813 ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
1814 edgeMapMaxR, faceMapMaxR, false);
1815 }
1816 else
1817 {
1818 // Get prismatic transformation matrix
1819 LocalRegions::MatrixKey PrismR(PreconR, ePrism, *PrismExp,
1820 linSysKey.GetConstFactors());
1821 maxRmat[ePrism] = PrismExp->GetLocMatrix(PrismR);
1822 PrismExp->DropLocMatrix(
1823 PrismR); // Need to delete reference from manager
1824
1825 if (Shapes[eTetrahedron]) // reset triangular faces from Tet
1826 {
1827 ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
1828 edgeMapMaxR, faceMapMaxR, true);
1829 }
1830 }
1831 }
1832}
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefHexGeom
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1545
Nektar::MultiRegions::PreconditionerLowEnergy::ReSetPrismMaxRMat
void ReSetPrismMaxRMat(int nummodesmax, LocalRegions::PrismExpSharedPtr &PirsmExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR, bool UseTetOnly)
Definition: PreconditionerLowEnergy.cpp:2057
Nektar::MultiRegions::PreconditionerLowEnergy::ReSetTetMaxRMat
void ReSetTetMaxRMat(int nummodesmax, LocalRegions::TetExpSharedPtr &TetExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
Definition: PreconditionerLowEnergy.cpp:1999
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefTetGeom
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom(void)
Sets up the reference tretrahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1469
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPrismGeom
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1282
Nektar::MultiRegions::PreconditionerLowEnergy::SetUpPyrMaxRMat
void SetUpPyrMaxRMat(int nummodesmax, LocalRegions::PyrExpSharedPtr &PyrExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
Definition: PreconditionerLowEnergy.cpp:1834
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:1379
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 Gs::Free(tmp);
321
322 cnt = 0;
323 for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
324 {
325 it->second = (int)FacetLen[cnt];
326 }
327
328 for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
329 {
330 it->second = (int)FacetLen[cnt];
331 }
332 }
333
334 // Loop over all the elements in the domain and compute total edge
335 // DOF and set up unique ordering.
336 map<int, int> nblks;
337 int matrixlocation = 0;
338
339 // First do periodic edges
340 for (auto &pIt : periodicEdges)
341 {
342 meshEdgeId = pIt.first;
343
344 if (edgeDirMap.count(meshEdgeId) == 0)
345 {
346 dof = EdgeSize[meshEdgeId];
347 if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
348 {
349 bool SetUpNewEdge = true;
350
351 for (i = 0; i < pIt.second.size(); ++i)
352 {
353 if (!pIt.second[i].isLocal)
354 {
355 continue;
356 }
357
358 int meshEdgeId2 = pIt.second[i].id;
359 if (edgeDirMap.count(meshEdgeId2) == 0)
360 {
361 if (uniqueEdgeMap.count(meshEdgeId2) != 0)
362 {
363 // set unique map to same location
364 uniqueEdgeMap[meshEdgeId] =
365 uniqueEdgeMap[meshEdgeId2];
366 SetUpNewEdge = false;
367 }
368 }
369 else
370 {
371 edgeDirMap.insert(meshEdgeId);
372 SetUpNewEdge = false;
373 }
374 }
375
376 if (SetUpNewEdge)
377 {
378 uniqueEdgeMap[meshEdgeId] = matrixlocation;
379 globaloffset[matrixlocation] += ntotalentries;
380 modeoffset[matrixlocation] = dof * dof;
381 ntotalentries += dof * dof;
382 nblks[matrixlocation++] = dof;
383 }
384 }
385 }
386 }
387
388 for (cnt = n = 0; n < n_exp; ++n)
389 {
390 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
391
392 for (j = 0; j < locExpansion->GetNedges(); ++j)
393 {
394 meshEdgeId = locExpansion->GetGeom()->GetEid(j);
395 dof = EdgeSize[meshEdgeId];
396 maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
397
398 if (edgeDirMap.count(meshEdgeId) == 0)
399 {
400 if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
401
402 {
403 uniqueEdgeMap[meshEdgeId] = matrixlocation;
404
405 globaloffset[matrixlocation] += ntotalentries;
406 modeoffset[matrixlocation] = dof * dof;
407 ntotalentries += dof * dof;
408 nblks[matrixlocation++] = dof;
409 }
410 nlocalNonDirEdges += dof * dof;
411 }
412 }
413 }
414
415 // Loop over all the elements in the domain and compute max face
416 // DOF. Reduce across all processes to get universal maximum.
417 // - Periodic faces
418 for (auto &pIt : periodicFaces)
419 {
420 meshFaceId = pIt.first;
421
422 if (faceDirMap.count(meshFaceId) == 0)
423 {
424 dof = FaceSize[meshFaceId];
425
426 if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
427 {
428 bool SetUpNewFace = true;
429
430 if (pIt.second[0].isLocal)
431 {
432 int meshFaceId2 = pIt.second[0].id;
433
434 if (faceDirMap.count(meshFaceId2) == 0)
435 {
436 if (uniqueFaceMap.count(meshFaceId2) != 0)
437 {
438 // set unique map to same location
439 uniqueFaceMap[meshFaceId] =
440 uniqueFaceMap[meshFaceId2];
441 SetUpNewFace = false;
442 }
443 }
444 else // set face to be a Dirichlet face
445 {
446 faceDirMap.insert(meshFaceId);
447 SetUpNewFace = false;
448 }
449 }
450
451 if (SetUpNewFace)
452 {
453 uniqueFaceMap[meshFaceId] = matrixlocation;
454
455 modeoffset[matrixlocation] = dof * dof;
456 globaloffset[matrixlocation] += ntotalentries;
457 ntotalentries += dof * dof;
458 nblks[matrixlocation++] = dof;
459 }
460 }
461 }
462 }
463
464 for (cnt = n = 0; n < n_exp; ++n)
465 {
466 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
467
468 for (j = 0; j < locExpansion->GetNtraces(); ++j)
469 {
470 meshFaceId = locExpansion->GetGeom()->GetFid(j);
471
472 dof = FaceSize[meshFaceId];
473 maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
474
475 if (faceDirMap.count(meshFaceId) == 0)
476 {
477 if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
478 {
479 uniqueFaceMap[meshFaceId] = matrixlocation;
480 modeoffset[matrixlocation] = dof * dof;
481 globaloffset[matrixlocation] += ntotalentries;
482 ntotalentries += dof * dof;
483 nblks[matrixlocation++] = dof;
484 }
485 nlocalNonDirFaces += dof * dof;
486 }
487 }
488 }
489
490 vComm->AllReduce(maxEdgeDof, ReduceMax);
491 vComm->AllReduce(maxFaceDof, ReduceMax);
492
493 // Allocate arrays for block to universal map (number of expansions * p^2)
494 Array<OneD, long> BlockToUniversalMap(ntotalentries, -1);
495 Array<OneD, int> localToGlobalMatrixMap(
496 nlocalNonDirEdges + nlocalNonDirFaces, -1);
497
498 // Allocate arrays to store matrices (number of expansions * p^2)
499 Array<OneD, NekDouble> BlockArray(nlocalNonDirEdges + nlocalNonDirFaces,
500 0.0);
501
502 int matrixoffset = 0;
503 int vGlobal;
504 int uniEdgeOffset = 0;
505
506 // Need to obtain a fixed offset for the universal number
507 // of the faces which come after the edge numbering
508 for (n = 0; n < n_exp; ++n)
509 {
510 for (j = 0; j < locExpansion->GetNedges(); ++j)
511 {
512 meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
513 ->GetGeom3D()
514 ->GetEid(j);
515
516 uniEdgeOffset = max(meshEdgeId, uniEdgeOffset);
517 }
518 }
519 uniEdgeOffset++;
520
521 vComm->AllReduce(uniEdgeOffset, ReduceMax);
522 uniEdgeOffset = uniEdgeOffset * maxEdgeDof * maxEdgeDof;
523
524 for (n = 0; n < n_exp; ++n)
525 {
526 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
527
528 // loop over the edges of the expansion
529 for (j = 0; j < locExpansion->GetNedges(); ++j)
530 {
531 // get mesh edge id
532 meshEdgeId = locExpansion->GetGeom()->GetEid(j);
533
534 nedgemodes = EdgeSize[meshEdgeId];
535
536 if (edgeDirMap.count(meshEdgeId) == 0 && nedgemodes > 0)
537 {
538 // Determine the Global edge offset
539 int edgeOffset = globaloffset[uniqueEdgeMap[meshEdgeId]];
540
541 // Determine a universal map offset
542 int uniOffset = meshEdgeId;
543 auto pIt = periodicEdges.find(meshEdgeId);
544 if (pIt != periodicEdges.end())
545 {
546 for (int l = 0; l < pIt->second.size(); ++l)
547 {
548 uniOffset = min(uniOffset, pIt->second[l].id);
549 }
550 }
551 uniOffset = uniOffset * maxEdgeDof * maxEdgeDof;
552
553 for (k = 0; k < nedgemodes * nedgemodes; ++k)
554 {
555 vGlobal = edgeOffset + k;
556 localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
557 BlockToUniversalMap[vGlobal] = uniOffset + k + 1;
558 }
559 matrixoffset += nedgemodes * nedgemodes;
560 }
561 }
562
563 Array<OneD, unsigned int> faceInteriorMap;
564 Array<OneD, int> faceInteriorSign;
565 // loop over the faces of the expansion
566 for (j = 0; j < locExpansion->GetNtraces(); ++j)
567 {
568 // get mesh face id
569 meshFaceId = locExpansion->GetGeom()->GetFid(j);
570
571 nfacemodes = FaceSize[meshFaceId];
572
573 // Check if face has dirichlet values
574 if (faceDirMap.count(meshFaceId) == 0 && nfacemodes > 0)
575 {
576 // Determine the Global edge offset
577 int faceOffset = globaloffset[uniqueFaceMap[meshFaceId]];
578 // Determine a universal map offset
579 int uniOffset = meshFaceId;
580 // use minimum face edge when periodic
581 auto pIt = periodicFaces.find(meshFaceId);
582 if (pIt != periodicFaces.end())
583 {
584 uniOffset = min(uniOffset, pIt->second[0].id);
585 }
586 uniOffset = uniOffset * maxFaceDof * maxFaceDof;
587
588 for (k = 0; k < nfacemodes * nfacemodes; ++k)
589 {
590 vGlobal = faceOffset + k;
591
592 localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
593
594 BlockToUniversalMap[vGlobal] =
595 uniOffset + uniEdgeOffset + k + 1;
596 }
597 matrixoffset += nfacemodes * nfacemodes;
598 }
599 }
600 }
601
602 matrixoffset = 0;
603
604 map<int, int>::iterator it;
605 Array<OneD, unsigned int> n_blks(nblks.size() + 1);
606 n_blks[0] = nNonDirVerts;
607 for (i = 1, it = nblks.begin(); it != nblks.end(); ++it)
608 {
609 n_blks[i++] = it->second;
610 }
611
612 m_BlkMat = MemoryManager<DNekBlkMat>::AllocateSharedPtr(n_blks, n_blks,
613 blkmatStorage);
614
615 // Here we loop over the expansion and build the block low energy
616 // preconditioner as well as the block versions of the transformation
617 // matrices.
618 for (cnt = n = 0; n < n_exp; ++n)
619 {
620 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
621 nCoeffs = locExpansion->NumBndryCoeffs();
622
623 // Get correct transformation matrix for element type
624 DNekMat &R = (*m_RBlk->GetBlock(n, n));
625
626 pRSRT = MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs,
627 zero, storage);
628 RSRT = (*pRSRT);
629
630 nVerts = locExpansion->GetGeom()->GetNumVerts();
631 nEdges = locExpansion->GetGeom()->GetNumEdges();
632 nFaces = locExpansion->GetGeom()->GetNumFaces();
633
634 // Get statically condensed matrix
635 loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
636
637 // Extract boundary block (elemental S1)
638 bnd_mat = loc_mat->GetBlock(0, 0);
639
640 // offset by number of rows
641 offset = bnd_mat->GetRows();
642
643 DNekScalMat &S = (*bnd_mat);
644
645 DNekMat Sloc(nCoeffs, nCoeffs);
646
647 // For variable p we need to just use the active modes
648 NekDouble val;
649
650 for (int i = 0; i < nCoeffs; ++i)
651 {
652 for (int j = 0; j < nCoeffs; ++j)
653 {
654 val = S(i, j) * m_variablePmask[cnt + i] *
655 m_variablePmask[cnt + j];
656 Sloc.SetValue(i, j, val);
657 }
658 }
659
660 // Calculate R*S*trans(R)
661 RSRT = R * Sloc * Transpose(R);
662
663 // loop over vertices of the element and return the vertex map
664 // for each vertex
665 for (v = 0; v < nVerts; ++v)
666 {
667 vMap1 = locExpansion->GetVertexMap(v);
668
669 // Get vertex map
670 globalrow = asmMap->GetLocalToGlobalBndMap(cnt + vMap1) - nDirBnd;
671
672 if (globalrow >= 0)
673 {
674 for (m = 0; m < nVerts; ++m)
675 {
676 vMap2 = locExpansion->GetVertexMap(m);
677
678 // global matrix location (without offset due to
679 // dirichlet values)
680 globalcol =
681 asmMap->GetLocalToGlobalBndMap(cnt + vMap2) - nDirBnd;
682
683 // offset for dirichlet conditions
684 if (globalcol == globalrow)
685 {
686 // modal connectivity between elements
687 sign1 = asmMap->GetLocalToGlobalBndSign(cnt + vMap1);
688 sign2 = asmMap->GetLocalToGlobalBndSign(cnt + vMap2);
689
690 vertArray[globalrow] +=
691 sign1 * sign2 * RSRT(vMap1, vMap2);
692
693 meshVertId = locExpansion->GetGeom()->GetVid(v);
694
695 auto pIt = periodicVerts.find(meshVertId);
696 if (pIt != periodicVerts.end())
697 {
698 for (k = 0; k < pIt->second.size(); ++k)
699 {
700 meshVertId = min(meshVertId, pIt->second[k].id);
701 }
702 }
703
704 VertBlockToUniversalMap[globalrow] = meshVertId + 1;
705 }
706 }
707 }
708 }
709
710 // loop over edges of the element and return the edge map
711 for (eid = 0; eid < nEdges; ++eid)
712 {
713 meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
714 ->GetGeom3D()
715 ->GetEid(eid);
716
717 nedgemodes = EdgeSize[meshEdgeId];
718 if (nedgemodes)
719 {
720 nedgemodesloc = locExpansion->GetEdgeNcoeffs(eid) - 2;
721 DNekMatSharedPtr m_locMat =
722 MemoryManager<DNekMat>::AllocateSharedPtr(
723 nedgemodes, nedgemodes, zero, storage);
724
725 Array<OneD, unsigned int> edgemodearray =
726 locExpansion->GetEdgeInverseBoundaryMap(eid);
727
728 if (edgeDirMap.count(meshEdgeId) == 0)
729 {
730 for (v = 0; v < nedgemodesloc; ++v)
731 {
732 eMap1 = edgemodearray[v];
733 sign1 = asmMap->GetLocalToGlobalBndSign(cnt + eMap1);
734
735 if (sign1 == 0)
736 {
737 continue;
738 }
739
740 for (m = 0; m < nedgemodesloc; ++m)
741 {
742 eMap2 = edgemodearray[m];
743
744 // modal connectivity between elements
745 sign2 =
746 asmMap->GetLocalToGlobalBndSign(cnt + eMap2);
747
748 NekDouble globalEdgeValue =
749 sign1 * sign2 * RSRT(eMap1, eMap2);
750
751 if (sign2 != 0)
752 {
753 // if(eMap1 == eMap2)
754 BlockArray[matrixoffset + v * nedgemodes + m] =
755 globalEdgeValue;
756 }
757 }
758 }
759 matrixoffset += nedgemodes * nedgemodes;
760 }
761 }
762 }
763
764 // loop over faces of the element and return the face map
765 for (fid = 0; fid < nFaces; ++fid)
766 {
767 meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
768 ->GetGeom3D()
769 ->GetFid(fid);
770
771 nfacemodes = FaceSize[meshFaceId];
772 if (nfacemodes > 0)
773 {
774 DNekMatSharedPtr m_locMat =
775 MemoryManager<DNekMat>::AllocateSharedPtr(
776 nfacemodes, nfacemodes, zero, storage);
777
778 if (faceDirMap.count(meshFaceId) == 0)
779 {
780 Array<OneD, unsigned int> facemodearray;
781 StdRegions::Orientation faceOrient =
782 locExpansion->GetTraceOrient(fid);
783
784 auto pIt = periodicFaces.find(meshFaceId);
785 if (pIt != periodicFaces.end())
786 {
787 if (meshFaceId == min(meshFaceId, pIt->second[0].id))
788 {
789 faceOrient = DeterminePeriodicFaceOrient(
790 faceOrient, pIt->second[0].orient);
791 }
792 }
793
794 facemodearray = locExpansion->GetTraceInverseBoundaryMap(
795 fid, faceOrient, FaceModes[meshFaceId].first,
796 FaceModes[meshFaceId].second);
797
798 for (v = 0; v < nfacemodes; ++v)
799 {
800 fMap1 = facemodearray[v];
801
802 sign1 = asmMap->GetLocalToGlobalBndSign(cnt + fMap1);
803
804 ASSERTL1(sign1 != 0,
805 "Something is wrong since we "
806 "shoudl only be extracting modes for "
807 "lowest order expansion");
808
809 for (m = 0; m < nfacemodes; ++m)
810 {
811 fMap2 = facemodearray[m];
812
813 // modal connectivity between elements
814 sign2 =
815 asmMap->GetLocalToGlobalBndSign(cnt + fMap2);
816
817 ASSERTL1(sign2 != 0,
818 "Something is wrong since "
819 "we shoudl only be extracting modes for "
820 "lowest order expansion");
821
822 // Get the face-face value from the
823 // low energy matrix (S2)
824 NekDouble globalFaceValue =
825 sign1 * sign2 * RSRT(fMap1, fMap2);
826
827 // local face value to global face value
828 // if(fMap1 == fMap2)
829 BlockArray[matrixoffset + v * nfacemodes + m] =
830 globalFaceValue;
831 }
832 }
833 matrixoffset += nfacemodes * nfacemodes;
834 }
835 }
836 }
837
838 // offset for the expansion
839 cnt += nCoeffs;
840 }
841
842 // Exchange vertex data over different processes
843 Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, vComm, verbose);
844 Gs::Gather(vertArray, Gs::gs_add, tmp);
845 Gs::Free(tmp);
846
847 Array<OneD, NekDouble> GlobalBlock(ntotalentries, 0.0);
848 if (ntotalentries)
849 {
850 // Assemble edge matrices of each process
851 Vmath::Assmb(BlockArray.size(), BlockArray, localToGlobalMatrixMap,
852 GlobalBlock);
853 }
854
855 // Exchange edge & face data over different processes
856 Gs::gs_data *tmp1 = Gs::Init(BlockToUniversalMap, vComm, verbose);
857 Gs::Gather(GlobalBlock, Gs::gs_add, tmp1);
858 Gs::Free(tmp1);
859
860 // Populate vertex block
861 for (int i = 0; i < nNonDirVerts; ++i)
862 {
863 VertBlk->SetValue(i, i, 1.0 / vertArray[i]);
864 }
865
866 // Set the first block to be the diagonal of the vertex space
867 m_BlkMat->SetBlock(0, 0, VertBlk);
868
869 // Build the edge and face matrices from the vector
870 DNekMatSharedPtr gmat;
871
872 offset = 0;
873 // -1 since we ignore vert block
874 for (int loc = 0; loc < n_blks.size() - 1; ++loc)
875 {
876 nmodes = n_blks[1 + loc];
877 gmat = MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
878 storage);
879
880 for (v = 0; v < nmodes; ++v)
881 {
882 for (m = 0; m < nmodes; ++m)
883 {
884 NekDouble Value = GlobalBlock[offset + v * nmodes + m];
885 gmat->SetValue(v, m, Value);
886 }
887 }
888 m_BlkMat->SetBlock(1 + loc, 1 + loc, gmat);
889 offset += modeoffset[loc];
890 }
891
892 // invert blocks.
893 int totblks = m_BlkMat->GetNumberOfBlockRows();
894 for (i = 1; i < totblks; ++i)
895 {
896 unsigned int nmodes = m_BlkMat->GetNumberOfRowsInBlockRow(i);
897 if (nmodes)
898 {
899 DNekMatSharedPtr tmp_mat =
900 MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
901 storage);
902
903 tmp_mat = m_BlkMat->GetBlock(i, i);
904 tmp_mat->Invert();
905
906 m_BlkMat->SetBlock(i, i, tmp_mat);
907 }
908 }
909}
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
Gs::Free
static void Free(gs_data *pGsh)
Deallocates GSLib mapping data without finalising MPI.
Definition: GsLib.hpp:263
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::Free(), 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 915 of file PreconditionerLowEnergy.cpp.

919{
920 ASSERTL0(isLocal == false, "PreconditionerLowEnergy is only currently "
921 "set up for Global iteratives sovles");
922 int nDir = m_locToGloMap.lock()->GetNumGlobalDirBndCoeffs();
923 int nGlobal = m_locToGloMap.lock()->GetNumGlobalBndCoeffs();
924 int nNonDir = nGlobal - nDir;
925 DNekBlkMat &M = (*m_BlkMat);
926
927 NekVector<NekDouble> r(nNonDir, pInput, eWrapper);
928
929 NekVector<NekDouble> z(nNonDir, pOutput, eWrapper);
930
931 z = M * r;
932}
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 1130 of file PreconditionerLowEnergy.cpp.

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

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

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

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

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

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

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

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 Oct 30 2024 16:07:44 for Nektar++ by doxygen 1.9.4