Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Attributes | Private Member Functions | List of all members
Nektar::MultiRegions::PreconditionerLowEnergy Class Reference

#include <PreconditionerLowEnergy.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerLowEnergy:
Inheritance graph
[legend]
Collaboration diagram for Nektar::MultiRegions::PreconditionerLowEnergy:
Collaboration graph
[legend]

Public Member Functions

 PreconditionerLowEnergy (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~PreconditionerLowEnergy ()
 
- Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
 Preconditioner (const boost::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, const boost::shared_ptr< DNekScalMat > &loc_mat)
 

Static Public Member Functions

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

Static Public Attributes

static std::string className
 Name of class. More...
 

Protected Attributes

const boost::weak_ptr
< GlobalLinSys > 		m_linsys
 
boost::shared_ptr< AssemblyMap > m_locToGloMap
 
DNekBlkMatSharedPtr m_BlkMat
 
DNekScalMatSharedPtr m_bnd_mat
 
DNekScalBlkMatSharedPtr m_RBlk
 
DNekScalBlkMatSharedPtr m_RTBlk
 
DNekScalBlkMatSharedPtr m_InvRBlk
 
DNekScalBlkMatSharedPtr m_InvRTBlk
 
DNekScalMatSharedPtr m_Rtet
 
DNekScalMatSharedPtr m_RTtet
 
DNekScalMatSharedPtr m_Rinvtet
 
DNekScalMatSharedPtr m_RTinvtet
 
DNekScalMatSharedPtr m_Rhex
 
DNekScalMatSharedPtr m_RThex
 
DNekScalMatSharedPtr m_Rinvhex
 
DNekScalMatSharedPtr m_RTinvhex
 
DNekScalMatSharedPtr m_Rprism
 
DNekScalMatSharedPtr m_RTprism
 
DNekScalMatSharedPtr m_Rinvprism
 
DNekScalMatSharedPtr m_RTinvprism
 
Array< OneD, NekDouble > m_locToGloSignMult
 
Array< OneD, NekDouble > m_multiplicity
 
Array< OneD, int > m_map
 
- Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const boost::weak_ptr
< GlobalLinSys > 		m_linsys
 
PreconditionerType m_preconType
 
DNekMatSharedPtr m_preconditioner
 
boost::shared_ptr< AssemblyMap > m_locToGloMap
 
LibUtilities::CommSharedPtr m_comm
 

Private Member Functions

void SetUpReferenceElements (void)
 Sets up the reference elements needed by the preconditioner. More...
 
void CreateMultiplicityMap (void)
 
void SetupBlockTransformationMatrix (void)
 
void ModifyPrismTransformationMatrix (LocalRegions::TetExpSharedPtr TetExp, LocalRegions::PrismExpSharedPtr PrismExp, DNekMatSharedPtr Rmodprism, DNekMatSharedPtr RTmodprism)
 Modify the prism transformation matrix to align with the tetrahedral modes. More...
 
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom (void)
 Sets up the reference tretrahedral element needed to construct a low energy basis. 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 vector to low energy. 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 offset, const boost::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.

Constructor & Destructor Documentation

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

Definition at line 64 of file PreconditionerLowEnergy.cpp.

67  : Preconditioner(plinsys, pLocToGloMap),
68  m_linsys(plinsys),
69  m_locToGloMap(pLocToGloMap)
70  {
71  }
Nektar::MultiRegions::Preconditioner::Preconditioner
Preconditioner(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
Definition: Preconditioner.cpp:76
Nektar::MultiRegions::PreconditionerLowEnergy::m_linsys
const boost::weak_ptr< GlobalLinSys > m_linsys
Definition: PreconditionerLowEnergy.h:79
virtual Nektar::MultiRegions::PreconditionerLowEnergy::~PreconditionerLowEnergy ( )
inlinevirtual

Definition at line 75 of file PreconditionerLowEnergy.h.

75 {}

Member Function Documentation

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

Creates an instance of this class.

Definition at line 57 of file PreconditionerLowEnergy.h.

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

61  {
62  PreconditionerSharedPtr p = MemoryManager<PreconditionerLowEnergy>::AllocateSharedPtr(plinsys,pLocToGloMap);
63  p->InitObject();
64  return p;
65  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::MultiRegions::PreconditionerSharedPtr
boost::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:62
void Nektar::MultiRegions::PreconditionerLowEnergy::CreateMultiplicityMap ( void  )
private

Create the inverse multiplicity map.

Definition at line 1272 of file PreconditionerLowEnergy.cpp.

References m_locToGloMap, m_locToGloSignMult, m_multiplicity, Vmath::Sdiv(), and sign.

Referenced by v_InitObject().

1273  {
1274  unsigned int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
1275  unsigned int nEntries = m_locToGloMap->GetNumLocalBndCoeffs();
1276  unsigned int i;
1277 
1278  const Array<OneD, const int> &vMap
1279  = m_locToGloMap->GetLocalToGlobalBndMap();
1280 
1281  const Array< OneD, const NekDouble > &sign
1282  = m_locToGloMap->GetLocalToGlobalBndSign();
1283 
1284  bool m_signChange=m_locToGloMap->GetSignChange();
1285 
1286  // Count the multiplicity of each global DOF on this process
1287  Array<OneD, NekDouble> vCounts(nGlobalBnd, 0.0);
1288  for (i = 0; i < nEntries; ++i)
1289  {
1290  vCounts[vMap[i]] += 1.0;
1291  }
1292 
1293  // Get universal multiplicity by globally assembling counts
1294  m_locToGloMap->UniversalAssembleBnd(vCounts);
1295 
1296  // Construct a map of 1/multiplicity
1297  m_locToGloSignMult = Array<OneD, NekDouble>(nEntries);
1298  for (i = 0; i < nEntries; ++i)
1299  {
1300  if(m_signChange)
1301  {
1302  m_locToGloSignMult[i] = sign[i]*1.0/vCounts[vMap[i]];
1303  }
1304  else
1305  {
1306  m_locToGloSignMult[i] = 1.0/vCounts[vMap[i]];
1307  }
1308  }
1309 
1310  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1311  int nGlobHomBnd = nGlobalBnd - nDirBnd;
1312  int nLocBnd = m_locToGloMap->GetNumLocalBndCoeffs();
1313 
1314  //Set up multiplicity array for inverse transposed transformation matrix
1315  Array<OneD,NekDouble> tmp(nGlobHomBnd,1.0);
1316  m_multiplicity = Array<OneD,NekDouble>(nGlobHomBnd,1.0);
1317  Array<OneD,NekDouble> loc(nLocBnd,1.0);
1318 
1319  m_locToGloMap->GlobalToLocalBnd(tmp,loc, nDirBnd);
1320  m_locToGloMap->AssembleBnd(loc,m_multiplicity, nDirBnd);
1321  Vmath::Sdiv(nGlobHomBnd,1.0,m_multiplicity,1,m_multiplicity,1);
1322 
1323  }
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
SpatialDomains::HexGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefHexGeom ( void  )
private

Sets up the reference hexahedral element needed to construct a low energy basis.

Definition at line 1511 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::eBackwards, and Nektar::StdRegions::eForwards.

Referenced by SetUpReferenceElements().

1512  {
1513  ////////////////////////////////
1514  // Set up Hexahedron vertices //
1515  ////////////////////////////////
1516 
1517  const int three=3;
1518 
1519  const int nVerts = 8;
1520  const double point[][3] = {
1521  {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1522  {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1523  };
1524 
1525  // Populate the list of verts
1526  SpatialDomains::PointGeomSharedPtr verts[8];
1527  for( int i = 0; i < nVerts; ++i ) {
1528  verts[i] = MemoryManager<SpatialDomains::PointGeom>
1529  ::AllocateSharedPtr(three, i, point[i][0],
1530  point[i][1], point[i][2]);
1531  }
1532 
1533  /////////////////////////////
1534  // Set up Hexahedron Edges //
1535  /////////////////////////////
1536 
1537  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1538  const int nEdges = 12;
1539  const int vertexConnectivity[][2] = {
1540  {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1541  {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1542  };
1543 
1544  // Populate the list of edges
1545  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1546  for( int i = 0; i < nEdges; ++i ) {
1547  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1548  for( int j = 0; j < 2; ++j ) {
1549  vertsArray[j] = verts[vertexConnectivity[i][j]];
1550  }
1551  edges[i] = MemoryManager<SpatialDomains::SegGeom>::
1552  AllocateSharedPtr( i, three, vertsArray);
1553  }
1554 
1555  /////////////////////////////
1556  // Set up Hexahedron faces //
1557  /////////////////////////////
1558 
1559  const int nFaces = 6;
1560  const int edgeConnectivity[][4] = {
1561  {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1562  {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1563  };
1564  const bool isEdgeFlipped[][4] = {
1565  {0,0,0,1}, {0,0,1,1}, {0,0,1,1},
1566  {0,0,1,1}, {0,0,1,1}, {0,0,0,1}
1567  };
1568 
1569  // Populate the list of faces
1570  SpatialDomains::QuadGeomSharedPtr faces[nFaces];
1571  for( int i = 0; i < nFaces; ++i ) {
1572  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1573  StdRegions::Orientation eorientArray[4];
1574  for( int j = 0; j < 4; ++j ) {
1575  edgeArray[j] = edges[edgeConnectivity[i][j]];
1576  eorientArray[j] = isEdgeFlipped[i][j] ?
1577  StdRegions::eBackwards : StdRegions::eForwards;
1578  }
1579  faces[i] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(i, edgeArray,
1580  eorientArray);
1581  }
1582 
1583  SpatialDomains::HexGeomSharedPtr geom =
1584  MemoryManager<SpatialDomains::HexGeom>::AllocateSharedPtr
1585  (faces);
1586 
1587  geom->SetOwnData();
1588 
1589  return geom;
1590  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::SpatialDomains::QuadGeomSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
Nektar::SpatialDomains::HexGeomSharedPtr
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Nektar::SpatialDomains::SegGeomSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Nektar::SpatialDomains::PointGeomSharedPtr
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
SpatialDomains::PrismGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPrismGeom ( void  )
private

Sets up the reference prismatic element needed to construct a low energy basis.

Definition at line 1329 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::eBackwards, and Nektar::StdRegions::eForwards.

Referenced by SetUpReferenceElements().

1330  {
1331  //////////////////////////
1332  // Set up Prism element //
1333  //////////////////////////
1334 
1335  const int three=3;
1336  const int nVerts = 6;
1337  const double point[][3] = {
1338  {-1,-1,0}, {1,-1,0}, {1,1,0},
1339  {-1,1,0}, {0,-1,sqrt(double(3))}, {0,1,sqrt(double(3))},
1340  };
1341 
1342  //boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
1343  SpatialDomains::PointGeomSharedPtr verts[6];
1344  for(int i=0; i < nVerts; ++i)
1345  {
1346  verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr
1347  ( three, i, point[i][0], point[i][1], point[i][2] );
1348  }
1349  const int nEdges = 9;
1350  const int vertexConnectivity[][2] = {
1351  {0,1}, {1,2}, {3,2}, {0,3}, {0,4},
1352  {1,4}, {2,5}, {3,5}, {4,5}
1353  };
1354 
1355  // Populate the list of edges
1356  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1357  for(int i=0; i < nEdges; ++i){
1358  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1359  for(int j=0; j<2; ++j)
1360  {
1361  vertsArray[j] = verts[vertexConnectivity[i][j]];
1362  }
1363  edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(i, three, vertsArray);
1364  }
1365 
1366  ////////////////////////
1367  // Set up Prism faces //
1368  ////////////////////////
1369 
1370  const int nFaces = 5;
1371  //quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
1372  const int quadEdgeConnectivity[][4] = { {0,1,2,3}, {1,6,8,5}, {3,7,8,4} };
1373  const bool isQuadEdgeFlipped[][4] = { {0,0,1,1}, {0,0,1,1}, {0,0,1,1} };
1374  // QuadId ordered as 0, 1, 2, otherwise return false
1375  const int quadId[] = { 0,-1,1,-1,2 };
1376 
1377  //triangle-edge connectivity side-triface-1, side triface-3
1378  const int triEdgeConnectivity[][3] = { {0,5,4}, {2,6,7} };
1379  const bool isTriEdgeFlipped[][3] = { {0,0,1}, {0,0,1} };
1380  // TriId ordered as 0, 1, otherwise return false
1381  const int triId[] = { -1,0,-1,1,-1 };
1382 
1383  // Populate the list of faces
1384  SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1385  for(int f = 0; f < nFaces; ++f){
1386  if(f == 1 || f == 3) {
1387  int i = triId[f];
1388  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1389  StdRegions::Orientation eorientArray[3];
1390  for(int j = 0; j < 3; ++j){
1391  edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1392  eorientArray[j] = isTriEdgeFlipped[i][j] ? StdRegions::eBackwards : StdRegions::eForwards;
1393  }
1394  faces[f] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(f, edgeArray, eorientArray);
1395  }
1396  else {
1397  int i = quadId[f];
1398  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1399  StdRegions::Orientation eorientArray[4];
1400  for(int j=0; j < 4; ++j){
1401  edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1402  eorientArray[j] = isQuadEdgeFlipped[i][j] ? StdRegions::eBackwards : StdRegions::eForwards;
1403  }
1404  faces[f] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(f, edgeArray, eorientArray);
1405  }
1406  }
1407 
1408  SpatialDomains::PrismGeomSharedPtr geom = MemoryManager<SpatialDomains::PrismGeom>::AllocateSharedPtr(faces);
1409 
1410  geom->SetOwnData();
1411 
1412  return geom;
1413  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::SpatialDomains::SegGeomSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Nektar::SpatialDomains::Geometry2DSharedPtr
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
Nektar::SpatialDomains::PrismGeomSharedPtr
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
Nektar::SpatialDomains::PointGeomSharedPtr
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
SpatialDomains::TetGeomSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefTetGeom ( void  )
private

Sets up the reference tretrahedral element needed to construct a low energy basis.

Definition at line 1419 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::eBackwards, and Nektar::StdRegions::eForwards.

Referenced by SetUpReferenceElements().

1420  {
1421  /////////////////////////////////
1422  // Set up Tetrahedron vertices //
1423  /////////////////////////////////
1424 
1425  int i,j;
1426  const int three=3;
1427  const int nVerts = 4;
1428  const double point[][3] = {
1429  {-1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1430  {1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1431  {0,2/sqrt(double(3)),-1/sqrt(double(6))},
1432  {0,0,3/sqrt(double(6))}};
1433 
1434  boost::shared_ptr<SpatialDomains::PointGeom> verts[4];
1435  for(i=0; i < nVerts; ++i)
1436  {
1437  verts[i] =
1438  MemoryManager<SpatialDomains::PointGeom>::
1439  AllocateSharedPtr
1440  ( three, i, point[i][0], point[i][1], point[i][2] );
1441  }
1442 
1443  //////////////////////////////
1444  // Set up Tetrahedron Edges //
1445  //////////////////////////////
1446 
1447  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1448  const int nEdges = 6;
1449  const int vertexConnectivity[][2] = {
1450  {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1451  };
1452 
1453  // Populate the list of edges
1454  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1455  for(i=0; i < nEdges; ++i)
1456  {
1457  boost::shared_ptr<SpatialDomains::PointGeom>
1458  vertsArray[2];
1459  for(j=0; j<2; ++j)
1460  {
1461  vertsArray[j] = verts[vertexConnectivity[i][j]];
1462  }
1463 
1464  edges[i] = MemoryManager<SpatialDomains::SegGeom>
1465  ::AllocateSharedPtr(i, three, vertsArray);
1466  }
1467 
1468  //////////////////////////////
1469  // Set up Tetrahedron faces //
1470  //////////////////////////////
1471 
1472  const int nFaces = 4;
1473  const int edgeConnectivity[][3] = {
1474  {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1475  };
1476  const bool isEdgeFlipped[][3] = {
1477  {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
1478  };
1479 
1480  // Populate the list of faces
1481  SpatialDomains::TriGeomSharedPtr faces[nFaces];
1482  for(i=0; i < nFaces; ++i)
1483  {
1484  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1485  StdRegions::Orientation eorientArray[3];
1486  for(j=0; j < 3; ++j)
1487  {
1488  edgeArray[j] = edges[edgeConnectivity[i][j]];
1489  eorientArray[j] = isEdgeFlipped[i][j] ?
1490  StdRegions::eBackwards : StdRegions::eForwards;
1491  }
1492 
1493 
1494  faces[i] = MemoryManager<SpatialDomains::TriGeom>
1495  ::AllocateSharedPtr(i, edgeArray, eorientArray);
1496  }
1497 
1498  SpatialDomains::TetGeomSharedPtr geom =
1499  MemoryManager<SpatialDomains::TetGeom>::AllocateSharedPtr
1500  (faces);
1501 
1502  geom->SetOwnData();
1503 
1504  return geom;
1505  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::SpatialDomains::SegGeomSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Nektar::SpatialDomains::TetGeomSharedPtr
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
Nektar::SpatialDomains::TriGeomSharedPtr
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
void Nektar::MultiRegions::PreconditionerLowEnergy::ModifyPrismTransformationMatrix ( LocalRegions::TetExpSharedPtr  TetExp,
LocalRegions::PrismExpSharedPtr  PrismExp,
DNekMatSharedPtr  Rmodprism,
DNekMatSharedPtr  RTmodprism 
)
private

Modify the prism transformation matrix to align with the tetrahedral modes.

This routine replaces the edge and triangular face components of the prismatic vertex transformation matrices $\mathbf{R}_{ve}$ and $\mathbf{R}_{vf}$ with the corresponding components from the tetrahedral transformation matrices. Additionally, triangular face components in the prismatic edge transformation matrix $\mathbf{R}_{ef}$ with the corresponding component from the tetrahedral transformation matrix.

Definition at line 1868 of file PreconditionerLowEnergy.cpp.

Referenced by SetUpReferenceElements().

1873  {
1874  NekDouble Rvalue, RTvalue;
1875  int i, j;
1876 
1877  //For a tet element the bottom face is made up of the following:
1878  //vertices: 0, 1 and 2 edges: 0, 1 and 2 face: 0. We first need to
1879  //determine the mode locations of these vertices, edges and face so
1880  //we can extract the correct values from the tetrahedral R matrix.
1881 
1882  //These are the vertex mode locations of R which need to be replaced
1883  //in the prism element
1884  int TetVertex0=TetExp->GetVertexMap(0);
1885  int TetVertex1=TetExp->GetVertexMap(1);
1886  int TetVertex2=TetExp->GetVertexMap(2);
1887  int TetVertex3=TetExp->GetVertexMap(3);
1888 
1889 
1890  //These are the edge mode locations of R which need to be replaced
1891  //in the prism element
1892  Array<OneD, unsigned int> TetEdge0=TetExp->GetEdgeInverseBoundaryMap(0);
1893  Array<OneD, unsigned int> TetEdge1=TetExp->GetEdgeInverseBoundaryMap(1);
1894  Array<OneD, unsigned int> TetEdge2=TetExp->GetEdgeInverseBoundaryMap(2);
1895  Array<OneD, unsigned int> TetEdge3=TetExp->GetEdgeInverseBoundaryMap(3);
1896  Array<OneD, unsigned int> TetEdge4=TetExp->GetEdgeInverseBoundaryMap(4);
1897  Array<OneD, unsigned int> TetEdge5=TetExp->GetEdgeInverseBoundaryMap(5);
1898 
1899  //These are the face mode locations of R which need to be replaced
1900  //in the prism element
1901  Array<OneD, unsigned int> TetFace=TetExp->GetFaceInverseBoundaryMap(1);
1902 
1903  //Prism vertex modes
1904  int PrismVertex0=PrismExp->GetVertexMap(0);
1905  int PrismVertex1=PrismExp->GetVertexMap(1);
1906  int PrismVertex2=PrismExp->GetVertexMap(2);
1907  int PrismVertex3=PrismExp->GetVertexMap(3);
1908  int PrismVertex4=PrismExp->GetVertexMap(4);
1909  int PrismVertex5=PrismExp->GetVertexMap(5);
1910 
1911  //Prism edge modes
1912  Array<OneD, unsigned int> PrismEdge0=
1913  PrismExp->GetEdgeInverseBoundaryMap(0);
1914  Array<OneD, unsigned int> PrismEdge1=
1915  PrismExp->GetEdgeInverseBoundaryMap(1);
1916  Array<OneD, unsigned int> PrismEdge2=
1917  PrismExp->GetEdgeInverseBoundaryMap(2);
1918  Array<OneD, unsigned int> PrismEdge3=
1919  PrismExp->GetEdgeInverseBoundaryMap(3);
1920  Array<OneD, unsigned int> PrismEdge4=
1921  PrismExp->GetEdgeInverseBoundaryMap(4);
1922  Array<OneD, unsigned int> PrismEdge5=
1923  PrismExp->GetEdgeInverseBoundaryMap(5);
1924  Array<OneD, unsigned int> PrismEdge6=
1925  PrismExp->GetEdgeInverseBoundaryMap(6);
1926  Array<OneD, unsigned int> PrismEdge7=
1927  PrismExp->GetEdgeInverseBoundaryMap(7);
1928  Array<OneD, unsigned int> PrismEdge8=
1929  PrismExp->GetEdgeInverseBoundaryMap(8);
1930 
1931  //Prism face 1 & 3 face modes
1932  Array<OneD, unsigned int> PrismFace1=
1933  PrismExp->GetFaceInverseBoundaryMap(1);
1934  Array<OneD, unsigned int> PrismFace3=
1935  PrismExp->GetFaceInverseBoundaryMap(3);
1936  Array<OneD, unsigned int> PrismFace0=
1937  PrismExp->GetFaceInverseBoundaryMap(0);
1938  Array<OneD, unsigned int> PrismFace2=
1939  PrismExp->GetFaceInverseBoundaryMap(2);
1940  Array<OneD, unsigned int> PrismFace4=
1941  PrismExp->GetFaceInverseBoundaryMap(4);
1942 
1943  //vertex 0 edge 0 3 & 4
1944  for(i=0; i< PrismEdge0.num_elements(); ++i)
1945  {
1946  Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
1947  Rmodprism->SetValue(PrismVertex0,PrismEdge0[i],Rvalue);
1948  Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
1949  Rmodprism->SetValue(PrismVertex0,PrismEdge3[i],Rvalue);
1950  Rvalue=(*m_Rtet)(TetVertex0,TetEdge3[i]);
1951  Rmodprism->SetValue(PrismVertex0,PrismEdge4[i],Rvalue);
1952 
1953  //transposed values
1954  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
1955  RTmodprism->SetValue(PrismEdge0[i],PrismVertex0,RTvalue);
1956  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
1957  RTmodprism->SetValue(PrismEdge3[i],PrismVertex0,RTvalue);
1958  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex0);
1959  RTmodprism->SetValue(PrismEdge4[i],PrismVertex0,RTvalue);
1960  }
1961 
1962  //vertex 1 edge 0 1 & 5
1963  for(i=0; i< PrismEdge1.num_elements(); ++i)
1964  {
1965  Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1966  Rmodprism->SetValue(PrismVertex1,PrismEdge0[i],Rvalue);
1967  Rvalue=(*m_Rtet)(TetVertex1,TetEdge1[i]);
1968  Rmodprism->SetValue(PrismVertex1,PrismEdge1[i],Rvalue);
1969  Rvalue=(*m_Rtet)(TetVertex1,TetEdge4[i]);
1970  Rmodprism->SetValue(PrismVertex1,PrismEdge5[i],Rvalue);
1971 
1972  //transposed values
1973  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1974  RTmodprism->SetValue(PrismEdge0[i],PrismVertex1,RTvalue);
1975  RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex1);
1976  RTmodprism->SetValue(PrismEdge1[i],PrismVertex1,RTvalue);
1977  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex1);
1978  RTmodprism->SetValue(PrismEdge5[i],PrismVertex1,RTvalue);
1979  }
1980 
1981  //vertex 2 edge 1 2 & 6
1982  for(i=0; i< PrismEdge2.num_elements(); ++i)
1983  {
1984  Rvalue=(*m_Rtet)(TetVertex2,TetEdge1[i]);
1985  Rmodprism->SetValue(PrismVertex2,PrismEdge1[i],Rvalue);
1986  Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1987  Rmodprism->SetValue(PrismVertex2,PrismEdge2[i],Rvalue);
1988  Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
1989  Rmodprism->SetValue(PrismVertex2,PrismEdge6[i],Rvalue);
1990 
1991  //transposed values
1992  RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex2);
1993  RTmodprism->SetValue(PrismEdge1[i],PrismVertex2,RTvalue);
1994  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1995  RTmodprism->SetValue(PrismEdge2[i],PrismVertex2,RTvalue);
1996  RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
1997  RTmodprism->SetValue(PrismEdge6[i],PrismVertex2,RTvalue);
1998  }
1999 
2000  //vertex 3 edge 3 2 & 7
2001  for(i=0; i< PrismEdge3.num_elements(); ++i)
2002  {
2003  Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2004  Rmodprism->SetValue(PrismVertex3,PrismEdge3[i],Rvalue);
2005  Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
2006  Rmodprism->SetValue(PrismVertex3,PrismEdge2[i],Rvalue);
2007  Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
2008  Rmodprism->SetValue(PrismVertex3,PrismEdge7[i],Rvalue);
2009 
2010  //transposed values
2011  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2012  RTmodprism->SetValue(PrismEdge3[i],PrismVertex3,RTvalue);
2013  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
2014  RTmodprism->SetValue(PrismEdge2[i],PrismVertex3,RTvalue);
2015  RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2016  RTmodprism->SetValue(PrismEdge7[i],PrismVertex3,RTvalue);
2017  }
2018 
2019  //vertex 4 edge 4 5 & 8
2020  for(i=0; i< PrismEdge4.num_elements(); ++i)
2021  {
2022  Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2023  Rmodprism->SetValue(PrismVertex4,PrismEdge4[i],Rvalue);
2024  Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2025  Rmodprism->SetValue(PrismVertex4,PrismEdge5[i],Rvalue);
2026  Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
2027  Rmodprism->SetValue(PrismVertex4,PrismEdge8[i],Rvalue);
2028 
2029  //transposed values
2030  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2031  RTmodprism->SetValue(PrismEdge4[i],PrismVertex4,RTvalue);
2032  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2033  RTmodprism->SetValue(PrismEdge5[i],PrismVertex4,RTvalue);
2034  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
2035  RTmodprism->SetValue(PrismEdge8[i],PrismVertex4,RTvalue);
2036  }
2037 
2038  //vertex 5 edge 6 7 & 8
2039  for(i=0; i< PrismEdge5.num_elements(); ++i)
2040  {
2041  Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2042  Rmodprism->SetValue(PrismVertex5,PrismEdge6[i],Rvalue);
2043  Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2044  Rmodprism->SetValue(PrismVertex5,PrismEdge7[i],Rvalue);
2045  Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2046  Rmodprism->SetValue(PrismVertex5,PrismEdge8[i],Rvalue);
2047 
2048  //transposed values
2049  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2050  RTmodprism->SetValue(PrismEdge6[i],PrismVertex5,RTvalue);
2051  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2052  RTmodprism->SetValue(PrismEdge7[i],PrismVertex5,RTvalue);
2053  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2054  RTmodprism->SetValue(PrismEdge8[i],PrismVertex5,RTvalue);
2055  }
2056 
2057  // face 1 vertices 0 1 4
2058  for(i=0; i< PrismFace1.num_elements(); ++i)
2059  {
2060  Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2061  Rmodprism->SetValue(PrismVertex0,PrismFace1[i],Rvalue);
2062  Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2063  Rmodprism->SetValue(PrismVertex1,PrismFace1[i],Rvalue);
2064  Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2065  Rmodprism->SetValue(PrismVertex4,PrismFace1[i],Rvalue);
2066 
2067  //transposed values
2068  RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2069  RTmodprism->SetValue(PrismFace1[i],PrismVertex0,RTvalue);
2070  RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2071  RTmodprism->SetValue(PrismFace1[i],PrismVertex1,RTvalue);
2072  RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2073  RTmodprism->SetValue(PrismFace1[i],PrismVertex4,RTvalue);
2074  }
2075 
2076  // face 3 vertices 2, 3 & 5
2077  for(i=0; i< PrismFace3.num_elements(); ++i)
2078  {
2079  Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2080  Rmodprism->SetValue(PrismVertex2,PrismFace3[i],Rvalue);
2081  Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2082  Rmodprism->SetValue(PrismVertex3,PrismFace3[i],Rvalue);
2083  Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2084  Rmodprism->SetValue(PrismVertex5,PrismFace3[i],Rvalue);
2085 
2086  //transposed values
2087  RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2088  RTmodprism->SetValue(PrismFace3[i],PrismVertex2,RTvalue);
2089  RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2090  RTmodprism->SetValue(PrismFace3[i],PrismVertex3,RTvalue);
2091  RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2092  RTmodprism->SetValue(PrismFace3[i],PrismVertex5,RTvalue);
2093  }
2094 
2095  // Face 1 edge 0 4 5
2096  for(i=0; i< PrismFace1.num_elements(); ++i)
2097  {
2098  for(j=0; j<PrismEdge0.num_elements(); ++j)
2099  {
2100  Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2101  Rmodprism->SetValue(PrismEdge0[j],PrismFace1[i],Rvalue);
2102  Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2103  Rmodprism->SetValue(PrismEdge4[j],PrismFace1[i],Rvalue);
2104  Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2105  Rmodprism->SetValue(PrismEdge5[j],PrismFace1[i],Rvalue);
2106 
2107  //transposed values
2108  RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2109  RTmodprism->SetValue(PrismFace1[i],PrismEdge0[j],RTvalue);
2110  RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2111  RTmodprism->SetValue(PrismFace1[i],PrismEdge4[j],RTvalue);
2112  RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2113  RTmodprism->SetValue(PrismFace1[i],PrismEdge5[j],RTvalue);
2114  }
2115  }
2116 
2117  // Face 3 edge 2 6 7
2118  for(i=0; i< PrismFace3.num_elements(); ++i)
2119  {
2120  for(j=0; j<PrismEdge2.num_elements(); ++j)
2121  {
2122  Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2123  Rmodprism->SetValue(PrismEdge2[j],PrismFace3[i],Rvalue);
2124  Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2125  Rmodprism->SetValue(PrismEdge6[j],PrismFace3[i],Rvalue);
2126  Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2127  Rmodprism->SetValue(PrismEdge7[j],PrismFace3[i],Rvalue);
2128 
2129  RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2130  RTmodprism->SetValue(PrismFace3[i],PrismEdge2[j],RTvalue);
2131  RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2132  RTmodprism->SetValue(PrismFace3[i],PrismEdge6[j],RTvalue);
2133  RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2134  RTmodprism->SetValue(PrismFace3[i],PrismEdge7[j],RTvalue);
2135  }
2136  }
2137  }
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
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 913 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eDIAGONAL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eTetrahedron, m_InvRBlk, m_InvRTBlk, m_linsys, m_locToGloMap, m_RBlk, m_Rhex, m_Rinvhex, m_Rinvprism, m_Rinvtet, m_Rprism, m_RTBlk, m_Rtet, m_RThex, m_RTinvhex, m_RTinvprism, m_RTinvtet, m_RTprism, and m_RTtet.

Referenced by v_InitObject().

914  {
915  boost::shared_ptr<MultiRegions::ExpList>
916  expList=((m_linsys.lock())->GetLocMat()).lock();
917  StdRegions::StdExpansionSharedPtr locExpansion;
918 
919  int n, nel;
920 
921  const Array<OneD,const unsigned int>& nbdry_size
922  = m_locToGloMap->GetNumLocalBndCoeffsPerPatch();
923 
924  int n_exp=expList->GetNumElmts();
925 
926  //maps for different element types
927  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
928  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
929  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransmatrixmap;
930  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransposedtransmatrixmap;
931 
932  //Transformation matrix map
933  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
934  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
935  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
936 
937  //Transposed transformation matrix map
938  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
939  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
940  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
941 
942  //Inverse transfomation map
943  invtransmatrixmap[LibUtilities::eTetrahedron]=m_Rinvtet;
944  invtransmatrixmap[LibUtilities::ePrism]=m_Rinvprism;
945  invtransmatrixmap[LibUtilities::eHexahedron]=m_Rinvhex;
946 
947  //Inverse transposed transformation map
948  invtransposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTinvtet;
949  invtransposedtransmatrixmap[LibUtilities::ePrism]=m_RTinvprism;
950  invtransposedtransmatrixmap[LibUtilities::eHexahedron]=m_RTinvhex;
951 
952  MatrixStorage blkmatStorage = eDIAGONAL;
953 
954  //Variants of R matrices required for low energy preconditioning
955  m_RBlk = MemoryManager<DNekScalBlkMat>
956  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
957  m_RTBlk = MemoryManager<DNekScalBlkMat>
958  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
959  m_InvRBlk = MemoryManager<DNekScalBlkMat>
960  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
961  m_InvRTBlk = MemoryManager<DNekScalBlkMat>
962  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
963 
964  for(n=0; n < n_exp; ++n)
965  {
966  nel = expList->GetOffset_Elmt_Id(n);
967 
968  locExpansion = expList->GetExp(nel);
969  LibUtilities::ShapeType eType=locExpansion->DetShapeType();
970 
971  //Block R matrix
972  m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
973 
974  //Block RT matrix
975  m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
976 
977  //Block inverse R matrix
978  m_InvRBlk->SetBlock(n,n, invtransmatrixmap[eType]);
979 
980  //Block inverse RT matrix
981  m_InvRTBlk->SetBlock(n,n, invtransposedtransmatrixmap[eType]);
982  }
983  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::MultiRegions::PreconditionerLowEnergy::m_linsys
const boost::weak_ptr< GlobalLinSys > m_linsys
Definition: PreconditionerLowEnergy.h:79
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
void Nektar::MultiRegions::PreconditionerLowEnergy::SetUpReferenceElements ( void  )
private

Sets up the reference elements needed by the preconditioner.

Sets up reference elements which are used to preconditioning the corresponding matrices. Currently we support tetrahedral, prismatic and hexahedral elements

Definition at line 1600 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), CreateRefHexGeom(), CreateRefPrismGeom(), CreateRefTetGeom(), Nektar::eFULL, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eGaussRadauMAlpha2Beta0, Nektar::LibUtilities::eHexahedron, Nektar::StdRegions::eMass, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, Nektar::StdRegions::ePreconRT, Nektar::StdRegions::ePreconRTMass, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eTetrahedron, Nektar::MultiRegions::GlobalMatrixKey::GetConstFactors(), Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), Nektar::MultiRegions::GlobalMatrixKey::GetNVarCoeffs(), Nektar::MultiRegions::GlobalMatrixKey::GetVarCoeffs(), m_linsys, m_Rhex, m_Rinvhex, m_Rinvprism, m_Rinvtet, m_Rprism, m_Rtet, m_RThex, m_RTinvhex, m_RTinvprism, m_RTinvtet, m_RTprism, m_RTtet, and ModifyPrismTransformationMatrix().

Referenced by v_InitObject().

1601  {
1602  int cnt,i,j;
1603  boost::shared_ptr<MultiRegions::ExpList>
1604  expList=((m_linsys.lock())->GetLocMat()).lock();
1605  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
1606  StdRegions::VarCoeffMap vVarCoeffMap;
1607  StdRegions::StdExpansionSharedPtr locExpansion;
1608  locExpansion = expList->GetExp(0);
1609 
1610  DNekScalBlkMatSharedPtr RtetBlk, RprismBlk;
1611  DNekScalBlkMatSharedPtr RTtetBlk, RTprismBlk;
1612 
1613  DNekScalMatSharedPtr Rprismoriginal;
1614  DNekScalMatSharedPtr RTprismoriginal;
1615  DNekMatSharedPtr Rtettmp, RTtettmp, Rhextmp, RThextmp, Rprismtmp, RTprismtmp ;
1616 
1617  /*
1618  * Set up a Tetrahral & prismatic element which comprises
1619  * equilateral triangles as all faces for the tet and the end faces
1620  * for the prism. Using these elements a new expansion is created
1621  * (which is the same as the expansion specified in the input
1622  * file).
1623  */
1624  SpatialDomains::TetGeomSharedPtr tetgeom=CreateRefTetGeom();
1625  SpatialDomains::PrismGeomSharedPtr prismgeom=CreateRefPrismGeom();
1626  SpatialDomains::HexGeomSharedPtr hexgeom=CreateRefHexGeom();
1627 
1628  //Expansion as specified in the input file - here we need to alter
1629  //this so we can read in different exapansions for different element
1630  //types
1631  int nummodes=locExpansion->GetBasisNumModes(0);
1632 
1633  //Bases for Tetrahedral element
1634  const LibUtilities::BasisKey TetBa(
1635  LibUtilities::eModified_A, nummodes,
1636  LibUtilities::PointsKey(nummodes+1,LibUtilities::eGaussLobattoLegendre));
1637  const LibUtilities::BasisKey TetBb(
1638  LibUtilities::eModified_B, nummodes,
1639  LibUtilities::PointsKey(nummodes,LibUtilities::eGaussRadauMAlpha1Beta0));
1640  const LibUtilities::BasisKey TetBc(
1641  LibUtilities::eModified_C, nummodes,
1642  LibUtilities::PointsKey(nummodes,LibUtilities::eGaussRadauMAlpha2Beta0));
1643 
1644  //Create reference tetrahedral expansion
1645  LocalRegions::TetExpSharedPtr TetExp;
1646 
1647  TetExp = MemoryManager<LocalRegions::TetExp>
1648  ::AllocateSharedPtr(TetBa,TetBb,TetBc,
1649  tetgeom);
1650 
1651  //Bases for prismatic element
1652  const LibUtilities::BasisKey PrismBa(
1653  LibUtilities::eModified_A, nummodes,
1654  LibUtilities::PointsKey(nummodes+1,LibUtilities::eGaussLobattoLegendre));
1655  const LibUtilities::BasisKey PrismBb(
1656  LibUtilities::eModified_A, nummodes,
1657  LibUtilities::PointsKey(nummodes+1,LibUtilities::eGaussLobattoLegendre));
1658  const LibUtilities::BasisKey PrismBc(
1659  LibUtilities::eModified_B, nummodes,
1660  LibUtilities::PointsKey(nummodes,LibUtilities::eGaussRadauMAlpha1Beta0));
1661 
1662  //Create reference prismatic expansion
1663  LocalRegions::PrismExpSharedPtr PrismExp;
1664 
1665  PrismExp = MemoryManager<LocalRegions::PrismExp>
1666  ::AllocateSharedPtr(PrismBa,PrismBb,PrismBc,
1667  prismgeom);
1668 
1669  //Bases for prismatic element
1670  const LibUtilities::BasisKey HexBa(
1671  LibUtilities::eModified_A, nummodes,
1672  LibUtilities::PointsKey(nummodes+1,LibUtilities::eGaussLobattoLegendre));
1673  const LibUtilities::BasisKey HexBb(
1674  LibUtilities::eModified_A, nummodes,
1675  LibUtilities::PointsKey(nummodes+1,LibUtilities::eGaussLobattoLegendre));
1676  const LibUtilities::BasisKey HexBc(
1677  LibUtilities::eModified_A, nummodes,
1678  LibUtilities::PointsKey(nummodes+1,LibUtilities::eGaussLobattoLegendre));
1679 
1680  //Create reference prismatic expansion
1681  LocalRegions::HexExpSharedPtr HexExp;
1682 
1683  HexExp = MemoryManager<LocalRegions::HexExp>
1684  ::AllocateSharedPtr(HexBa,HexBb,HexBc,
1685  hexgeom);
1686 
1687 
1688  // retrieve variable coefficient
1689  if(m_linSysKey.GetNVarCoeffs() > 0)
1690  {
1691  StdRegions::VarCoeffMap::const_iterator x;
1692  cnt = expList->GetPhys_Offset(0);
1693  for (x = m_linSysKey.GetVarCoeffs().begin();
1694  x != m_linSysKey.GetVarCoeffs().end(); ++x)
1695  {
1696  vVarCoeffMap[x->first] = x->second + cnt;
1697  }
1698  }
1699 
1700  StdRegions::MatrixType PreconR,PreconRT;
1701 
1702  if(m_linSysKey.GetMatrixType() == StdRegions::eMass)
1703  {
1704  PreconR = StdRegions::ePreconRMass;
1705  PreconRT = StdRegions::ePreconRTMass;
1706  }
1707  else
1708  {
1709  PreconR = StdRegions::ePreconR;
1710  PreconRT = StdRegions::ePreconRT;
1711  }
1712 
1713 
1714 
1715  /*
1716  * Matrix keys - for each element type there are two matrix keys
1717  * corresponding to the transformation matrix R and its transpose
1718  */
1719 
1720 
1721  //Matrix keys for tetrahedral element transformation matrix
1722  LocalRegions::MatrixKey TetR
1723  (PreconR, LibUtilities::eTetrahedron,
1724  *TetExp, m_linSysKey.GetConstFactors(),
1725  vVarCoeffMap);
1726 
1727  //Matrix keys for tetrahedral transposed transformation matrix
1728  LocalRegions::MatrixKey TetRT
1729  (PreconRT, LibUtilities::eTetrahedron,
1730  *TetExp, m_linSysKey.GetConstFactors(),
1731  vVarCoeffMap);
1732 
1733  //Matrix keys for prismatic element transformation matrix
1734  LocalRegions::MatrixKey PrismR
1735  (PreconR, LibUtilities::ePrism,
1736  *PrismExp, m_linSysKey.GetConstFactors(),
1737  vVarCoeffMap);
1738 
1739  //Matrix keys for prismatic element transposed transformation matrix
1740  LocalRegions::MatrixKey PrismRT
1741  (PreconRT, LibUtilities::ePrism,
1742  *PrismExp, m_linSysKey.GetConstFactors(),
1743  vVarCoeffMap);
1744 
1745  //Matrix keys for hexahedral element transformation matrix
1746  LocalRegions::MatrixKey HexR
1747  (PreconR, LibUtilities::eHexahedron,
1748  *HexExp, m_linSysKey.GetConstFactors(),
1749  vVarCoeffMap);
1750 
1751  //Matrix keys for hexahedral element transposed transformation
1752  //matrix
1753  LocalRegions::MatrixKey HexRT
1754  (PreconRT, LibUtilities::eHexahedron,
1755  *HexExp, m_linSysKey.GetConstFactors(),
1756  vVarCoeffMap);
1757 
1758  /*
1759  * Create transformation matrices for the tetrahedral element
1760  */
1761 
1762  //Get tetrahedral transformation matrix
1763  m_Rtet = TetExp->GetLocMatrix(TetR);
1764 
1765  //Get tetrahedral transposed transformation matrix
1766  m_RTtet = TetExp->GetLocMatrix(TetRT);
1767 
1768  // Using the transformation matrix and the inverse transformation
1769  // matrix create the inverse matrices
1770  Rtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1771 
1772  //Inverse transposed transformation matrix
1773  RTtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1774  RTtettmp->Transpose();
1775 
1776  m_Rinvtet = MemoryManager<DNekScalMat>
1777  ::AllocateSharedPtr(1.0,Rtettmp);
1778  m_RTinvtet = MemoryManager<DNekScalMat>
1779  ::AllocateSharedPtr(1.0,RTtettmp);
1780 
1781  /*
1782  * Create transformation matrices for the hexahedral element
1783  */
1784 
1785  //Get hexahedral transformation matrix
1786  m_Rhex = HexExp->GetLocMatrix(HexR);
1787  //Get hexahedral transposed transformation matrix
1788  m_RThex = HexExp->GetLocMatrix(HexRT);
1789 
1790  // Using the transformation matrix and the inverse transformation
1791  // matrix create the inverse matrices
1792  Rhextmp=HexExp->BuildInverseTransformationMatrix(m_Rhex);
1793  //Inverse transposed transformation matrix
1794  RThextmp=HexExp->BuildInverseTransformationMatrix(m_Rhex);
1795  RThextmp->Transpose();
1796 
1797  m_Rinvhex = MemoryManager<DNekScalMat>
1798  ::AllocateSharedPtr(1.0,Rhextmp);
1799  m_RTinvhex = MemoryManager<DNekScalMat>
1800  ::AllocateSharedPtr(1.0,RThextmp);
1801 
1802  /*
1803  * Create transformation matrices for the prismatic element
1804  */
1805 
1806  //Get prism transformation matrix
1807  Rprismoriginal = PrismExp->GetLocMatrix(PrismR);
1808  //Get prism transposed transformation matrix
1809  RTprismoriginal = PrismExp->GetLocMatrix(PrismRT);
1810 
1811  unsigned int nRows=Rprismoriginal->GetRows();
1812  NekDouble zero=0.0;
1813  DNekMatSharedPtr Rtmpprism = MemoryManager<DNekMat>::
1814  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1815  DNekMatSharedPtr RTtmpprism = MemoryManager<DNekMat>::
1816  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1817  NekDouble Rvalue, RTvalue;
1818 
1819  //Copy values from the prism transformation matrix
1820  for(i=0; i<nRows; ++i)
1821  {
1822  for(j=0; j<nRows; ++j)
1823  {
1824  Rvalue=(*Rprismoriginal)(i,j);
1825  RTvalue=(*RTprismoriginal)(i,j);
1826  Rtmpprism->SetValue(i,j,Rvalue);
1827  RTtmpprism->SetValue(i,j,RTvalue);
1828  }
1829  }
1830 
1831  //Replace triangular faces and edges of the prims transformation
1832  //matrix with the corresponding values of the tetrahedral
1833  //transformation matrix.
1834  ModifyPrismTransformationMatrix(TetExp,PrismExp,Rtmpprism,RTtmpprism);
1835 
1836  m_Rprism = MemoryManager<DNekScalMat>
1837  ::AllocateSharedPtr(1.0,Rtmpprism);
1838 
1839  m_RTprism = MemoryManager<DNekScalMat>
1840  ::AllocateSharedPtr(1.0,RTtmpprism);
1841 
1842  //Inverse transformation matrix
1843  Rprismtmp=PrismExp->BuildInverseTransformationMatrix(m_Rprism);
1844 
1845  //Inverse transposed transformation matrix
1846  RTprismtmp=PrismExp->BuildInverseTransformationMatrix(m_Rprism);
1847  RTprismtmp->Transpose();
1848 
1849  m_Rinvprism = MemoryManager<DNekScalMat>
1850  ::AllocateSharedPtr(1.0,Rprismtmp);
1851 
1852  m_RTinvprism = MemoryManager<DNekScalMat>
1853  ::AllocateSharedPtr(1.0,RTprismtmp);
1854  }
Nektar::LibUtilities::eModified_C
Principle Modified Functions .
Definition: BasisType.h:51
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::LibUtilities::eModified_A
Principle Modified Functions .
Definition: BasisType.h:49
Nektar::LocalRegions::HexExpSharedPtr
boost::shared_ptr< HexExp > HexExpSharedPtr
Definition: HexExp.h:57
Nektar::SpatialDomains::HexGeomSharedPtr
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Nektar::LocalRegions::TetExpSharedPtr
boost::shared_ptr< TetExp > TetExpSharedPtr
Definition: TetExp.h:223
Nektar::DNekScalMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:73
Nektar::LibUtilities::eGaussRadauMAlpha1Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefHexGeom
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1511
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:226
Nektar::LibUtilities::eModified_B
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::DNekScalBlkMatSharedPtr
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LocalRegions::PrismExpSharedPtr
boost::shared_ptr< PrismExp > PrismExpSharedPtr
Definition: PrismExp.h:217
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefTetGeom
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom(void)
Sets up the reference tretrahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1419
Nektar::MultiRegions::PreconditionerLowEnergy::ModifyPrismTransformationMatrix
void ModifyPrismTransformationMatrix(LocalRegions::TetExpSharedPtr TetExp, LocalRegions::PrismExpSharedPtr PrismExp, DNekMatSharedPtr Rmodprism, DNekMatSharedPtr RTmodprism)
Modify the prism transformation matrix to align with the tetrahedral modes.
Definition: PreconditionerLowEnergy.cpp:1868
Nektar::SpatialDomains::PrismGeomSharedPtr
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
Nektar::LibUtilities::eGaussRadauMAlpha2Beta0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Nektar::SpatialDomains::TetGeomSharedPtr
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
Nektar::MultiRegions::PreconditionerLowEnergy::CreateRefPrismGeom
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1329
Nektar::MultiRegions::PreconditionerLowEnergy::m_linsys
const boost::weak_ptr< GlobalLinSys > m_linsys
Definition: PreconditionerLowEnergy.h:79
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
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 edges and face sizes

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 119 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::as(), Vmath::Assmb(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eTetrahedron, Gs::Gather(), Nektar::StdRegions::StdExpansion::GetEdgeNcoeffs(), Nektar::StdRegions::StdExpansion::GetFaceIntNcoeffs(), Gs::gs_add, Gs::Init(), m_BlkMat, Nektar::MultiRegions::Preconditioner::m_comm, m_linsys, m_locToGloMap, m_RBlk, m_Rhex, m_Rprism, m_RTBlk, m_Rtet, m_RThex, m_RTprism, m_RTtet, and Nektar::LibUtilities::ReduceMax.

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

References ASSERTL1, Nektar::eWrapper, Vmath::Gathr(), m_InvRBlk, m_locToGloMap, m_locToGloSignMult, m_map, and Vmath::Vcopy().

1131  {
1132  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1133  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1134  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1135  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1136 
1137  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1138  "Input array is greater than the nGlobHomBndDofs");
1139  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1140  "Output array is greater than the nGlobHomBndDofs");
1141 
1142  //vectors of length number of non-dirichlet boundary dofs
1143  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1144  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1145  eWrapper);
1146  //Block inverse transformation matrix
1147  DNekScalBlkMat &invR = *m_InvRBlk;
1148 
1149  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1150  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1151  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1152 
1153  // Allocated array of size number of global boundary dofs and copy
1154  // the input array to the tmp array offset by Dirichlet boundary
1155  // conditions.
1156  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1157  Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1158 
1159  //Global boundary dofs (with zeroed dirichlet values) to local boundary dofs
1160  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1161 
1162  //Multiply by block inverse transformation matrix
1163  F_LocBnd=invR*F_LocBnd;
1164 
1165  //Assemble local boundary to global non-dirichlet boundary
1166  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
1167 
1168  }
Vmath::Gathr
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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 1173 of file PreconditionerLowEnergy.cpp.

References ASSERTL1, Nektar::eWrapper, m_InvRTBlk, m_locToGloMap, m_map, m_multiplicity, and Vmath::Vmul().

1176  {
1177  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1178  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1179  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1180  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1181 
1182  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1183  "Input array is greater than the nGlobHomBndDofs");
1184  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1185  "Output array is greater than the nGlobHomBndDofs");
1186 
1187  //vectors of length number of non-dirichlet boundary dofs
1188  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1189  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1190  eWrapper);
1191  //Block inverse transformation matrix
1192  DNekScalBlkMat &invRT = *m_InvRTBlk;
1193 
1194  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1195  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1196  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1197 
1198  m_locToGloMap->GlobalToLocalBnd(pInput,pLocal, nDirBndDofs);
1199 
1200  //Multiply by the block transposed transformation matrix
1201  F_LocBnd=invRT*F_LocBnd;
1202 
1203  m_locToGloMap->AssembleBnd(pLocal,pOutput, nDirBndDofs);
1204 
1205  Vmath::Vmul(nGlobHomBndDofs,pOutput,1,m_multiplicity,1,pOutput,1);
1206  }
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
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 892 of file PreconditionerLowEnergy.cpp.

References Nektar::eWrapper, and m_locToGloMap.

895  {
896  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
897  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
898  int nNonDir = nGlobal-nDir;
899  DNekBlkMat &M = (*m_BlkMat);
900 
901  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
902  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
903 
904  z = M * r;
905  }
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 1086 of file PreconditionerLowEnergy.cpp.

References ASSERTL1, Vmath::Assmb(), Nektar::eWrapper, m_locToGloMap, m_locToGloSignMult, m_map, m_RTBlk, and Vmath::Vcopy().

1088  {
1089  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1090  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1091  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1092  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1093 
1094  ASSERTL1(pInOut.num_elements() >= nGlobBndDofs,
1095  "Output array is greater than the nGlobBndDofs");
1096 
1097  //Block transposed transformation matrix
1098  DNekScalBlkMat &RT = *m_RTBlk;
1099 
1100  NekVector<NekDouble> V_GlobHomBnd(nGlobHomBndDofs,pInOut+nDirBndDofs,
1101  eWrapper);
1102 
1103  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1104  NekVector<NekDouble> V_LocBnd(nLocBndDofs,pLocal,eWrapper);
1105  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1106  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1107 
1108  //Global boundary (less dirichlet) to local boundary
1109  m_locToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
1110 
1111  //Multiply by the block transposed transformation matrix
1112  V_LocBnd=RT*V_LocBnd;
1113 
1114 
1115  //Assemble local boundary to global boundary
1116  Vmath::Assmb(nLocBndDofs, m_locToGloSignMult.get(),pLocal.get(), m_map.get(), tmp.get());
1117 
1118  //Universal assemble across processors
1119  m_locToGloMap->UniversalAssembleBnd(tmp);
1120 
1121  //copy non-dirichlet boundary values
1122  Vmath::Vcopy(nGlobBndDofs-nDirBndDofs, tmp.get() + nDirBndDofs, 1, pInOut.get() + nDirBndDofs, 1);
1123  }
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:695
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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 994 of file PreconditionerLowEnergy.cpp.

References Nektar::eWrapper, Vmath::Gathr(), m_locToGloMap, m_locToGloSignMult, m_map, m_RBlk, and Vmath::Vcopy().

997  {
998  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
999  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1000  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1001  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1002 
1003  //Non-dirichlet boundary dofs
1004  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pInOut+offset,
1005  eWrapper);
1006 
1007  //Block transformation matrix
1008  DNekScalBlkMat &R = *m_RBlk;
1009 
1010  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1011  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1012  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1013 
1014  //Not actually needed but we should only work with the Global boundary dofs
1015  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1016  Vmath::Vcopy(nGlobBndDofs, pInOut.get(), 1, tmp.get(), 1);
1017 
1018  //Global boundary (with dirichlet values) to local boundary with multiplicity
1019  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1020 
1021  //Multiply by the block transformation matrix
1022  F_LocBnd=R*F_LocBnd;
1023 
1024  //Assemble local boundary to global non-dirichlet Dofs
1025  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd, nDirBndDofs);
1026  }
Vmath::Gathr
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformToLowEnergy ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
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 1035 of file PreconditionerLowEnergy.cpp.

References ASSERTL1, Nektar::eWrapper, Vmath::Gathr(), m_locToGloMap, m_locToGloSignMult, m_map, m_RBlk, and Vmath::Vcopy().

1038  {
1039  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1040  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1041  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1042  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1043 
1044  //Input/output vectors should be length nGlobHomBndDofs
1045  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1046  "Input array is greater than the nGlobHomBndDofs");
1047  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1048  "Output array is greater than the nGlobHomBndDofs");
1049 
1050  //vectors of length number of non-dirichlet boundary dofs
1051  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1052  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1053  eWrapper);
1054  //Block transformation matrix
1055  DNekScalBlkMat &R = *m_RBlk;
1056 
1057  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1058  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1059  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1060 
1061  // Allocated array of size number of global boundary dofs and copy
1062  // the input array to the tmp array offset by Dirichlet boundary
1063  // conditions.
1064  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1065  Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1066 
1067  //Global boundary dofs (with zeroed dirichlet values) to local boundary dofs
1068  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1069 
1070  //Multiply by the block transformation matrix
1071  F_LocBnd=R*F_LocBnd;
1072 
1073  //Assemble local boundary to global non-dirichlet boundary
1074  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
1075  }
Vmath::Gathr
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::MultiRegions::PreconditionerLowEnergy::v_InitObject ( )
privatevirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 73 of file PreconditionerLowEnergy.cpp.

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

74  {
75  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
76  ASSERTL0(solvertype == eIterativeStaticCond ||
77  solvertype == ePETScStaticCond, "Solver type not valid");
78 
79  boost::shared_ptr<MultiRegions::ExpList>
80  expList=((m_linsys.lock())->GetLocMat()).lock();
81 
82  StdRegions::StdExpansionSharedPtr locExpansion;
83 
84  locExpansion = expList->GetExp(0);
85 
86  int nDim = locExpansion->GetShapeDimension();
87 
88  ASSERTL0(nDim==3,
89  "Preconditioner type only valid in 3D");
90 
91  //Sets up reference element and builds transformation matrix
92  SetUpReferenceElements();
93 
94  //Set up block transformation matrix
95  SetupBlockTransformationMatrix();
96 
97  //Sets up multiplicity map for transformation from global to local
98  CreateMultiplicityMap();
99  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::MultiRegions::PreconditionerLowEnergy::SetUpReferenceElements
void SetUpReferenceElements(void)
Sets up the reference elements needed by the preconditioner.
Definition: PreconditionerLowEnergy.cpp:1600
Nektar::MultiRegions::PreconditionerLowEnergy::m_linsys
const boost::weak_ptr< GlobalLinSys > m_linsys
Definition: PreconditionerLowEnergy.h:79
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::v_TransformedSchurCompl ( int  offset,
const boost::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 1217 of file PreconditionerLowEnergy.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eTetrahedron, m_linsys, m_Rhex, m_Rprism, m_Rtet, m_RThex, m_RTprism, and m_RTtet.

1220  {
1221  boost::shared_ptr<MultiRegions::ExpList>
1222  expList=((m_linsys.lock())->GetLocMat()).lock();
1223 
1224  StdRegions::StdExpansionSharedPtr locExpansion;
1225  locExpansion = expList->GetExp(offset);
1226  unsigned int nbnd=locExpansion->NumBndryCoeffs();
1227 
1228  //This is the SC elemental matrix in the orginal basis (S1)
1229  DNekScalMatSharedPtr pS1=loc_mat;
1230 
1231  //Transformation matrices
1232  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
1233  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
1234  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
1235  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
1236  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
1237  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
1238  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
1239  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
1240 
1241  DNekScalMat &S1 = (*pS1);
1242 
1243  MatrixStorage storage = eFULL;
1244 
1245  DNekMatSharedPtr pS2 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1246  DNekMatSharedPtr pRS1 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1247 
1248  LibUtilities::ShapeType eType=
1249  (expList->GetExp(offset))->DetShapeType();
1250 
1251  //transformation matrices
1252  DNekScalMat &R = (*(transmatrixmap[eType]));
1253  DNekScalMat &RT = (*(transposedtransmatrixmap[eType]));
1254 
1255  //create low energy matrix
1256  DNekMat &RS1 = (*pRS1);
1257  DNekMat &S2 = (*pS2);
1258 
1259  //setup S2
1260  RS1=R*S1;
1261  S2=RS1*RT;
1262 
1263  DNekScalMatSharedPtr tmp_mat;
1264  tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
1265 
1266  return tmp_mat;
1267  }
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::DNekMatSharedPtr
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Nektar::DNekScalMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:73
Nektar::MultiRegions::PreconditionerLowEnergy::m_linsys
const boost::weak_ptr< GlobalLinSys > m_linsys
Definition: PreconditionerLowEnergy.h:79
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873

Member Data Documentation

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

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_BlkMat
protected

Definition at line 82 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_bnd_mat
protected

Definition at line 83 of file PreconditionerLowEnergy.h.

DNekScalBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_InvRBlk
protected

Definition at line 87 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and v_DoMultiplybyInverseTransformationMatrix().

DNekScalBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_InvRTBlk
protected

Definition at line 88 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and v_DoMultiplybyInverseTransposedTransformationMatrix().

const boost::weak_ptr<GlobalLinSys> Nektar::MultiRegions::PreconditionerLowEnergy::m_linsys
protected

Definition at line 79 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), v_InitObject(), and v_TransformedSchurCompl().

boost::shared_ptr<AssemblyMap> Nektar::MultiRegions::PreconditionerLowEnergy::m_locToGloMap
protected

Definition at line 80 of file PreconditionerLowEnergy.h.

Referenced by CreateMultiplicityMap(), SetupBlockTransformationMatrix(), v_BuildPreconditioner(), v_DoMultiplybyInverseTransformationMatrix(), v_DoMultiplybyInverseTransposedTransformationMatrix(), v_DoPreconditioner(), v_DoTransformFromLowEnergy(), v_DoTransformToLowEnergy(), and v_InitObject().

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

Definition at line 105 of file PreconditionerLowEnergy.h.

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

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

Definition at line 107 of file PreconditionerLowEnergy.h.

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

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

Definition at line 106 of file PreconditionerLowEnergy.h.

Referenced by CreateMultiplicityMap(), and v_DoMultiplybyInverseTransposedTransformationMatrix().

DNekScalBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RBlk
protected

Definition at line 85 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_BuildPreconditioner(), and v_DoTransformToLowEnergy().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_Rhex
protected

Definition at line 95 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_Rinvhex
protected

Definition at line 97 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_Rinvprism
protected

Definition at line 102 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_Rinvtet
protected

Definition at line 92 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_Rprism
protected

Definition at line 100 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

DNekScalBlkMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RTBlk
protected

Definition at line 86 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), v_BuildPreconditioner(), and v_DoTransformFromLowEnergy().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_Rtet
protected

Definition at line 90 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RThex
protected

Definition at line 96 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RTinvhex
protected

Definition at line 98 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RTinvprism
protected

Definition at line 103 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RTinvtet
protected

Definition at line 93 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RTprism
protected

Definition at line 101 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

DNekScalMatSharedPtr Nektar::MultiRegions::PreconditionerLowEnergy::m_RTtet
protected

Definition at line 91 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), SetUpReferenceElements(), v_BuildPreconditioner(), and v_TransformedSchurCompl().

Generated on Wed Mar 23 2016 20:14:35 for Nektar++ by   doxygen 1.8.8