#include <PreconditionerLowEnergy.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerLowEnergy:

Public Member Functions
	PreconditionerLowEnergy (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)

virtual	~PreconditionerLowEnergy ()

Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
	Preconditioner (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)

virtual	~Preconditioner ()

void	DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)

void	DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)

void	DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)

void	DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)

void	DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)

void	DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)

void	BuildPreconditioner ()

void	InitObject ()

Array< OneD, NekDouble >	AssembleStaticCondGlobalDiagonals ()
	Performs global assembly of diagonal entries to global Schur complement matrix. More...

const DNekScalBlkMatSharedPtr &	GetBlockTransformedSchurCompl () const

const DNekScalBlkMatSharedPtr &	GetBlockCMatrix () const

const DNekScalBlkMatSharedPtr &	GetBlockInvDMatrix () const

const DNekScalBlkMatSharedPtr &	GetBlockSchurCompl () const

const DNekScalBlkMatSharedPtr &	GetBlockTransformationMatrix () const

const DNekScalBlkMatSharedPtr &	GetBlockTransposedTransformationMatrix () const

DNekScalMatSharedPtr	TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat)

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

Static Public Attributes
static std::string	className
	Name of class. More...

Protected Member Functions
virtual void	v_InitObject () override

virtual void	v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override

virtual void	v_BuildPreconditioner () override
	Construct the low energy preconditioner from \(\mathbf{S}_{2}\). More...

virtual DNekScalMatSharedPtr	v_TransformedSchurCompl (int n, int offset, const std::shared_ptr< DNekScalMat > &loc_mat) override
	Set up the transformed block matrix system. More...

virtual void	v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut) override
	Transform the basis vector to low energy space. More...

virtual void	v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut) override
	transform the solution coeffiicents from low energy back to the original coefficient space. More...

virtual void	v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override
	Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis. More...

virtual void	v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) 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 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...

Protected Attributes
DNekBlkMatSharedPtr	m_BlkMat

DNekBlkMatSharedPtr	m_RBlk

DNekBlkMatSharedPtr	m_InvRBlk

Array< OneD, NekDouble >	m_variablePmask

bool	m_signChange

std::vector< std::pair< int, int > >	m_sameBlock

Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const std::weak_ptr< GlobalLinSys >	m_linsys

PreconditionerType	m_preconType

DNekMatSharedPtr	m_preconditioner

std::weak_ptr< AssemblyMap >	m_locToGloMap

LibUtilities::CommSharedPtr	m_comm

Private Types
typedef std::map< LibUtilities::ShapeType, DNekScalMatSharedPtr >	ShapeToDNekMap

typedef std::map< LibUtilities::ShapeType, LocalRegions::Expansion3DSharedPtr >	ShapeToExpMap

typedef std::map< LibUtilities::ShapeType, Array< OneD, unsigned int > >	ShapeToIntArrayMap

typedef std::map< LibUtilities::ShapeType, Array< OneD, Array< OneD, unsigned int > > >	ShapeToIntArrayArrayMap

Private Member Functions
void	SetupBlockTransformationMatrix (void)

void	SetUpReferenceElements (ShapeToDNekMap &maxRmat, ShapeToExpMap &maxElmt, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR)
	Loop expansion and determine different variants of the transformation matrix. More...

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

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

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

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

void	CreateVariablePMask (void)

SpatialDomains::TetGeomSharedPtr	CreateRefTetGeom (void)
	Sets up the reference tretrahedral element needed to construct a low energy basis. More...

SpatialDomains::PyrGeomSharedPtr	CreateRefPyrGeom (void)
	Sets up the reference prismatic element needed to construct a low energy basis mapping arrays. More...

SpatialDomains::PrismGeomSharedPtr	CreateRefPrismGeom (void)
	Sets up the reference prismatic element needed to construct a low energy basis. More...

SpatialDomains::HexGeomSharedPtr	CreateRefHexGeom (void)
	Sets up the reference hexahedral element needed to construct a low energy basis. More...

Detailed Description

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

Definition at line 53 of file PreconditionerLowEnergy.h.

Member Typedef Documentation

◆ ShapeToDNekMap

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

private

Definition at line 122 of file PreconditionerLowEnergy.h.

◆ ShapeToExpMap

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

private

Definition at line 125 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayArrayMap

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

private

Definition at line 130 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayMap

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

private

Definition at line 127 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

◆ PreconditionerLowEnergy()

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

Definition at line 66 of file PreconditionerLowEnergy.cpp.

     : Preconditioner(plinsys, pLocToGloMap)
 {
 }

◆ ~PreconditionerLowEnergy()

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

inlinevirtual

Definition at line 76 of file PreconditionerLowEnergy.h.

77 {

78 }

Member Function Documentation

◆ create()

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

inlinestatic

Creates an instance of this class.

Definition at line 57 of file PreconditionerLowEnergy.h.

     {
         PreconditionerSharedPtr p =
             MemoryManager<PreconditionerLowEnergy>::AllocateSharedPtr(
                 plinsys, pLocToGloMap);
         p->InitObject();
         return p;
     }

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.

 {
     ////////////////////////////////
     // Set up Hexahedron vertices //
     ////////////////////////////////
  
     const int three = 3;
  
     const int nVerts        = 8;
     const double point[][3] = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
                                {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
  
     // Populate the list of verts
     SpatialDomains::PointGeomSharedPtr verts[8];
     for (int i = 0; i < nVerts; ++i)
     {
         verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
             three, i, point[i][0], point[i][1], point[i][2]);
     }
  
     /////////////////////////////
     // Set up Hexahedron Edges //
     /////////////////////////////
  
     // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
     const int nEdges                  = 12;
     const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {2, 3}, {0, 3},
                                          {0, 4}, {1, 5}, {2, 6}, {3, 7},
                                          {4, 5}, {5, 6}, {6, 7}, {4, 7}};
  
     // Populate the list of edges
     SpatialDomains::SegGeomSharedPtr edges[nEdges];
     for (int i = 0; i < nEdges; ++i)
     {
         SpatialDomains::PointGeomSharedPtr vertsArray[2];
         for (int j = 0; j < 2; ++j)
         {
             vertsArray[j] = verts[vertexConnectivity[i][j]];
         }
         edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
             i, three, vertsArray);
     }
  
     /////////////////////////////
     // Set up Hexahedron faces //
     /////////////////////////////
  
     const int nFaces                = 6;
     const int edgeConnectivity[][4] = {{0, 1, 2, 3},  {0, 5, 8, 4},
                                        {1, 6, 9, 5},  {2, 7, 10, 6},
                                        {3, 7, 11, 4}, {8, 9, 10, 11}};
  
     // Populate the list of faces
     SpatialDomains::QuadGeomSharedPtr faces[nFaces];
     for (int i = 0; i < nFaces; ++i)
     {
         SpatialDomains::SegGeomSharedPtr edgeArray[4];
         for (int j = 0; j < 4; ++j)
         {
             edgeArray[j] = edges[edgeConnectivity[i][j]];
         }
         faces[i] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(
             i, edgeArray);
     }
  
     SpatialDomains::HexGeomSharedPtr geom =
         MemoryManager<SpatialDomains::HexGeom>::AllocateSharedPtr(0, faces);
  
     return geom;
 }

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.

 {
     //////////////////////////
     // Set up Prism element //
     //////////////////////////
  
     const int three         = 3;
     const int nVerts        = 6;
     const double point[][3] = {
         {-1, -1, 0},
         {1, -1, 0},
         {1, 1, 0},
         {-1, 1, 0},
         {0, -1, sqrt(double(3))},
         {0, 1, sqrt(double(3))},
     };
  
     // std::shared_ptr<SpatialDomains::PointGeom> verts[6];
     SpatialDomains::PointGeomSharedPtr verts[6];
     for (int i = 0; i < nVerts; ++i)
     {
         verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
             three, i, point[i][0], point[i][1], point[i][2]);
     }
     const int nEdges                  = 9;
     const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {3, 2}, {0, 3}, {0, 4},
                                          {1, 4}, {2, 5}, {3, 5}, {4, 5}};
  
     // Populate the list of edges
     SpatialDomains::SegGeomSharedPtr edges[nEdges];
     for (int i = 0; i < nEdges; ++i)
     {
         SpatialDomains::PointGeomSharedPtr vertsArray[2];
         for (int j = 0; j < 2; ++j)
         {
             vertsArray[j] = verts[vertexConnectivity[i][j]];
         }
         edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
             i, three, vertsArray);
     }
  
     ////////////////////////
     // Set up Prism faces //
     ////////////////////////
  
     const int nFaces = 5;
     // quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
     const int quadEdgeConnectivity[][4] = {
         {0, 1, 2, 3}, {1, 6, 8, 5}, {3, 7, 8, 4}};
     // QuadId ordered as 0, 1, 2, otherwise return false
     const int quadId[] = {0, -1, 1, -1, 2};
  
     // triangle-edge connectivity side-triface-1, side triface-3
     const int triEdgeConnectivity[][3] = {{0, 5, 4}, {2, 6, 7}};
     // TriId ordered as 0, 1, otherwise return false
     const int triId[] = {-1, 0, -1, 1, -1};
  
     // Populate the list of faces
     SpatialDomains::Geometry2DSharedPtr faces[nFaces];
     for (int f = 0; f < nFaces; ++f)
     {
         if (f == 1 || f == 3)
         {
             int i = triId[f];
             SpatialDomains::SegGeomSharedPtr edgeArray[3];
             for (int j = 0; j < 3; ++j)
             {
                 edgeArray[j] = edges[triEdgeConnectivity[i][j]];
             }
             faces[f] =
                 MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(
                     f, edgeArray);
         }
         else
         {
             int i = quadId[f];
             SpatialDomains::SegGeomSharedPtr edgeArray[4];
             for (int j = 0; j < 4; ++j)
             {
                 edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
             }
             faces[f] =
                 MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(
                     f, edgeArray);
         }
     }
  
     SpatialDomains::PrismGeomSharedPtr geom =
         MemoryManager<SpatialDomains::PrismGeom>::AllocateSharedPtr(0, faces);
  
     return geom;
 }

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.

 {
     //////////////////////////
     // Set up Pyramid element //
     //////////////////////////
  
     const int nVerts        = 5;
     const double point[][3] = {{-1, -1, 0},
                                {1, -1, 0},
                                {1, 1, 0},
                                {-1, 1, 0},
                                {0, 0, sqrt(double(2))}};
  
     // boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
     const int three = 3;
     SpatialDomains::PointGeomSharedPtr verts[5];
     for (int i = 0; i < nVerts; ++i)
     {
         verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
             three, i, point[i][0], point[i][1], point[i][2]);
     }
     const int nEdges                  = 8;
     const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0},
                                          {0, 4}, {1, 4}, {2, 4}, {3, 4}};
  
     // Populate the list of edges
     SpatialDomains::SegGeomSharedPtr edges[nEdges];
     for (int i = 0; i < nEdges; ++i)
     {
         SpatialDomains::PointGeomSharedPtr vertsArray[2];
         for (int j = 0; j < 2; ++j)
         {
             vertsArray[j] = verts[vertexConnectivity[i][j]];
         }
         edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
             i, three, vertsArray);
     }
  
     ////////////////////////
     // Set up Pyramid faces //
     ////////////////////////
  
     const int nFaces = 5;
     // quad-edge connectivity base-face0,
     const int quadEdgeConnectivity[][4] = {{0, 1, 2, 3}};
  
     // triangle-edge connectivity side-triface-1, side triface-2
     const int triEdgeConnectivity[][3] = {
         {0, 5, 4}, {1, 6, 5}, {2, 7, 6}, {3, 4, 7}};
  
     // Populate the list of faces
     SpatialDomains::Geometry2DSharedPtr faces[nFaces];
     for (int f = 0; f < nFaces; ++f)
     {
         if (f == 0)
         {
             SpatialDomains::SegGeomSharedPtr edgeArray[4];
             for (int j = 0; j < 4; ++j)
             {
                 edgeArray[j] = edges[quadEdgeConnectivity[f][j]];
             }
  
             faces[f] =
                 MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(
                     f, edgeArray);
         }
         else
         {
             int i = f - 1;
             SpatialDomains::SegGeomSharedPtr edgeArray[3];
             for (int j = 0; j < 3; ++j)
             {
                 edgeArray[j] = edges[triEdgeConnectivity[i][j]];
             }
             faces[f] =
                 MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(
                     f, edgeArray);
         }
     }
  
     SpatialDomains::PyrGeomSharedPtr geom =
         MemoryManager<SpatialDomains::PyrGeom>::AllocateSharedPtr(0, faces);
  
     return geom;
 }

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.

 {
     /////////////////////////////////
     // Set up Tetrahedron vertices //
     /////////////////////////////////
  
     int i, j;
     const int three         = 3;
     const int nVerts        = 4;
     const double point[][3] = {{-1, -1 / sqrt(double(3)), -1 / sqrt(double(6))},
                                {1, -1 / sqrt(double(3)), -1 / sqrt(double(6))},
                                {0, 2 / sqrt(double(3)), -1 / sqrt(double(6))},
                                {0, 0, 3 / sqrt(double(6))}};
  
     std::shared_ptr<SpatialDomains::PointGeom> verts[4];
     for (i = 0; i < nVerts; ++i)
     {
         verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
             three, i, point[i][0], point[i][1], point[i][2]);
     }
  
     //////////////////////////////
     // Set up Tetrahedron Edges //
     //////////////////////////////
  
     // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
     const int nEdges                  = 6;
     const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {0, 2},
                                          {0, 3}, {1, 3}, {2, 3}};
  
     // Populate the list of edges
     SpatialDomains::SegGeomSharedPtr edges[nEdges];
     for (i = 0; i < nEdges; ++i)
     {
         std::shared_ptr<SpatialDomains::PointGeom> vertsArray[2];
         for (j = 0; j < 2; ++j)
         {
             vertsArray[j] = verts[vertexConnectivity[i][j]];
         }
  
         edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
             i, three, vertsArray);
     }
  
     //////////////////////////////
     // Set up Tetrahedron faces //
     //////////////////////////////
  
     const int nFaces                = 4;
     const int edgeConnectivity[][3] = {
         {0, 1, 2}, {0, 4, 3}, {1, 5, 4}, {2, 5, 3}};
  
     // Populate the list of faces
     SpatialDomains::TriGeomSharedPtr faces[nFaces];
     for (i = 0; i < nFaces; ++i)
     {
         SpatialDomains::SegGeomSharedPtr edgeArray[3];
         for (j = 0; j < 3; ++j)
         {
             edgeArray[j] = edges[edgeConnectivity[i][j]];
         }
  
         faces[i] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(
             i, edgeArray);
     }
  
     SpatialDomains::TetGeomSharedPtr geom =
         MemoryManager<SpatialDomains::TetGeom>::AllocateSharedPtr(0, faces);
  
     return geom;
 }

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

