#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	DoTransformToLowEnergy (Array< OneD, NekDouble > &pInOut, int offset)

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

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

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

void	DoMultiplybyInverseTransposedTransformationMatrix (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 Attributes
DNekBlkMatSharedPtr	m_BlkMat

DNekBlkMatSharedPtr	m_RBlk

DNekBlkMatSharedPtr	m_InvRBlk

Array< OneD, NekDouble >	m_locToGloSignMult

Array< OneD, NekDouble >	m_multiplicity

Array< OneD, int >	m_map

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::ExpansionSharedPtr >	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 (StdRegions::StdExpansionSharedPtr &locExp, DNekScalMatSharedPtr &maxRmat, LocalRegions::ExpansionSharedPtr &expMax, Array< OneD, unsigned int > &vertMapMaxR, Array< OneD, Array< OneD, unsigned int > > &edgeMapMaxR)

void	CreateMultiplicityMap (void)

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

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

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

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

virtual void	v_InitObject ()

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

virtual void	v_DoTransformToLowEnergy (Array< OneD, NekDouble > &pInOut, int offset)
	Transform the solution vector vector to low energy. More...

virtual void	v_DoTransformToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
	Transform the solution vector to low energy form. More...

virtual void	v_DoTransformFromLowEnergy (Array< OneD, NekDouble > &pInOut)
	transform the solution vector from low energy back to the original basis. More...

virtual void	v_DoMultiplybyInverseTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
	Multiply by the block inverse transformation matrix. More...

virtual void	v_DoMultiplybyInverseTransposedTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
	Multiply by the block tranposed inverse transformation matrix. More...

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

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

Additional Inherited Members

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

◆ ShapeToExpMap

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

private

Definition at line 100 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayArrayMap

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

private

Definition at line 105 of file PreconditionerLowEnergy.h.

◆ ShapeToIntArrayMap

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

private

Definition at line 102 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

◆ PreconditionerLowEnergy()

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

Definition at line 67 of file PreconditionerLowEnergy.cpp.

             : Preconditioner(plinsys, pLocToGloMap)
         {
         }

◆ ~PreconditionerLowEnergy()

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

inlinevirtual

Definition at line 74 of file PreconditionerLowEnergy.h.

74 {}

Member Function Documentation

◆ create()

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

inlinestatic

Creates an instance of this class.

Definition at line 57 of file PreconditionerLowEnergy.h.

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

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

◆ CreateMultiplicityMap()

void Nektar::MultiRegions::PreconditionerLowEnergy::CreateMultiplicityMap ( void )

private

Create the inverse multiplicity map.

Definition at line 1441 of file PreconditionerLowEnergy.cpp.

References CG_Iterations::loc, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_locToGloSignMult, m_multiplicity, m_signChange, Vmath::Sdiv(), and sign.

Referenced by v_InitObject().

         {
             auto asmMap = m_locToGloMap.lock();
 
             unsigned int nGlobalBnd = asmMap->GetNumGlobalBndCoeffs();
             unsigned int nEntries   = asmMap->GetNumLocalBndCoeffs();
             unsigned int i;
             
             const Array<OneD, const int> &vMap
                 = asmMap->GetLocalToGlobalBndMap();
 
             const Array< OneD, const NekDouble > &sign 
                 = asmMap->GetLocalToGlobalBndSign();
 
             m_signChange=asmMap->GetSignChange();
 
             // Count the multiplicity of each global DOF on this process
             Array<OneD, NekDouble> vCounts(nGlobalBnd, 0.0);
             if(m_signChange)
             {
                 for (i = 0; i < nEntries; ++i)
                 {
                     if(fabs(sign[i]) > 0.0)
                     {
                         vCounts[vMap[i]] += 1.0;
                     }
                     else // set zero modes to 1 so inverse below does not cause problems
                     {
                         vCounts[vMap[i]] = 1.0;
                     }
                 }
             }
             else
             {
                 for (i = 0; i < nEntries; ++i)
                 {
                     vCounts[vMap[i]] += 1.0;
                 }
             }
                 
             // Get universal multiplicity by globally assembling counts
             asmMap->UniversalAssembleBnd(vCounts);
 
             // Construct a map of 1/multiplicity
             m_locToGloSignMult = Array<OneD, NekDouble>(nEntries);
             for (i = 0; i < nEntries; ++i)
             {
                 if(m_signChange)
                 {
                     m_locToGloSignMult[i] = sign[i]*1.0/vCounts[vMap[i]];
                 }
                 else
                 {
                     m_locToGloSignMult[i] = 1.0/vCounts[vMap[i]];
                 }
             }
             
             int nDirBnd     = asmMap->GetNumGlobalDirBndCoeffs();
             int nGlobHomBnd = nGlobalBnd - nDirBnd;
             int nLocBnd     = asmMap->GetNumLocalBndCoeffs();
 
             //Set up multiplicity array for inverse transposed transformation matrix
             Array<OneD,NekDouble> tmp(nGlobHomBnd,1.0);
             m_multiplicity = Array<OneD,NekDouble>(nGlobHomBnd,1.0);
             Array<OneD,NekDouble> loc(nLocBnd,1.0);
 
             asmMap->GlobalToLocalBnd(tmp,loc, nDirBnd);
             asmMap->AssembleBnd(loc,m_multiplicity, nDirBnd);
             Vmath::Sdiv(nGlobHomBnd,1.0,m_multiplicity,1,m_multiplicity,1);
 
         }

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

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

Referenced by SetUpReferenceElements().

         {
             ////////////////////////////////
             // 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;
         }

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

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

Referenced by SetUpReferenceElements().

         {
             //////////////////////////
             // 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;
         }

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

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

Referenced by SetUpReferenceElements().

         {
             //////////////////////////
             // 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;
         }

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

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

Referenced by SetUpReferenceElements().

         {
             /////////////////////////////////
             // 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;
         }

◆ ExtractLocMat()

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

private

Definition at line 2458 of file PreconditionerLowEnergy.cpp.

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

Referenced by ReSetPrismMaxRMat(), and SetupBlockTransformationMatrix().

         {
             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->GetNfaces();
             
             // 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->GetFaceIntNcoeffs(f);
                 
                 locExp->GetFaceNumModes(f,FwdOrient,m0,m1);
                 
                 Array<OneD, unsigned int> fmapRmat = maxExp->
                     GetFaceInverseBoundaryMap(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;
         }

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

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

Referenced by SetUpReferenceElements().

         {
             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;
         }

◆ ReSetTetMaxRMat()

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

private

Definition at line 2194 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

         {
             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;
         }

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), 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().

         {
             std::shared_ptr<MultiRegions::ExpList> 
                 expList=((m_linsys.lock())->GetLocMat()).lock();
             StdRegions::StdExpansionSharedPtr locExp;
             StdRegions::StdExpansionSharedPtr locExpSav;
             map<int,int> EdgeSize;
             
             int n;
 
             std::map<ShapeType, DNekScalMatSharedPtr>         maxRmat;
             map<ShapeType, LocalRegions::ExpansionSharedPtr > 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);
                 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; 
                     }
                 }
             }
         }