Referenced by SetUpReferenceElements().

◆ CreateVariablePMask()

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

private

Create the inverse multiplicity map.

Definition at line 1252 of file PreconditionerLowEnergy.cpp.

 {
     unsigned int nLocBnd = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
     unsigned int i;
     auto asmMap = m_locToGloMap.lock();
  
     const Array<OneD, const NekDouble> &sign =
         asmMap->GetLocalToGlobalBndSign();
  
     m_signChange = asmMap->GetSignChange();
  
     // Construct a map of 1/multiplicity
     m_variablePmask = Array<OneD, NekDouble>(nLocBnd);
     for (i = 0; i < nLocBnd; ++i)
     {
         if (m_signChange)
         {
             m_variablePmask[i] = fabs(sign[i]);
         }
         else
         {
             m_variablePmask[i] = 1.0;
         }
     }
 }

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

 {
     NekDouble val;
     NekDouble zero = 0.0;
  
     int nRows = locExp->NumBndryCoeffs();
     DNekMatSharedPtr newmat =
         MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
  
     Array<OneD, unsigned int> vlocmap;
     Array<OneD, Array<OneD, unsigned int>> elocmap;
     Array<OneD, Array<OneD, unsigned int>> flocmap;
  
     locExp->GetInverseBoundaryMaps(vlocmap, elocmap, flocmap);
  
     // fill diagonal
     for (int i = 0; i < nRows; ++i)
     {
         val = 1.0;
         newmat->SetValue(i, i, val);
     }
  
     int nverts = locExp->GetNverts();
     int nedges = locExp->GetNedges();
     int nfaces = locExp->GetNtraces();
  
     // fill vertex to edge coupling
     for (int e = 0; e < nedges; ++e)
     {
         int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
  
         for (int v = 0; v < nverts; ++v)
         {
             for (int i = 0; i < nEdgeInteriorCoeffs; ++i)
             {
                 val = (*maxRmat)(vmap[v], emap[e][i]);
                 newmat->SetValue(vlocmap[v], elocmap[e][i], val);
             }
         }
     }
  
     for (int f = 0; f < nfaces; ++f)
     {
         // Get details to extrac this face from max reference matrix
         StdRegions::Orientation FwdOrient =
             StdRegions::eDir1FwdDir1_Dir2FwdDir2;
         int m0, m1; // Local face expansion orders.
  
         int nFaceInteriorCoeffs = locExp->GetTraceIntNcoeffs(f);
  
         locExp->GetTraceNumModes(f, m0, m1, FwdOrient);
  
         Array<OneD, unsigned int> fmapRmat =
             maxExp->GetTraceInverseBoundaryMap(f, FwdOrient, m0, m1);
  
         // fill in vertex to face coupling
         for (int v = 0; v < nverts; ++v)
         {
             for (int i = 0; i < nFaceInteriorCoeffs; ++i)
             {
                 val = (*maxRmat)(vmap[v], fmapRmat[i]);
                 newmat->SetValue(vlocmap[v], flocmap[f][i], val);
             }
         }
  
         // fill in edges to face coupling
         for (int e = 0; e < nedges; ++e)
         {
             int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
  
             for (int j = 0; j < nEdgeInteriorCoeffs; ++j)
             {
  
                 for (int i = 0; i < nFaceInteriorCoeffs; ++i)
                 {
                     val = (*maxRmat)(emap[e][j], fmapRmat[i]);
                     newmat->SetValue(elocmap[e][j], flocmap[f][i], val);
                 }
             }
         }
     }
  
     return newmat;
 }

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

 {
     int nRows = PrismExp->NumBndryCoeffs();
     NekDouble val;
     NekDouble zero = 0.0;
     DNekMatSharedPtr newmat =
         MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
  
     int nfacemodes;
  
     if (UseTetOnly)
     {
         // copy existing system
         for (int i = 0; i < nRows; ++i)
         {
             for (int j = 0; j < nRows; ++j)
             {
                 val = (*maxRmat[ePrism])(i, j);
                 newmat->SetValue(i, j, val);
             }
         }
  
         // Reset vertex to edge mapping from tet.
         const int VETetVert[][2]   = {{0, 0}, {1, 1}, {1, 1},
                                     {0, 0}, {3, 3}, {3, 3}};
         const int VETetEdge[][2]   = {{0, 3}, {0, 4}, {0, 4},
                                     {0, 3}, {3, 4}, {4, 3}};
         const int VEPrismEdge[][2] = {{0, 4}, {0, 5}, {2, 6},
                                       {2, 7}, {4, 5}, {6, 7}};
  
         // fill vertex to edge coupling
         for (int v = 0; v < 6; ++v)
         {
             for (int e = 0; e < 2; ++e)
             {
                 for (int i = 0; i < nummodesmax - 2; ++i)
                 {
                     // Note this is using wrong shape but gives
                     // answer that seems to have correct error!
                     val = (*maxRmat[eTetrahedron])(
                         vertMapMaxR[eTetrahedron][VETetVert[v][e]],
                         edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i]);
                     newmat->SetValue(vertMapMaxR[ePrism][v],
                                      edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
                                      val);
                 }
             }
         }
     }
     else
     {
  
         // set diagonal to 1
         for (int i = 0; i < nRows; ++i)
         {
             newmat->SetValue(i, i, 1.0);
         }
  
         // Set vertex to edge mapping from Hex.
  
         // The following lists specify the number of adjacent
         // edges to each vertex (nadj) then the Hex vert to
         // use for each prism ver in the vert-edge map (VEVert)
         // followed by the hex edge to use for each prism edge
         // in the vert-edge map (VEEdge)
         const int VEHexVert[][3]   = {{0, 0, 0}, {1, 1, 1}, {2, 2, 2},
                                     {3, 3, 3}, {4, 5, 5}, {6, 7, 7}};
         const int VEHexEdge[][3]   = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6},
                                     {2, 3, 7}, {4, 5, 9}, {6, 7, 11}};
         const int VEPrismEdge[][3] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6},
                                       {2, 3, 7}, {4, 5, 8}, {6, 7, 8}};
  
         // fill vertex to edge coupling
         for (int v = 0; v < 6; ++v)
         {
             for (int e = 0; e < 3; ++e)
             {
                 for (int i = 0; i < nummodesmax - 2; ++i)
                 {
                     // Note this is using wrong shape but gives
                     // answer that seems to have correct error!
                     val = (*maxRmat[eHexahedron])(
                         vertMapMaxR[eHexahedron][VEHexVert[v][e]],
                         edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
                     newmat->SetValue(vertMapMaxR[ePrism][v],
                                      edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
                                      val);
                 }
             }
         }
  
         // Setup vertex to face mapping from Hex
         const int VFHexVert[][2] = {{0, 0}, {1, 1}, {4, 5},
                                     {2, 2}, {3, 3}, {6, 7}};
         const int VFHexFace[][2] = {{0, 4}, {0, 2}, {4, 2},
                                     {0, 2}, {0, 4}, {2, 4}};
  
         const int VQFPrismVert[][2] = {{0, 0}, {1, 1}, {4, 4},
                                        {2, 2}, {3, 3}, {5, 5}};
         const int VQFPrismFace[][2] = {{0, 4}, {0, 2}, {4, 2},
                                        {0, 2}, {0, 4}, {2, 4}};
  
         nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
         // Replace two Quad faces  on every vertex
         for (int v = 0; v < 6; ++v)
         {
             for (int f = 0; f < 2; ++f)
             {
                 for (int i = 0; i < nfacemodes; ++i)
                 {
                     val = (*maxRmat[eHexahedron])(
                         vertMapMaxR[eHexahedron][VFHexVert[v][f]],
                         faceMapMaxR[eHexahedron][VFHexFace[v][f]][i]);
                     newmat->SetValue(vertMapMaxR[ePrism][VQFPrismVert[v][f]],
                                      faceMapMaxR[ePrism][VQFPrismFace[v][f]][i],
                                      val);
                 }
             }
         }
  
         // Mapping of Hex Edge-Face mappings to Prism Edge-Face Mappings
         const int nadjface[]     = {1, 2, 1, 2, 1, 1, 1, 1, 2};
         const int EFHexEdge[][2] = {{0, -1}, {1, 1},  {2, -1}, {3, 3}, {4, -1},
                                     {5, -1}, {6, -1}, {7, -1}, {9, 11}};
         const int EFHexFace[][2] = {{0, -1}, {0, 2},  {0, -1}, {0, 4}, {4, -1},
                                     {2, -1}, {2, -1}, {4, -1}, {2, 4}};
         const int EQFPrismEdge[][2] = {{0, -1}, {1, 1},  {2, -1},
                                        {3, 3},  {4, -1}, {5, -1},
                                        {6, -1}, {7, -1}, {8, 8}};
         const int EQFPrismFace[][2] = {{0, -1}, {0, 2},  {0, -1},
                                        {0, 4},  {4, -1}, {2, -1},
                                        {2, -1}, {4, -1}, {2, 4}};
  
         // all base edges are coupled to face 0
         nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
         for (int e = 0; e < 9; ++e)
         {
             for (int f = 0; f < nadjface[e]; ++f)
             {
                 for (int i = 0; i < nummodesmax - 2; ++i)
                 {
                     for (int j = 0; j < nfacemodes; ++j)
                     {
                         int edgemapid =
                             edgeMapMaxR[eHexahedron][EFHexEdge[e][f]][i];
                         int facemapid =
                             faceMapMaxR[eHexahedron][EFHexFace[e][f]][j];
  
                         val = (*maxRmat[eHexahedron])(edgemapid, facemapid);
  
                         int edgemapid1 =
                             edgeMapMaxR[ePrism][EQFPrismEdge[e][f]][i];
                         int facemapid1 =
                             faceMapMaxR[ePrism][EQFPrismFace[e][f]][j];
                         newmat->SetValue(edgemapid1, facemapid1, val);
                     }
                 }
             }
         }
     }
  
     const int VFTetVert[]    = {0, 1, 3, 1, 0, 3};
     const int VFTetFace[]    = {1, 1, 1, 1, 1, 1};
     const int VTFPrismVert[] = {0, 1, 4, 2, 3, 5};
     const int VTFPrismFace[] = {1, 1, 1, 3, 3, 3};
  
     //  Handle all triangular faces from tetrahedron
     nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
     for (int v = 0; v < 6; ++v)
     {
         for (int i = 0; i < nfacemodes; ++i)
         {
             val = (*maxRmat[eTetrahedron])(
                 vertMapMaxR[eTetrahedron][VFTetVert[v]],
                 faceMapMaxR[eTetrahedron][VFTetFace[v]][i]);
  
             newmat->SetValue(vertMapMaxR[ePrism][VTFPrismVert[v]],
                              faceMapMaxR[ePrism][VTFPrismFace[v]][i], val);
         }
     }
  
     // Mapping of Tet Edge-Face mappings to Prism Edge-Face Mappings
     const int EFTetEdge[]    = {0, 3, 4, 0, 4, 3};
     const int EFTetFace[]    = {1, 1, 1, 1, 1, 1};
     const int ETFPrismEdge[] = {0, 4, 5, 2, 6, 7};
     const int ETFPrismFace[] = {1, 1, 1, 3, 3, 3};
  
     // handle all edge to triangular faces from tetrahedron
     // (only 6 this time)
     nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
     for (int e = 0; e < 6; ++e)
     {
         for (int i = 0; i < nummodesmax - 2; ++i)
         {
             for (int j = 0; j < nfacemodes; ++j)
             {
                 int edgemapid = edgeMapMaxR[eTetrahedron][EFTetEdge[e]][i];
                 int facemapid = faceMapMaxR[eTetrahedron][EFTetFace[e]][j];
                 val           = (*maxRmat[eTetrahedron])(edgemapid, facemapid);
  
                 newmat->SetValue(edgeMapMaxR[ePrism][ETFPrismEdge[e]][i],
                                  faceMapMaxR[ePrism][ETFPrismFace[e]][j], val);
             }
         }
     }
  
     DNekScalMatSharedPtr PrismR =
         MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
     maxRmat[ePrism] = PrismR;
 }

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

 {
     boost::ignore_unused(faceMapMaxR);
  
     int nRows = TetExp->NumBndryCoeffs();
     NekDouble val;
     NekDouble zero = 0.0;
     DNekMatSharedPtr newmat =
         MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
  
     // copy existing system
     for (int i = 0; i < nRows; ++i)
     {
         for (int j = 0; j < nRows; ++j)
         {
             val = (*maxRmat[eTetrahedron])(i, j);
             newmat->SetValue(i, j, val);
         }
     }
  
     // The following lists specify the number of adjacent
     // edges to each vertex (nadj) then the Hex vert to
     // use for each pyramid ver in the vert-edge map (VEVert)
     // followed by the hex edge to use for each Tet edge
     // in the vert-edge map (VEEdge)
     const int VEHexVert[][4] = {{0, 0, 0}, {1, 1, 1}, {2, 2, 2}, {4, 5, 6}};
     const int VEHexEdge[][4] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6}, {4, 5, 6}};
     const int VETetEdge[][4] = {{0, 2, 3}, {0, 1, 4}, {1, 2, 5}, {3, 4, 5}};
  
     // fill vertex to edge coupling
     for (int v = 0; v < 4; ++v)
     {
         for (int e = 0; e < 3; ++e)
         {
             for (int i = 0; i < nummodesmax - 2; ++i)
             {
                 // Note this is using wrong shape but gives
                 // answer that seems to have correct error!
                 val = (*maxRmat[eHexahedron])(
                     vertMapMaxR[eHexahedron][VEHexVert[v][e]],
                     edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
                 newmat->SetValue(vertMapMaxR[eTetrahedron][v],
                                  edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i],
                                  val);
             }
         }
     }
  
     DNekScalMatSharedPtr TetR =
         MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
  
     maxRmat[eTetrahedron] = TetR;
 }

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