◆ SetUpPyrMaxRMat()

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

private

Definition at line 2058 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

         {
             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;
         }

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

Referenced by SetupBlockTransformationMatrix().

         {
 
             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]) // if Pyramids 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);
                     }
                 }
             }
         }

◆ v_BuildPreconditioner()

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

privatevirtual

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

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

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

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

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 121 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::as(), ASSERTL1, Vmath::Assmb(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Gs::Gather(), 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_locToGloSignMult, m_RBlk, m_signChange, Nektar::LibUtilities::ReduceMax, and Nektar::Transpose().

         {
             std::shared_ptr<MultiRegions::ExpList> 
                 expList=((m_linsys.lock())->GetLocMat()).lock();
             LocalRegions::ExpansionSharedPtr 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.num_elements(); ++i)
             {
                 meshEdgeId=extradiredges[i];
                 edgeDirMap.insert(meshEdgeId);
             }
             
             //Determine which boundary edges and faces have dirichlet values
             for(i = 0; i < bndCondExp.num_elements(); 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->GetNedges(); k++)
                         {
                             meshEdgeId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(k);
                             if(edgeDirMap.count(meshEdgeId) == 0)
                             {
                                 edgeDirMap.insert(meshEdgeId);
                             }
                         }
                         meshFaceId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->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);
 
                 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->GetNfaces();
                 for(j = 0; j < nFaces; ++j)
                 {
                     int nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
                     meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
                         ->GetGeom3D()->GetFid(j);
                     if(FaceSize.count(meshFaceId) == 0)
                     {
                         FaceSize[meshFaceId] = nFaceInteriorCoeffs;
 
                         int m0,m1;
                         locExpansion->GetFaceNumModes(j,locExpansion->GetForient(j),m0,m1);
                         FaceModes[meshFaceId] = pair<int,int>(m0,m1);
                     }
                     else
                     {
                         if(nFaceInteriorCoeffs < FaceSize[meshFaceId])
                         {
                             FaceSize[meshFaceId] =  nFaceInteriorCoeffs;
                             int m0,m1;
                             locExpansion->GetFaceNumModes(j,locExpansion->GetForient(j),m0,m1);
                             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);
 
                 for (j = 0; j < locExpansion->GetNedges(); ++j)
                 {
                     meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->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);
 
                 for (j = 0; j < locExpansion->GetNfaces(); ++j)
                 {
                     meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->
                         GetGeom3D()->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);
                 
                 //loop over the edges of the expansion
                 for(j = 0; j < locExpansion->GetNedges(); ++j)
                 {
                     //get mesh edge id
                     meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
                         ->GetGeom3D()->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->GetNfaces(); ++j)
                 {
                     //get mesh face id
                     meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
                         ->GetGeom3D()->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);
                 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 mask1 = 1.0;
                 NekDouble mask2 = 1.0;
                 NekDouble val;
                 
                 for(int i = 0; i < nCoeffs; ++i)
                 {
                     if(m_signChange)
                     {
                         mask1 = (m_locToGloSignMult[cnt+i] == 0.0)? 0.0:1.0;
                     }
                     for(int j = 0; j < nCoeffs; ++j)
                     {
                         if(m_signChange)
                         {
                             mask2 = (m_locToGloSignMult[cnt+j] == 0.0)? 0.0:1.0;
                         }
                         val = S(i,j)*mask1*mask2;
                         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->as<LocalRegions::Expansion3D>()->GetGeom3D()->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->GetForient(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->GetFaceInverseBoundaryMap
                                 (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.num_elements(),  
                              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.num_elements()-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);
                 }
             }
         }

◆ v_DoMultiplybyInverseTransformationMatrix()

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

privatevirtual

Multiply by the block inverse transformation matrix.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1258 of file PreconditionerLowEnergy.cpp.

References ASSERTL1, Blas::Dgemm(), Nektar::eWrapper, Vmath::Gathr(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_locToGloSignMult, m_map, m_sameBlock, and Vmath::Vcopy().

         {
             auto asmMap = m_locToGloMap.lock();
 
             int nGlobBndDofs       = asmMap->GetNumGlobalBndCoeffs();
             int nDirBndDofs        = asmMap->GetNumGlobalDirBndCoeffs();
             int nGlobHomBndDofs    = nGlobBndDofs - nDirBndDofs;
             int nLocBndDofs        = asmMap->GetNumLocalBndCoeffs();
 
             ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
                      "Input array is greater than the nGlobHomBndDofs");
             ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
                      "Output array is greater than the nGlobHomBndDofs");
 
             //vectors of length number of non-dirichlet boundary dofs
             NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
             NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
                                           eWrapper);
             //Block inverse transformation matrix
             DNekBlkMat &invR = *m_InvRBlk;
 
             Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
             NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
             m_map = asmMap->GetLocalToGlobalBndMap();
             Array<OneD, NekDouble> pLocalIn(nLocBndDofs, 0.0);
 
             // Allocated array of size number of global boundary dofs and copy
             // the input array to the tmp array offset by Dirichlet boundary
             // conditions.
             Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
             Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
 
             //Global boundary dofs (with zeroed dirichlet values) to
             //local boundary dofs
             Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(),
                          tmp.get(), m_map.get(), pLocalIn.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, pLocal.get() + cnt, nbndcoeffs);
         cnt  += nbndcoeffs*nexp;
         cnt1 += nexp;
         }
             
 
             //Assemble local boundary to global non-dirichlet boundary
             asmMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
 
     }