Referenced by SetUpReferenceElements().

◆ SetupBlockTransformationMatrix()

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

private

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

Definition at line 936 of file PreconditionerLowEnergy.cpp.

 {
     std::shared_ptr<MultiRegions::ExpList> expList =
         ((m_linsys.lock())->GetLocMat()).lock();
     LocalRegions::Expansion3DSharedPtr locExp;
     LocalRegions::Expansion3DSharedPtr locExpSav;
     map<int, int> EdgeSize;
  
     int n;
  
     std::map<ShapeType, DNekScalMatSharedPtr> maxRmat;
     map<ShapeType, LocalRegions::Expansion3DSharedPtr> maxElmt;
     map<ShapeType, Array<OneD, unsigned int>> vertMapMaxR;
     map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> edgeMapMaxR;
  
     // Sets up reference element and builds transformation matrix for
     // maximum polynomial order meshes
     SetUpReferenceElements(maxRmat, maxElmt, vertMapMaxR, edgeMapMaxR);
  
     const Array<OneD, const unsigned int> &nbdry_size =
         m_locToGloMap.lock()->GetNumLocalBndCoeffsPerPatch();
  
     int n_exp = expList->GetNumElmts();
  
     MatrixStorage blkmatStorage = eDIAGONAL;
  
     // Variants of R matrices required for low energy preconditioning
     m_RBlk = MemoryManager<DNekBlkMat>::AllocateSharedPtr(
         nbdry_size, nbdry_size, blkmatStorage);
     m_InvRBlk = MemoryManager<DNekBlkMat>::AllocateSharedPtr(
         nbdry_size, nbdry_size, blkmatStorage);
  
     DNekMatSharedPtr rmat, invrmat;
  
     int offset = 0;
  
     // Set up transformation matrices whilst checking to see if
     // consecutive matrices are the same and if so reuse the
     // matrices and store how many consecutive offsets there
     // are
     for (n = 0; n < n_exp; ++n)
     {
         locExp = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
  
         ShapeType eltype = locExp->DetShapeType();
  
         int nbndcoeffs = locExp->NumBndryCoeffs();
  
         if (m_sameBlock.size() == 0)
         {
             rmat = ExtractLocMat(locExp, maxRmat[eltype], maxElmt[eltype],
                                  vertMapMaxR[eltype], edgeMapMaxR[eltype]);
             // Block R matrix
             m_RBlk->SetBlock(n, n, rmat);
  
             invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
             invrmat->Invert();
  
             // Block inverse R matrix
             m_InvRBlk->SetBlock(n, n, invrmat);
  
             m_sameBlock.push_back(pair<int, int>(1, nbndcoeffs));
             locExpSav = locExp;
         }
         else
         {
             bool reuse = true;
  
             // check to see if same as previous matrix and
             // reuse if we can
             for (int i = 0; i < 3; ++i)
             {
                 if (locExpSav->GetBasis(i) != locExp->GetBasis(i))
                 {
                     reuse = false;
                     break;
                 }
             }
  
             if (reuse)
             {
                 m_RBlk->SetBlock(n, n, rmat);
                 m_InvRBlk->SetBlock(n, n, invrmat);
  
                 m_sameBlock[offset] =
                     (pair<int, int>(m_sameBlock[offset].first + 1, nbndcoeffs));
             }
             else
             {
                 rmat = ExtractLocMat(locExp, maxRmat[eltype], maxElmt[eltype],
                                      vertMapMaxR[eltype], edgeMapMaxR[eltype]);
  
                 // Block R matrix
                 m_RBlk->SetBlock(n, n, rmat);
  
                 invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
                 invrmat->Invert();
                 // Block inverse R matrix
                 m_InvRBlk->SetBlock(n, n, invrmat);
  
                 m_sameBlock.push_back(pair<int, int>(1, nbndcoeffs));
                 offset++;
                 locExpSav = locExp;
             }
         }
     }
 }

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

 {
     int nRows = PyrExp->NumBndryCoeffs();
     NekDouble val;
     NekDouble zero = 0.0;
     DNekMatSharedPtr newmat =
         MemoryManager<DNekMat>::AllocateSharedPtr(nRows, nRows, zero, eFULL);
  
     // set diagonal to 1
     for (int i = 0; i < nRows; ++i)
     {
         newmat->SetValue(i, i, 1.0);
     }
  
     // The following lists specify the number of adjacent
     // edges to each vertex (nadj) then the Hex vert to
     // use for each pyramid ver in the vert-edge map (VEVert)
     // followed by the hex edge to use for each pyramid edge
     // in the vert-edge map (VEEdge)
     const int nadjedge[]     = {3, 3, 3, 3, 4};
     const int VEHexVert[][4] = {{0, 0, 0, -1},
                                 {1, 1, 1, -1},
                                 {2, 2, 2, -1},
                                 {3, 3, 3, -1},
                                 {4, 5, 6, 7}};
     const int VEHexEdge[][4] = {{0, 3, 4, -1},
                                 {0, 1, 5, -1},
                                 {1, 2, 6, -1},
                                 {2, 3, 7, -1},
                                 {4, 5, 6, 7}};
     const int VEPyrEdge[][4] = {{0, 3, 4, -1},
                                 {0, 1, 5, -1},
                                 {1, 2, 6, -1},
                                 {2, 3, 7, -1},
                                 {4, 5, 6, 7}};
  
     // fill vertex to edge coupling
     for (int v = 0; v < 5; ++v)
     {
         for (int e = 0; e < nadjedge[v]; ++e)
         {
             for (int i = 0; i < nummodesmax - 2; ++i)
             {
                 // Note this is using wrong shape but gives
                 // answer that seems to have correct error!
                 val = (*maxRmat[eHexahedron])(
                     vertMapMaxR[eHexahedron][VEHexVert[v][e]],
                     edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
                 newmat->SetValue(vertMapMaxR[ePyramid][v],
                                  edgeMapMaxR[ePyramid][VEPyrEdge[v][e]][i],
                                  val);
             }
         }
     }
  
     int nfacemodes;
     nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
     // First four verties are all adjacent to base face
     for (int v = 0; v < 4; ++v)
     {
         for (int i = 0; i < nfacemodes; ++i)
         {
             val = (*maxRmat[eHexahedron])(vertMapMaxR[eHexahedron][v],
                                           faceMapMaxR[eHexahedron][0][i]);
             newmat->SetValue(vertMapMaxR[ePyramid][v],
                              faceMapMaxR[ePyramid][0][i], val);
         }
     }
  
     const int nadjface[]     = {2, 2, 2, 2, 4};
     const int VFTetVert[][4] = {{0, 0, -1, -1},
                                 {1, 1, -1, -1},
                                 {2, 2, -1, -1},
                                 {0, 2, -1, -1},
                                 {3, 3, 3, 3}};
     const int VFTetFace[][4] = {{1, 3, -1, -1},
                                 {1, 2, -1, -1},
                                 {2, 3, -1, -1},
                                 {1, 3, -1, -1},
                                 {1, 2, 1, 2}};
     const int VFPyrFace[][4] = {{1, 4, -1, -1},
                                 {1, 2, -1, -1},
                                 {2, 3, -1, -1},
                                 {3, 4, -1, -1},
                                 {1, 2, 3, 4}};
  
     // next handle all triangular faces from tetrahedron
     nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
     for (int v = 0; v < 5; ++v)
     {
         for (int f = 0; f < nadjface[v]; ++f)
         {
             for (int i = 0; i < nfacemodes; ++i)
             {
                 val = (*maxRmat[eTetrahedron])(
                     vertMapMaxR[eTetrahedron][VFTetVert[v][f]],
                     faceMapMaxR[eTetrahedron][VFTetFace[v][f]][i]);
                 newmat->SetValue(vertMapMaxR[ePyramid][v],
                                  faceMapMaxR[ePyramid][VFPyrFace[v][f]][i],
                                  val);
             }
         }
     }
  
     // Edge to face coupling
     // all base edges are coupled to face 0
     nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
     for (int e = 0; e < 4; ++e)
     {
         for (int i = 0; i < nummodesmax - 2; ++i)
         {
             for (int j = 0; j < nfacemodes; ++j)
             {
                 int edgemapid = edgeMapMaxR[eHexahedron][e][i];
                 int facemapid = faceMapMaxR[eHexahedron][0][j];
  
                 val = (*maxRmat[eHexahedron])(edgemapid, facemapid);
                 newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
                                  faceMapMaxR[ePyramid][0][j], val);
             }
         }
     }
  
     const int nadjface1[]    = {1, 1, 1, 1, 2, 2, 2, 2};
     const int EFTetEdge[][2] = {{0, -1}, {1, -1}, {0, -1}, {2, -1},
                                 {3, 3},  {4, 4},  {5, 5},  {3, 5}};
     const int EFTetFace[][2] = {{1, -1}, {2, -1}, {1, -1}, {3, -1},
                                 {1, 3},  {1, 2},  {2, 3},  {1, 3}};
     const int EFPyrFace[][2] = {{1, -1}, {2, -1}, {3, -1}, {4, -1},
                                 {1, 4},  {1, 2},  {2, 3},  {3, 4}};
  
     // next handle all triangular faces from tetrahedron
     nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
     for (int e = 0; e < 8; ++e)
     {
         for (int f = 0; f < nadjface1[e]; ++f)
         {
             for (int i = 0; i < nummodesmax - 2; ++i)
             {
                 for (int j = 0; j < nfacemodes; ++j)
                 {
                     int edgemapid =
                         edgeMapMaxR[eTetrahedron][EFTetEdge[e][f]][i];
                     int facemapid =
                         faceMapMaxR[eTetrahedron][EFTetFace[e][f]][j];
  
                     val = (*maxRmat[eTetrahedron])(edgemapid, facemapid);
                     newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
                                      faceMapMaxR[ePyramid][EFPyrFace[e][f]][j],
                                      val);
                 }
             }
         }
     }
  
     DNekScalMatSharedPtr PyrR;
     PyrR = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
     maxRmat[ePyramid] = PyrR;
 }

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.

 {
  
     std::shared_ptr<MultiRegions::ExpList> expList =
         ((m_linsys.lock())->GetLocMat()).lock();
     GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
  
     LocalRegions::ExpansionSharedPtr locExp;
  
     // face maps for pyramid and hybrid meshes - not need to return.
     map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> faceMapMaxR;
  
     /* Determine the maximum expansion order for all elements */
     int nummodesmax = 0;
     Array<OneD, int> Shapes(LibUtilities::SIZE_ShapeType, 0);
  
     for (int n = 0; n < expList->GetNumElmts(); ++n)
     {
         locExp = expList->GetExp(n);
  
         nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(0));
         nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(1));
         nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(2));
  
         Shapes[locExp->DetShapeType()] = 1;
     }
  
     m_comm->AllReduce(nummodesmax, ReduceMax);
     m_comm->AllReduce(Shapes, ReduceMax);
  
     if (Shapes[ePyramid] || Shapes[ePrism]) // if Pyramids or Prisms used then
                                             // need Tet and Hex expansion
     {
         Shapes[eTetrahedron] = 1;
         Shapes[eHexahedron]  = 1;
     }
  
     StdRegions::MatrixType PreconR;
     if (linSysKey.GetMatrixType() == StdRegions::eMass)
     {
         PreconR = StdRegions::ePreconRMass;
     }
     else
     {
         PreconR = StdRegions::ePreconR;
     }
  
     Array<OneD, unsigned int> vmap;
     Array<OneD, Array<OneD, unsigned int>> emap;
     Array<OneD, Array<OneD, unsigned int>> fmap;
  
     /*
      * Set up a transformation matrices for equal max order
      * polynomial meshes
      */
  
     if (Shapes[eHexahedron])
     {
         SpatialDomains::HexGeomSharedPtr hexgeom = CreateRefHexGeom();
         // Bases for Hex element
         const BasisKey HexBa(eModified_A, nummodesmax,
                              PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey HexBb(eModified_A, nummodesmax,
                              PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey HexBc(eModified_A, nummodesmax,
                              PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
  
         // Create reference Hexahdedral expansion
         LocalRegions::HexExpSharedPtr HexExp;
  
         HexExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
             HexBa, HexBb, HexBc, hexgeom);
  
         maxElmt[eHexahedron] = HexExp;
  
         // Hex:
         HexExp->GetInverseBoundaryMaps(vmap, emap, fmap);
         vertMapMaxR[eHexahedron] = vmap;
         edgeMapMaxR[eHexahedron] = emap;
         faceMapMaxR[eHexahedron] = fmap;
  
         // Get hexahedral transformation matrix
         LocalRegions::MatrixKey HexR(PreconR, eHexahedron, *HexExp,
                                      linSysKey.GetConstFactors());
         maxRmat[eHexahedron] = HexExp->GetLocMatrix(HexR);
     }
  
     if (Shapes[eTetrahedron])
     {
         SpatialDomains::TetGeomSharedPtr tetgeom = CreateRefTetGeom();
         // Bases for Tetrahedral element
         const BasisKey TetBa(eModified_A, nummodesmax,
                              PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey TetBb(eModified_B, nummodesmax,
                              PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
         const BasisKey TetBc(eModified_C, nummodesmax,
                              PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
  
         // Create reference tetrahedral expansion
         LocalRegions::TetExpSharedPtr TetExp;
  
         TetExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
             TetBa, TetBb, TetBc, tetgeom);
  
         maxElmt[eTetrahedron] = TetExp;
  
         TetExp->GetInverseBoundaryMaps(vmap, emap, fmap);
         vertMapMaxR[eTetrahedron] = vmap;
         edgeMapMaxR[eTetrahedron] = emap;
         faceMapMaxR[eTetrahedron] = fmap;
  
         // Get tetrahedral transformation matrix
         LocalRegions::MatrixKey TetR(PreconR, eTetrahedron, *TetExp,
                                      linSysKey.GetConstFactors());
         maxRmat[eTetrahedron] = TetExp->GetLocMatrix(TetR);
  
         if ((Shapes[ePyramid]) || (Shapes[eHexahedron]))
         {
             ReSetTetMaxRMat(nummodesmax, TetExp, maxRmat, vertMapMaxR,
                             edgeMapMaxR, faceMapMaxR);
         }
     }
  
     if (Shapes[ePyramid])
     {
         SpatialDomains::PyrGeomSharedPtr pyrgeom = CreateRefPyrGeom();
  
         // Bases for Pyramid element
         const BasisKey PyrBa(eModified_A, nummodesmax,
                              PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey PyrBb(eModified_A, nummodesmax,
                              PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey PyrBc(eModifiedPyr_C, nummodesmax,
                              PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
  
         // Create reference pyramid expansion
         LocalRegions::PyrExpSharedPtr PyrExp;
  
         PyrExp = MemoryManager<LocalRegions::PyrExp>::AllocateSharedPtr(
             PyrBa, PyrBb, PyrBc, pyrgeom);
  
         maxElmt[ePyramid] = PyrExp;
  
         // Pyramid:
         PyrExp->GetInverseBoundaryMaps(vmap, emap, fmap);
         vertMapMaxR[ePyramid] = vmap;
         edgeMapMaxR[ePyramid] = emap;
         faceMapMaxR[ePyramid] = fmap;
  
         // Set up Pyramid Transformation Matrix based on Tet
         // and Hex expansion
         SetUpPyrMaxRMat(nummodesmax, PyrExp, maxRmat, vertMapMaxR, edgeMapMaxR,
                         faceMapMaxR);
     }
  
     if (Shapes[ePrism])
     {
         SpatialDomains::PrismGeomSharedPtr prismgeom = CreateRefPrismGeom();
         // Bases for Prism element
         const BasisKey PrismBa(
             eModified_A, nummodesmax,
             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey PrismBb(
             eModified_A, nummodesmax,
             PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
         const BasisKey PrismBc(eModified_B, nummodesmax,
                                PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
  
         // Create reference prismatic expansion
         LocalRegions::PrismExpSharedPtr PrismExp;
  
         PrismExp = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
             PrismBa, PrismBb, PrismBc, prismgeom);
         maxElmt[ePrism] = PrismExp;
  
         // Prism:
         PrismExp->GetInverseBoundaryMaps(vmap, emap, fmap);
         vertMapMaxR[ePrism] = vmap;
         edgeMapMaxR[ePrism] = emap;
  
         faceMapMaxR[ePrism] = fmap;
  
         if ((Shapes[ePyramid]) || (Shapes[eHexahedron]))
         {
             ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
                               edgeMapMaxR, faceMapMaxR, false);
         }
         else
         {
  
             // Get prismatic transformation matrix
             LocalRegions::MatrixKey PrismR(PreconR, ePrism, *PrismExp,
                                            linSysKey.GetConstFactors());
             maxRmat[ePrism] = PrismExp->GetLocMatrix(PrismR);
  
             if (Shapes[eTetrahedron]) // reset triangular faces from Tet
             {
                 ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
                                   edgeMapMaxR, faceMapMaxR, true);
             }
         }
     }
 }

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

 {
     std::shared_ptr<MultiRegions::ExpList> expList =
         ((m_linsys.lock())->GetLocMat()).lock();
     LocalRegions::Expansion3DSharedPtr locExpansion;
     GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
  
     int i, j, k;
     int nVerts, nEdges, nFaces;
     int eid, fid, n, cnt, nmodes, nedgemodes, nfacemodes;
     int nedgemodesloc;
     NekDouble zero = 0.0;
  
     int vMap1, vMap2, sign1, sign2;
     int m, v, eMap1, eMap2, fMap1, fMap2;
     int offset, globalrow, globalcol, nCoeffs;
  
     // Periodic information
     PeriodicMap periodicVerts;
     PeriodicMap periodicEdges;
     PeriodicMap periodicFaces;
     expList->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
  
     // matrix storage
     MatrixStorage storage       = eFULL;
     MatrixStorage vertstorage   = eDIAGONAL;
     MatrixStorage blkmatStorage = eDIAGONAL;
  
     // local element static condensed matrices
     DNekScalBlkMatSharedPtr loc_mat;
     DNekScalMatSharedPtr bnd_mat;
  
     DNekMatSharedPtr pRSRT;
  
     DNekMat RS;
     DNekMat RSRT;
  
     auto asmMap = m_locToGloMap.lock();
  
     int nDirBnd      = asmMap->GetNumGlobalDirBndCoeffs();
     int nNonDirVerts = asmMap->GetNumNonDirVertexModes();
  
     // Vertex, edge and face preconditioner matrices
     DNekMatSharedPtr VertBlk = MemoryManager<DNekMat>::AllocateSharedPtr(
         nNonDirVerts, nNonDirVerts, zero, vertstorage);
  
     Array<OneD, NekDouble> vertArray(nNonDirVerts, 0.0);
     Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts, -1);
  
     // maps for different element types
     int n_exp          = expList->GetNumElmts();
     int nNonDirEdgeIDs = asmMap->GetNumNonDirEdges();
     int nNonDirFaceIDs = asmMap->GetNumNonDirFaces();
  
     set<int> edgeDirMap;
     set<int> faceDirMap;
     map<int, int> uniqueEdgeMap;
     map<int, int> uniqueFaceMap;
  
     // this should be of size total number of local edges + faces
     Array<OneD, int> modeoffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs, 0);
     Array<OneD, int> globaloffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs, 0);
  
     const Array<OneD, const ExpListSharedPtr> &bndCondExp =
         expList->GetBndCondExpansions();
     LocalRegions::Expansion2DSharedPtr bndCondFaceExp;
     const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
         &bndConditions = expList->GetBndConditions();
  
     int meshVertId;
     int meshEdgeId;
     int meshFaceId;
  
     const Array<OneD, const int> &extradiredges = asmMap->GetExtraDirEdges();
     for (i = 0; i < extradiredges.size(); ++i)
     {
         meshEdgeId = extradiredges[i];
         edgeDirMap.insert(meshEdgeId);
     }
  
     // Determine which boundary edges and faces have dirichlet values
     for (i = 0; i < bndCondExp.size(); i++)
     {
         cnt = 0;
         for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
         {
             bndCondFaceExp =
                 std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
                     bndCondExp[i]->GetExp(j));
             if (bndConditions[i]->GetBoundaryConditionType() ==
                 SpatialDomains::eDirichlet)
             {
                 for (k = 0; k < bndCondFaceExp->GetNtraces(); k++)
                 {
                     meshEdgeId = bndCondFaceExp->GetGeom()->GetEid(k);
                     if (edgeDirMap.count(meshEdgeId) == 0)
                     {
                         edgeDirMap.insert(meshEdgeId);
                     }
                 }
                 meshFaceId = bndCondFaceExp->GetGeom()->GetGlobalID();
                 faceDirMap.insert(meshFaceId);
             }
         }
     }
  
     int dof               = 0;
     int maxFaceDof        = 0;
     int maxEdgeDof        = 0;
     int nlocalNonDirEdges = 0;
     int nlocalNonDirFaces = 0;
     int ntotalentries     = 0;
  
     map<int, int> EdgeSize;
     map<int, int> FaceSize;
     map<int, pair<int, int>> FaceModes;
  
     /// -  Count  edges, face and add up min edges and min face sizes
     for (n = 0; n < n_exp; ++n)
     {
         locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
  
         nEdges = locExpansion->GetNedges();
         for (j = 0; j < nEdges; ++j)
         {
             int nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
             meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                              ->GetGeom3D()
                              ->GetEid(j);
             if (EdgeSize.count(meshEdgeId) == 0)
             {
                 EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
             }
             else
             {
                 EdgeSize[meshEdgeId] =
                     min(EdgeSize[meshEdgeId], nEdgeInteriorCoeffs);
             }
         }
  
         nFaces = locExpansion->GetNtraces();
         for (j = 0; j < nFaces; ++j)
         {
             int nFaceInteriorCoeffs = locExpansion->GetTraceIntNcoeffs(j);
             meshFaceId              = locExpansion->GetGeom3D()->GetFid(j);
  
             if (FaceSize.count(meshFaceId) == 0)
             {
                 FaceSize[meshFaceId] = nFaceInteriorCoeffs;
  
                 int m0, m1;
                 locExpansion->GetTraceNumModes(j, m0, m1,
                                                locExpansion->GetTraceOrient(j));
                 FaceModes[meshFaceId] = pair<int, int>(m0, m1);
             }
             else
             {
                 if (nFaceInteriorCoeffs < FaceSize[meshFaceId])
                 {
                     FaceSize[meshFaceId] = nFaceInteriorCoeffs;
                     int m0, m1;
                     locExpansion->GetTraceNumModes(
                         j, m0, m1, locExpansion->GetTraceOrient(j));
                     FaceModes[meshFaceId] = pair<int, int>(m0, m1);
                 }
             }
         }
     }
  
     bool verbose = expList->GetSession()->DefinesCmdLineArgument("verbose");
  
     // For parallel runs need to check have minimum of edges and faces over
     // partition boundaries
     if (m_comm->GetSize() > 1)
     {
         int EdgeSizeLen = EdgeSize.size();
         int FaceSizeLen = FaceSize.size();
         Array<OneD, long> FacetMap(EdgeSizeLen + FaceSizeLen, -1);
         Array<OneD, NekDouble> FacetLen(EdgeSizeLen + FaceSizeLen, -1);
  
         map<int, int>::iterator it;
  
         cnt       = 0;
         int maxid = 0;
         for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
         {
             FacetMap[cnt] = it->first;
             maxid         = max(it->first, maxid);
             FacetLen[cnt] = it->second;
         }
         maxid++;
  
         m_comm->AllReduce(maxid, ReduceMax);
  
         for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
         {
             FacetMap[cnt] = it->first + maxid;
             FacetLen[cnt] = it->second;
         }
  
         // Exchange vertex data over different processes
         Gs::gs_data *tmp = Gs::Init(FacetMap, m_comm, verbose);
         Gs::Gather(FacetLen, Gs::gs_min, tmp);
  
         cnt = 0;
         for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
         {
             it->second = (int)FacetLen[cnt];
         }
  
         for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
         {
             it->second = (int)FacetLen[cnt];
         }
     }
  
     // Loop over all the elements in the domain and compute total edge
     // DOF and set up unique ordering.
     map<int, int> nblks;
     int matrixlocation = 0;
  
     // First do periodic edges
     for (auto &pIt : periodicEdges)
     {
         meshEdgeId = pIt.first;
  
         if (edgeDirMap.count(meshEdgeId) == 0)
         {
             dof = EdgeSize[meshEdgeId];
             if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
             {
                 bool SetUpNewEdge = true;
  
                 for (i = 0; i < pIt.second.size(); ++i)
                 {
                     if (!pIt.second[i].isLocal)
                     {
                         continue;
                     }
  
                     int meshEdgeId2 = pIt.second[i].id;
                     if (edgeDirMap.count(meshEdgeId2) == 0)
                     {
                         if (uniqueEdgeMap.count(meshEdgeId2) != 0)
                         {
                             // set unique map to same location
                             uniqueEdgeMap[meshEdgeId] =
                                 uniqueEdgeMap[meshEdgeId2];
                             SetUpNewEdge = false;
                         }
                     }
                     else
                     {
                         edgeDirMap.insert(meshEdgeId);
                         SetUpNewEdge = false;
                     }
                 }
  
                 if (SetUpNewEdge)
                 {
                     uniqueEdgeMap[meshEdgeId] = matrixlocation;
                     globaloffset[matrixlocation] += ntotalentries;
                     modeoffset[matrixlocation] = dof * dof;
                     ntotalentries += dof * dof;
                     nblks[matrixlocation++] = dof;
                 }
             }
         }
     }
  
     for (cnt = n = 0; n < n_exp; ++n)
     {
         locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
  
         for (j = 0; j < locExpansion->GetNedges(); ++j)
         {
             meshEdgeId = locExpansion->GetGeom()->GetEid(j);
             dof        = EdgeSize[meshEdgeId];
             maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
  
             if (edgeDirMap.count(meshEdgeId) == 0)
             {
                 if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
  
                 {
                     uniqueEdgeMap[meshEdgeId] = matrixlocation;
  
                     globaloffset[matrixlocation] += ntotalentries;
                     modeoffset[matrixlocation] = dof * dof;
                     ntotalentries += dof * dof;
                     nblks[matrixlocation++] = dof;
                 }
                 nlocalNonDirEdges += dof * dof;
             }
         }
     }
  
     // Loop over all the elements in the domain and compute max face
     // DOF. Reduce across all processes to get universal maximum.
     // - Periodic faces
     for (auto &pIt : periodicFaces)
     {
         meshFaceId = pIt.first;
  
         if (faceDirMap.count(meshFaceId) == 0)
         {
             dof = FaceSize[meshFaceId];
  
             if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
             {
                 bool SetUpNewFace = true;
  
                 if (pIt.second[0].isLocal)
                 {
                     int meshFaceId2 = pIt.second[0].id;
  
                     if (faceDirMap.count(meshFaceId2) == 0)
                     {
                         if (uniqueFaceMap.count(meshFaceId2) != 0)
                         {
                             // set unique map to same location
                             uniqueFaceMap[meshFaceId] =
                                 uniqueFaceMap[meshFaceId2];
                             SetUpNewFace = false;
                         }
                     }
                     else // set face to be a Dirichlet face
                     {
                         faceDirMap.insert(meshFaceId);
                         SetUpNewFace = false;
                     }
                 }
  
                 if (SetUpNewFace)
                 {
                     uniqueFaceMap[meshFaceId] = matrixlocation;
  
                     modeoffset[matrixlocation] = dof * dof;
                     globaloffset[matrixlocation] += ntotalentries;
                     ntotalentries += dof * dof;
                     nblks[matrixlocation++] = dof;
                 }
             }
         }
     }
  
     for (cnt = n = 0; n < n_exp; ++n)
     {
         locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
  
         for (j = 0; j < locExpansion->GetNtraces(); ++j)
         {
             meshFaceId = locExpansion->GetGeom()->GetFid(j);
  
             dof        = FaceSize[meshFaceId];
             maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
  
             if (faceDirMap.count(meshFaceId) == 0)
             {
                 if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
                 {
                     uniqueFaceMap[meshFaceId]  = matrixlocation;
                     modeoffset[matrixlocation] = dof * dof;
                     globaloffset[matrixlocation] += ntotalentries;
                     ntotalentries += dof * dof;
                     nblks[matrixlocation++] = dof;
                 }
                 nlocalNonDirFaces += dof * dof;
             }
         }
     }
  
     m_comm->AllReduce(maxEdgeDof, ReduceMax);
     m_comm->AllReduce(maxFaceDof, ReduceMax);
  
     // Allocate arrays for block to universal map (number of expansions * p^2)
     Array<OneD, long> BlockToUniversalMap(ntotalentries, -1);
     Array<OneD, int> localToGlobalMatrixMap(
         nlocalNonDirEdges + nlocalNonDirFaces, -1);
  
     // Allocate arrays to store matrices (number of expansions * p^2)
     Array<OneD, NekDouble> BlockArray(nlocalNonDirEdges + nlocalNonDirFaces,
                                       0.0);
  
     int matrixoffset = 0;
     int vGlobal;
     int uniEdgeOffset = 0;
  
     // Need to obtain a fixed offset for the universal number
     // of the faces which come after the edge numbering
     for (n = 0; n < n_exp; ++n)
     {
         for (j = 0; j < locExpansion->GetNedges(); ++j)
         {
             meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                              ->GetGeom3D()
                              ->GetEid(j);
  
             uniEdgeOffset = max(meshEdgeId, uniEdgeOffset);
         }
     }
     uniEdgeOffset++;
  
     m_comm->AllReduce(uniEdgeOffset, ReduceMax);
     uniEdgeOffset = uniEdgeOffset * maxEdgeDof * maxEdgeDof;
  
     for (n = 0; n < n_exp; ++n)
     {
         locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
  
         // loop over the edges of the expansion
         for (j = 0; j < locExpansion->GetNedges(); ++j)
         {
             // get mesh edge id
             meshEdgeId = locExpansion->GetGeom()->GetEid(j);
  
             nedgemodes = EdgeSize[meshEdgeId];
  
             if (edgeDirMap.count(meshEdgeId) == 0 && nedgemodes > 0)
             {
                 // Determine the Global edge offset
                 int edgeOffset = globaloffset[uniqueEdgeMap[meshEdgeId]];
  
                 // Determine a universal map offset
                 int uniOffset = meshEdgeId;
                 auto pIt      = periodicEdges.find(meshEdgeId);
                 if (pIt != periodicEdges.end())
                 {
                     for (int l = 0; l < pIt->second.size(); ++l)
                     {
                         uniOffset = min(uniOffset, pIt->second[l].id);
                     }
                 }
                 uniOffset = uniOffset * maxEdgeDof * maxEdgeDof;
  
                 for (k = 0; k < nedgemodes * nedgemodes; ++k)
                 {
                     vGlobal                                  = edgeOffset + k;
                     localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
                     BlockToUniversalMap[vGlobal] = uniOffset + k + 1;
                 }
                 matrixoffset += nedgemodes * nedgemodes;
             }
         }
  
         Array<OneD, unsigned int> faceInteriorMap;
         Array<OneD, int> faceInteriorSign;
         // loop over the faces of the expansion
         for (j = 0; j < locExpansion->GetNtraces(); ++j)
         {
             // get mesh face id
             meshFaceId = locExpansion->GetGeom()->GetFid(j);
  
             nfacemodes = FaceSize[meshFaceId];
  
             // Check if face has dirichlet values
             if (faceDirMap.count(meshFaceId) == 0 && nfacemodes > 0)
             {
                 // Determine the Global edge offset
                 int faceOffset = globaloffset[uniqueFaceMap[meshFaceId]];
                 // Determine a universal map offset
                 int uniOffset = meshFaceId;
                 // use minimum face edge when periodic
                 auto pIt = periodicFaces.find(meshFaceId);
                 if (pIt != periodicFaces.end())
                 {
                     uniOffset = min(uniOffset, pIt->second[0].id);
                 }
                 uniOffset = uniOffset * maxFaceDof * maxFaceDof;
  
                 for (k = 0; k < nfacemodes * nfacemodes; ++k)
                 {
                     vGlobal = faceOffset + k;
  
                     localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
  
                     BlockToUniversalMap[vGlobal] =
                         uniOffset + uniEdgeOffset + k + 1;
                 }
                 matrixoffset += nfacemodes * nfacemodes;
             }
         }
     }
  
     matrixoffset = 0;
  
     map<int, int>::iterator it;
     Array<OneD, unsigned int> n_blks(nblks.size() + 1);
     n_blks[0] = nNonDirVerts;
     for (i = 1, it = nblks.begin(); it != nblks.end(); ++it)
     {
         n_blks[i++] = it->second;
     }
  
     m_BlkMat = MemoryManager<DNekBlkMat>::AllocateSharedPtr(n_blks, n_blks,
                                                             blkmatStorage);
  
     // Here we loop over the expansion and build the block low energy
     // preconditioner as well as the block versions of the transformation
     // matrices.
     for (cnt = n = 0; n < n_exp; ++n)
     {
         locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
         nCoeffs      = locExpansion->NumBndryCoeffs();
  
         // Get correct transformation matrix for element type
         DNekMat &R = (*m_RBlk->GetBlock(n, n));
  
         pRSRT = MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs,
                                                           zero, storage);
         RSRT  = (*pRSRT);
  
         nVerts = locExpansion->GetGeom()->GetNumVerts();
         nEdges = locExpansion->GetGeom()->GetNumEdges();
         nFaces = locExpansion->GetGeom()->GetNumFaces();
  
         // Get statically condensed matrix
         loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
  
         // Extract boundary block (elemental S1)
         bnd_mat = loc_mat->GetBlock(0, 0);
  
         // offset by number of rows
         offset = bnd_mat->GetRows();
  
         DNekScalMat &S = (*bnd_mat);
  
         DNekMat Sloc(nCoeffs, nCoeffs);
  
         // For variable p we need to just use the active modes
         NekDouble val;
  
         for (int i = 0; i < nCoeffs; ++i)
         {
             for (int j = 0; j < nCoeffs; ++j)
             {
                 val = S(i, j) * m_variablePmask[cnt + i] *
                       m_variablePmask[cnt + j];
                 Sloc.SetValue(i, j, val);
             }
         }
  
         // Calculate R*S*trans(R)
         RSRT = R * Sloc * Transpose(R);
  
         // loop over vertices of the element and return the vertex map
         // for each vertex
         for (v = 0; v < nVerts; ++v)
         {
             vMap1 = locExpansion->GetVertexMap(v);
  
             // Get vertex map
             globalrow = asmMap->GetLocalToGlobalBndMap(cnt + vMap1) - nDirBnd;
  
             if (globalrow >= 0)
             {
                 for (m = 0; m < nVerts; ++m)
                 {
                     vMap2 = locExpansion->GetVertexMap(m);
  
                     // global matrix location (without offset due to
                     // dirichlet values)
                     globalcol =
                         asmMap->GetLocalToGlobalBndMap(cnt + vMap2) - nDirBnd;
  
                     // offset for dirichlet conditions
                     if (globalcol == globalrow)
                     {
                         // modal connectivity between elements
                         sign1 = asmMap->GetLocalToGlobalBndSign(cnt + vMap1);
                         sign2 = asmMap->GetLocalToGlobalBndSign(cnt + vMap2);
  
                         vertArray[globalrow] +=
                             sign1 * sign2 * RSRT(vMap1, vMap2);
  
                         meshVertId = locExpansion->GetGeom()->GetVid(v);
  
                         auto pIt = periodicVerts.find(meshVertId);
                         if (pIt != periodicVerts.end())
                         {
                             for (k = 0; k < pIt->second.size(); ++k)
                             {
                                 meshVertId = min(meshVertId, pIt->second[k].id);
                             }
                         }
  
                         VertBlockToUniversalMap[globalrow] = meshVertId + 1;
                     }
                 }
             }
         }
  
         // loop over edges of the element and return the edge map
         for (eid = 0; eid < nEdges; ++eid)
         {
  
             meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                              ->GetGeom3D()
                              ->GetEid(eid);
  
             nedgemodes = EdgeSize[meshEdgeId];
             if (nedgemodes)
             {
                 nedgemodesloc = locExpansion->GetEdgeNcoeffs(eid) - 2;
                 DNekMatSharedPtr m_locMat =
                     MemoryManager<DNekMat>::AllocateSharedPtr(
                         nedgemodes, nedgemodes, zero, storage);
  
                 Array<OneD, unsigned int> edgemodearray =
                     locExpansion->GetEdgeInverseBoundaryMap(eid);
  
                 if (edgeDirMap.count(meshEdgeId) == 0)
                 {
                     for (v = 0; v < nedgemodesloc; ++v)
                     {
                         eMap1 = edgemodearray[v];
                         sign1 = asmMap->GetLocalToGlobalBndSign(cnt + eMap1);
  
                         if (sign1 == 0)
                         {
                             continue;
                         }
  
                         for (m = 0; m < nedgemodesloc; ++m)
                         {
                             eMap2 = edgemodearray[m];
  
                             // modal connectivity between elements
                             sign2 =
                                 asmMap->GetLocalToGlobalBndSign(cnt + eMap2);
  
                             NekDouble globalEdgeValue =
                                 sign1 * sign2 * RSRT(eMap1, eMap2);
  
                             if (sign2 != 0)
                             {
                                 // if(eMap1 == eMap2)
                                 BlockArray[matrixoffset + v * nedgemodes + m] =
                                     globalEdgeValue;
                             }
                         }
                     }
                     matrixoffset += nedgemodes * nedgemodes;
                 }
             }
         }
  
         // loop over faces of the element and return the face map
         for (fid = 0; fid < nFaces; ++fid)
         {
             meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
                              ->GetGeom3D()
                              ->GetFid(fid);
  
             nfacemodes = FaceSize[meshFaceId];
             if (nfacemodes > 0)
             {
                 DNekMatSharedPtr m_locMat =
                     MemoryManager<DNekMat>::AllocateSharedPtr(
                         nfacemodes, nfacemodes, zero, storage);
  
                 if (faceDirMap.count(meshFaceId) == 0)
                 {
                     Array<OneD, unsigned int> facemodearray;
                     StdRegions::Orientation faceOrient =
                         locExpansion->GetTraceOrient(fid);
  
                     auto pIt = periodicFaces.find(meshFaceId);
                     if (pIt != periodicFaces.end())
                     {
                         if (meshFaceId == min(meshFaceId, pIt->second[0].id))
                         {
                             faceOrient = DeterminePeriodicFaceOrient(
                                 faceOrient, pIt->second[0].orient);
                         }
                     }
  
                     facemodearray = locExpansion->GetTraceInverseBoundaryMap(
                         fid, faceOrient, FaceModes[meshFaceId].first,
                         FaceModes[meshFaceId].second);
  
                     for (v = 0; v < nfacemodes; ++v)
                     {
                         fMap1 = facemodearray[v];
  
                         sign1 = asmMap->GetLocalToGlobalBndSign(cnt + fMap1);
  
                         ASSERTL1(sign1 != 0,
                                  "Something is wrong since we "
                                  "shoudl only be extracting modes for "
                                  "lowest order expansion");
  
                         for (m = 0; m < nfacemodes; ++m)
                         {
                             fMap2 = facemodearray[m];
  
                             // modal connectivity between elements
                             sign2 =
                                 asmMap->GetLocalToGlobalBndSign(cnt + fMap2);
  
                             ASSERTL1(sign2 != 0,
                                      "Something is wrong since "
                                      "we shoudl only be extracting modes for "
                                      "lowest order expansion");
  
                             // Get the face-face value from the
                             // low energy matrix (S2)
                             NekDouble globalFaceValue =
                                 sign1 * sign2 * RSRT(fMap1, fMap2);
  
                             // local face value to global face value
                             // if(fMap1 == fMap2)
                             BlockArray[matrixoffset + v * nfacemodes + m] =
                                 globalFaceValue;
                         }
                     }
                     matrixoffset += nfacemodes * nfacemodes;
                 }
             }
         }
  
         // offset for the expansion
         cnt += nCoeffs;
     }
  
     if (nNonDirVerts != 0)
     {
         // Exchange vertex data over different processes
         Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm, verbose);
         Gs::Gather(vertArray, Gs::gs_add, tmp);
     }
  
     Array<OneD, NekDouble> GlobalBlock(ntotalentries, 0.0);
     if (ntotalentries)
     {
         // Assemble edge matrices of each process
         Vmath::Assmb(BlockArray.size(), BlockArray, localToGlobalMatrixMap,
                      GlobalBlock);
     }
  
     // Exchange edge & face data over different processes
     Gs::gs_data *tmp1 = Gs::Init(BlockToUniversalMap, m_comm, verbose);
     Gs::Gather(GlobalBlock, Gs::gs_add, tmp1);
  
     // Populate vertex block
     for (int i = 0; i < nNonDirVerts; ++i)
     {
         VertBlk->SetValue(i, i, 1.0 / vertArray[i]);
     }
  
     // Set the first block to be the diagonal of the vertex space
     m_BlkMat->SetBlock(0, 0, VertBlk);
  
     // Build the edge and face matrices from the vector
     DNekMatSharedPtr gmat;
  
     offset = 0;
     // -1 since we ignore vert block
     for (int loc = 0; loc < n_blks.size() - 1; ++loc)
     {
         nmodes = n_blks[1 + loc];
         gmat   = MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
                                                          storage);
  
         for (v = 0; v < nmodes; ++v)
         {
             for (m = 0; m < nmodes; ++m)
             {
                 NekDouble Value = GlobalBlock[offset + v * nmodes + m];
                 gmat->SetValue(v, m, Value);
             }
         }
         m_BlkMat->SetBlock(1 + loc, 1 + loc, gmat);
         offset += modeoffset[loc];
     }
  
     // invert blocks.
     int totblks = m_BlkMat->GetNumberOfBlockRows();
     for (i = 1; i < totblks; ++i)
     {
         unsigned int nmodes = m_BlkMat->GetNumberOfRowsInBlockRow(i);
         if (nmodes)
         {
             DNekMatSharedPtr tmp_mat =
                 MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
                                                           storage);
  
             tmp_mat = m_BlkMat->GetBlock(i, i);
             tmp_mat->Invert();
  
             m_BlkMat->SetBlock(i, i, tmp_mat);
         }
     }
 }

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

◆ v_DoPreconditioner()

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

overrideprotectedvirtual

Apply the low energy preconditioner during the conjugate gradient routine

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 917 of file PreconditionerLowEnergy.cpp.

 {
     int nDir      = m_locToGloMap.lock()->GetNumGlobalDirBndCoeffs();
     int nGlobal   = m_locToGloMap.lock()->GetNumGlobalBndCoeffs();
     int nNonDir   = nGlobal - nDir;
     DNekBlkMat &M = (*m_BlkMat);
  
     NekVector<NekDouble> r(nNonDir, pInput, eWrapper);
     NekVector<NekDouble> z(nNonDir, pOutput, eWrapper);
  
     z = M * r;
 }

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

◆ v_DoTransformBasisFromLowEnergy()

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

overrideprotectedvirtual

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

Note; not typically required

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1127 of file PreconditionerLowEnergy.cpp.

 {
     int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
  
     ASSERTL1(pInput.size() >= nLocBndDofs,
              "Input array is smaller than nLocBndDofs");
     ASSERTL1(pOutput.size() >= nLocBndDofs,
              "Output array is smaller than nLocBndDofs");
  
     // Block inverse transformation matrix
     DNekBlkMat &invR = *m_InvRBlk;
  
     Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
  
     // Multiply by the inverse transformation matrix
     int cnt  = 0;
     int cnt1 = 0;
     for (int i = 0; i < m_sameBlock.size(); ++i)
     {
         int nexp       = m_sameBlock[i].first;
         int nbndcoeffs = m_sameBlock[i].second;
         Blas::Dgemm('N', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                     &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                     pLocalIn.get() + cnt, nbndcoeffs, 0.0, pOutput.get() + cnt,
                     nbndcoeffs);
         cnt += nbndcoeffs * nexp;
         cnt1 += nexp;
     }
 }

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