◆ v_DoMultiplybyInverseTransposedTransformationMatrix()

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

privatevirtual

Multiply by the block tranposed inverse transformation matrix.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1321 of file PreconditionerLowEnergy.cpp.

References ASSERTL1, Blas::Dgemm(), Nektar::eWrapper, m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_map, m_multiplicity, m_sameBlock, v_TransformedSchurCompl(), and Vmath::Vmul().

         {
             auto asmMap = m_locToGloMap.lock();
 
             int nGlobBndDofs       = asmMap->GetNumGlobalBndCoeffs();
             int nDirBndDofs        = asmMap->GetNumGlobalDirBndCoeffs();
             int nGlobHomBndDofs    = nGlobBndDofs - nDirBndDofs;
             int nLocBndDofs        = asmMap->GetNumLocalBndCoeffs();
 
             ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
                      "Input array is greater than the nGlobHomBndDofs");
             ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
                      "Output array is greater than the nGlobHomBndDofs");
 
             //vectors of length number of non-dirichlet boundary dofs
             NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
             NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
                                           eWrapper);
             //Block inverse transformation matrix
             DNekBlkMat &invR = *m_InvRBlk;
 
             Array<OneD, NekDouble> pLocalIn(nLocBndDofs, 0.0);
             NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocalIn,eWrapper);
             m_map = asmMap->GetLocalToGlobalBndMap();
             Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
 
             asmMap->GlobalToLocalBnd(pInput,pLocalIn, nDirBndDofs);
 
 
             //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, pLocal.get() + cnt, nbndcoeffs);
         cnt  += nbndcoeffs*nexp;
         cnt1 += nexp;
         }
 
             
             asmMap->AssembleBnd(pLocal,pOutput, nDirBndDofs);
 
             Vmath::Vmul(nGlobHomBndDofs,pOutput,1,m_multiplicity,1,pOutput,1);
     }

◆ v_DoPreconditioner()

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

privatevirtual

Apply the low energy preconditioner during the conjugate gradient routine

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 936 of file PreconditionerLowEnergy.cpp.

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

         {
             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;
     }

◆ v_DoTransformFromLowEnergy()

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

privatevirtual

transform the solution vector from low energy back to the original basis.

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

References ASSERTL1, Vmath::Assmb(), Blas::Dgemm(), Nektar::eWrapper, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_locToGloSignMult, m_map, m_RBlk, m_sameBlock, and Vmath::Vcopy().

         {
             auto asmMap = m_locToGloMap.lock();
 
             int nGlobBndDofs       = asmMap->GetNumGlobalBndCoeffs();
             int nDirBndDofs        = asmMap->GetNumGlobalDirBndCoeffs();
             int nGlobHomBndDofs    = nGlobBndDofs - nDirBndDofs;
             int nLocBndDofs        = asmMap->GetNumLocalBndCoeffs();
 
             ASSERTL1(pInOut.num_elements() >= nGlobBndDofs,
                      "Output array is greater than the nGlobBndDofs");
 
             //Block  transformation matrix
             DNekBlkMat &R = *m_RBlk;
 
             Array<OneD, NekDouble> tmpOffset = pInOut + nDirBndDofs;
             NekVector<NekDouble> V_GlobHomBnd(nGlobHomBndDofs, tmpOffset, eWrapper);
 
             Array<OneD, NekDouble> pLocalIn(nLocBndDofs, 0.0);
             NekVector<NekDouble> V_LocBnd(nLocBndDofs,pLocalIn,eWrapper);
             m_map = asmMap->GetLocalToGlobalBndMap();
             Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
             Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
 
             //Global boundary (less dirichlet) to local boundary
             asmMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
 
             //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, pLocal.get() + cnt, nbndcoeffs);
         cnt  += nbndcoeffs*nexp;
         cnt1 += nexp;
         }
 
             //Assemble local boundary to global boundary
             Vmath::Assmb(nLocBndDofs, m_locToGloSignMult.get(),pLocal.get(), m_map.get(), tmp.get());
 
             //Universal assemble across processors
             asmMap->UniversalAssembleBnd(tmp);
             
             //copy non-dirichlet boundary values
             Vmath::Vcopy(nGlobBndDofs-nDirBndDofs, tmp.get() + nDirBndDofs, 1, pInOut.get() + nDirBndDofs, 1);
         }