◆ v_DoTransformBasisToLowEnergy()

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

overrideprotectedvirtual

Transform the basis vector to low energy space.

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1051 of file PreconditionerLowEnergy.cpp.

 {
     int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
  
     // Block transformation matrix
     DNekBlkMat &R = *m_RBlk;
  
     Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
  
     // Apply mask in case of variable P
     Vmath::Vmul(nLocBndDofs, pLocalIn, 1, m_variablePmask, 1, pLocalIn, 1);
  
     // Multiply by the block transformation matrix
     int cnt  = 0;
     int cnt1 = 0;
     for (int i = 0; i < m_sameBlock.size(); ++i)
     {
         int nexp       = m_sameBlock[i].first;
         int nbndcoeffs = m_sameBlock[i].second;
         Blas::Dgemm('N', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                     &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                     pLocalIn.get() + cnt, nbndcoeffs, 0.0, pInOut.get() + cnt,
                     nbndcoeffs);
         cnt += nbndcoeffs * nexp;
         cnt1 += nexp;
     }
 }

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

◆ v_DoTransformCoeffsFromLowEnergy()

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

overrideprotectedvirtual

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1089 of file PreconditionerLowEnergy.cpp.

 {
     int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
  
     ASSERTL1(pInOut.size() >= nLocBndDofs,
              "Output array is not greater than the nLocBndDofs");
  
     // Block transposed transformation matrix
     DNekBlkMat &R = *m_RBlk;
  
     Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
  
     // Multiply by the transpose of block transformation matrix
     int cnt  = 0;
     int cnt1 = 0;
     for (int i = 0; i < m_sameBlock.size(); ++i)
     {
         int nexp       = m_sameBlock[i].first;
         int nbndcoeffs = m_sameBlock[i].second;
         Blas::Dgemm('T', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                     &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                     pLocalIn.get() + cnt, nbndcoeffs, 0.0, pInOut.get() + cnt,
                     nbndcoeffs);
         cnt += nbndcoeffs * nexp;
         cnt1 += nexp;
     }
  
     Vmath::Vmul(nLocBndDofs, pInOut, 1, m_variablePmask, 1, pInOut, 1);
 }

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