◆ v_DoTransformToLowEnergy() [1/2]

void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformToLowEnergy	(	Array< OneD, NekDouble > &	pInOut,
		int	offset
	)

privatevirtual

Transform the solution vector vector to low energy.

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

References Blas::Dgemm(), Nektar::eWrapper, Vmath::Gathr(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_locToGloSignMult, m_map, m_RBlk, m_sameBlock, and Vmath::Vcopy().

         {
             auto asmMap = m_locToGloMap.lock();
             int nGlobBndDofs       = asmMap->GetNumGlobalBndCoeffs();
             int nDirBndDofs        = asmMap->GetNumGlobalDirBndCoeffs();
             int nGlobHomBndDofs    = nGlobBndDofs - nDirBndDofs;
             int nLocBndDofs        = asmMap->GetNumLocalBndCoeffs();
 
             //Non-dirichlet boundary dofs
             Array<OneD, NekDouble> tmpOffset = pInOut + offset;
             NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs, tmpOffset, eWrapper);
 
             //Block transformation matrix
             DNekBlkMat &R = *m_RBlk;
 
             Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
             NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
             m_map = asmMap->GetLocalToGlobalBndMap();
             Array<OneD, NekDouble> pLocalIn(nLocBndDofs, 0.0);
 
             //Not actually needed but we should only work with the
             //Global boundary dofs
             Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
             Vmath::Vcopy(nGlobBndDofs, pInOut.get(), 1, tmp.get(), 1);
 
             //Global boundary (with dirichlet values) to local
             //boundary with multiplicity
             Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(),
                          tmp.get(), m_map.get(), pLocalIn.get());
 
         //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, pLocal.get() + cnt, nbndcoeffs);
         cnt  += nbndcoeffs*nexp;
         cnt1 += nexp;
         }
 
             //Assemble local boundary to global non-dirichlet Dofs
             asmMap->AssembleBnd(F_LocBnd,F_HomBnd, nDirBndDofs);
         }

◆ v_DoTransformToLowEnergy() [2/2]

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

privatevirtual

Transform the solution vector to low energy form.

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

References ASSERTL1, Blas::Dgemm(), Nektar::eWrapper, Vmath::Gathr(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_locToGloSignMult, m_map, m_RBlk, m_sameBlock, and Vmath::Vcopy().

         {
             auto asmMap = m_locToGloMap.lock();
 
             int nGlobBndDofs       = asmMap->GetNumGlobalBndCoeffs();
             int nDirBndDofs        = asmMap->GetNumGlobalDirBndCoeffs();
             int nGlobHomBndDofs    = nGlobBndDofs - nDirBndDofs;
             int nLocBndDofs        = asmMap->GetNumLocalBndCoeffs();
 
             //Input/output vectors should be length nGlobHomBndDofs
             ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
                      "Input array is greater than the nGlobHomBndDofs");
             ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
                      "Output array is greater than the nGlobHomBndDofs");
 
             //vectors of length number of non-dirichlet boundary dofs
             NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
                                           eWrapper);
             //Block transformation matrix
             DNekBlkMat &R = *m_RBlk;
 
             Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
             NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
             m_map = asmMap->GetLocalToGlobalBndMap();
             Array<OneD, NekDouble> pLocalIn(nLocBndDofs, 0.0);
 
             // Allocated array of size number of global boundary dofs and copy
             // the input array to the tmp array offset by Dirichlet boundary
             // conditions.
             Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
             Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get()
                          + nDirBndDofs, 1);
             
             //Global boundary dofs (with zeroed dirichlet values) to
             //local boundary dofs - This also divides by the mulplicity
             Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(),
                          tmp.get(), m_map.get(), pLocalIn.get());
 
             //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, pLocal.get() + cnt, nbndcoeffs);
         cnt  += nbndcoeffs*nexp;
         cnt1 += nexp;
         }
             
             //Assemble local boundary to global non-dirichlet boundary
             asmMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
         }

◆ v_InitObject()

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

privatevirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 74 of file PreconditionerLowEnergy.cpp.

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

         {
             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 multiplicity map for transformation from global to local
             CreateMultiplicityMap();
     }

◆ v_TransformedSchurCompl()

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

privatevirtual

Set up the transformed block matrix system.

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1382 of file PreconditionerLowEnergy.cpp.

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

Referenced by v_DoMultiplybyInverseTransposedTransformationMatrix().

     {
             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 mask1 = 1.0;
             NekDouble mask2 = 1.0;
             NekDouble val;
             
             for(int i = 0; i < nbnd; ++i)
             {
                 if(m_signChange)
                 {
                     mask1 = (m_locToGloSignMult[bndoffset+i] == 0.0)? 0.0:1.0;
                 }
                 for(int j = 0; j < nbnd; ++j)
                 {
                     if(m_signChange)
                     {
                         mask2 = (m_locToGloSignMult[bndoffset+j] == 0.0)? 0.0:1.0;
                     }
                     val = S1(i,j)*mask1*mask2;
                     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;
     }

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

◆ m_BlkMat

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_BlkMat

protected

Definition at line 78 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

◆ m_InvRBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_InvRBlk

protected

Definition at line 80 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_DoMultiplybyInverseTransformationMatrix(), and v_DoMultiplybyInverseTransposedTransformationMatrix().

◆ m_locToGloSignMult

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

protected

Definition at line 83 of file PreconditionerLowEnergy.h.

Referenced by CreateMultiplicityMap(), v_BuildPreconditioner(), v_DoMultiplybyInverseTransformationMatrix(), v_DoTransformFromLowEnergy(), v_DoTransformToLowEnergy(), and v_TransformedSchurCompl().

◆ m_map

Array<OneD, int> Nektar::MultiRegions::PreconditionerLowEnergy::m_map

protected

Definition at line 85 of file PreconditionerLowEnergy.h.

Referenced by v_DoMultiplybyInverseTransformationMatrix(), v_DoMultiplybyInverseTransposedTransformationMatrix(), v_DoTransformFromLowEnergy(), and v_DoTransformToLowEnergy().

◆ m_multiplicity

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

protected

Definition at line 84 of file PreconditionerLowEnergy.h.

Referenced by CreateMultiplicityMap(), and v_DoMultiplybyInverseTransposedTransformationMatrix().

◆ m_RBlk

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RBlk

protected

Definition at line 79 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_BuildPreconditioner(), v_DoTransformFromLowEnergy(), v_DoTransformToLowEnergy(), 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_DoMultiplybyInverseTransformationMatrix(), v_DoMultiplybyInverseTransposedTransformationMatrix(), v_DoTransformFromLowEnergy(), and v_DoTransformToLowEnergy().

◆ m_signChange

bool Nektar::MultiRegions::PreconditionerLowEnergy::m_signChange

protected

Definition at line 87 of file PreconditionerLowEnergy.h.

Referenced by CreateMultiplicityMap(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Attributes

Private Types

Private Member Functions

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ ShapeToDNekMap

◆ ShapeToExpMap

◆ ShapeToIntArrayArrayMap

◆ ShapeToIntArrayMap

Constructor & Destructor Documentation

◆ PreconditionerLowEnergy()

◆ ~PreconditionerLowEnergy()

Member Function Documentation

◆ create()

◆ CreateMultiplicityMap()

◆ CreateRefHexGeom()

◆ CreateRefPrismGeom()

◆ CreateRefPyrGeom()

◆ CreateRefTetGeom()

◆ ExtractLocMat()

◆ ReSetPrismMaxRMat()

◆ ReSetTetMaxRMat()

◆ SetupBlockTransformationMatrix()

◆ SetUpPyrMaxRMat()

◆ SetUpReferenceElements()

◆ v_BuildPreconditioner()

◆ v_DoMultiplybyInverseTransformationMatrix()

◆ v_DoMultiplybyInverseTransposedTransformationMatrix()

◆ v_DoPreconditioner()

◆ v_DoTransformFromLowEnergy()

◆ v_DoTransformToLowEnergy() [1/2]

◆ v_DoTransformToLowEnergy() [2/2]

◆ v_InitObject()

◆ v_TransformedSchurCompl()

Member Data Documentation

◆ className

◆ m_BlkMat

◆ m_InvRBlk

◆ m_locToGloSignMult

◆ m_map

◆ m_multiplicity

◆ m_RBlk

◆ m_sameBlock

◆ m_signChange