◆ v_DoTransformCoeffsToLowEnergy()

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

overrideprotectedvirtual

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1165 of file PreconditionerLowEnergy.cpp.

 {
     int nLocBndDofs = m_locToGloMap.lock()->GetNumLocalBndCoeffs();
  
     ASSERTL1(pInput.size() >= nLocBndDofs,
              "Input array is less than nLocBndDofs");
     ASSERTL1(pOutput.size() >= nLocBndDofs,
              "Output array is less than nLocBndDofs");
  
     // Block inverse transformation matrix
     DNekBlkMat &invR = *m_InvRBlk;
  
     Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
  
     // Multiply by the transpose of block transformation matrix
     int cnt  = 0;
     int cnt1 = 0;
     for (int i = 0; i < m_sameBlock.size(); ++i)
     {
         int nexp       = m_sameBlock[i].first;
         int nbndcoeffs = m_sameBlock[i].second;
         Blas::Dgemm('T', 'N', nbndcoeffs, nexp, nbndcoeffs, 1.0,
                     &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
                     pLocalIn.get() + cnt, nbndcoeffs, 0.0, pOutput.get() + cnt,
                     nbndcoeffs);
         cnt += nbndcoeffs * nexp;
         cnt1 += nexp;
     }
 }

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

 {
     GlobalSysSolnType solvertype = m_locToGloMap.lock()->GetGlobalSysSolnType();
  
     ASSERTL0(solvertype == eIterativeStaticCond ||
                  solvertype == ePETScStaticCond,
              "Solver type not valid");
  
     std::shared_ptr<MultiRegions::ExpList> expList =
         ((m_linsys.lock())->GetLocMat()).lock();
  
     m_comm = expList->GetComm();
  
     LocalRegions::ExpansionSharedPtr locExpansion;
  
     locExpansion = expList->GetExp(0);
  
     int nDim = locExpansion->GetShapeDimension();
  
     ASSERTL0(nDim == 3, "Preconditioner type only valid in 3D");
  
     // Set up block transformation matrix
     SetupBlockTransformationMatrix();
  
     // Sets up variable p mask
     CreateVariablePMask();
 }

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

◆ v_TransformedSchurCompl()

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

overrideprotectedvirtual

Set up the transformed block matrix system.

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1203 of file PreconditionerLowEnergy.cpp.

 {
     std::shared_ptr<MultiRegions::ExpList> expList =
         ((m_linsys.lock())->GetLocMat()).lock();
  
     LocalRegions::ExpansionSharedPtr locExpansion;
     locExpansion      = expList->GetExp(n);
     unsigned int nbnd = locExpansion->NumBndryCoeffs();
  
     MatrixStorage storage = eFULL;
     DNekMatSharedPtr pS2 =
         MemoryManager<DNekMat>::AllocateSharedPtr(nbnd, nbnd, 0.0, storage);
  
     // transformation matrices
     DNekMat &R = (*m_RBlk->GetBlock(n, n));
  
     // Original Schur Complement
     DNekScalMat &S1 = (*loc_mat);
  
     DNekMat Sloc(nbnd, nbnd);
  
     // For variable p we need to just use the active modes
     NekDouble val;
  
     for (int i = 0; i < nbnd; ++i)
     {
         for (int j = 0; j < nbnd; ++j)
         {
             val = S1(i, j) * m_variablePmask[bndoffset + i] *
                   m_variablePmask[bndoffset + j];
             Sloc.SetValue(i, j, val);
         }
     }
  
     // create low energy matrix
     DNekMat &S2 = (*pS2);
  
     S2 = R * Sloc * Transpose(R);
  
     DNekScalMatSharedPtr return_val;
     return_val = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
  
     return return_val;
 }

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

Name of class.

Registers the class with the Factory.

Definition at line 69 of file PreconditionerLowEnergy.h.

◆ m_BlkMat

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_BlkMat

protected

Definition at line 81 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

◆ m_InvRBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_InvRBlk

protected

Definition at line 83 of file PreconditionerLowEnergy.h.

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

◆ m_RBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RBlk

protected

Definition at line 82 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 91 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 87 of file PreconditionerLowEnergy.h.

Referenced by CreateVariablePMask().

◆ m_variablePmask

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

protected

Definition at line 85 of file PreconditionerLowEnergy.h.

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

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Private Types

Private Member Functions

Detailed Description

Member Typedef Documentation

◆ ShapeToDNekMap

◆ ShapeToExpMap

◆ ShapeToIntArrayArrayMap

◆ ShapeToIntArrayMap

Constructor & Destructor Documentation

◆ PreconditionerLowEnergy()

◆ ~PreconditionerLowEnergy()

Member Function Documentation

◆ create()

◆ CreateRefHexGeom()

◆ CreateRefPrismGeom()

◆ CreateRefPyrGeom()

◆ CreateRefTetGeom()

◆ CreateVariablePMask()

◆ ExtractLocMat()

◆ ReSetPrismMaxRMat()

◆ ReSetTetMaxRMat()

◆ SetupBlockTransformationMatrix()

◆ SetUpPyrMaxRMat()

◆ SetUpReferenceElements()

◆ v_BuildPreconditioner()

◆ v_DoPreconditioner()

◆ v_DoTransformBasisFromLowEnergy()

◆ v_DoTransformBasisToLowEnergy()

◆ v_DoTransformCoeffsFromLowEnergy()

◆ v_DoTransformCoeffsToLowEnergy()

◆ v_InitObject()

◆ v_TransformedSchurCompl()

Member Data Documentation

◆ className

◆ m_BlkMat

◆ m_InvRBlk

◆ m_RBlk

◆ m_sameBlock

◆ m_signChange

◆ m_variablePmask