Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 66 of file PreconditionerLowEnergy.cpp.

69  : Preconditioner(plinsys, pLocToGloMap),
70  m_linsys(plinsys),
71  m_locToGloMap(pLocToGloMap)
72  {
73  }
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(), and CellMLToNektar.cellml_metadata::p.

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
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
void Nektar::MultiRegions::PreconditionerLowEnergy::CreateMultiplicityMap ( void  )
private

Create the inverse multiplicity map.

Definition at line 1277 of file PreconditionerLowEnergy.cpp.

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

Referenced by v_InitObject().

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

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

Definition at line 1516 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

1517  {
1518  ////////////////////////////////
1519  // Set up Hexahedron vertices //
1520  ////////////////////////////////
1521 
1522  const int three=3;
1523 
1524  const int nVerts = 8;
1525  const double point[][3] = {
1526  {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1527  {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1528  };
1529 
1530  // Populate the list of verts
1531  SpatialDomains::PointGeomSharedPtr verts[8];
1532  for( int i = 0; i < nVerts; ++i ) {
1533  verts[i] = MemoryManager<SpatialDomains::PointGeom>
1534  ::AllocateSharedPtr(three, i, point[i][0],
1535  point[i][1], point[i][2]);
1536  }
1537 
1538  /////////////////////////////
1539  // Set up Hexahedron Edges //
1540  /////////////////////////////
1541 
1542  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1543  const int nEdges = 12;
1544  const int vertexConnectivity[][2] = {
1545  {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1546  {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1547  };
1548 
1549  // Populate the list of edges
1550  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1551  for( int i = 0; i < nEdges; ++i ) {
1552  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1553  for( int j = 0; j < 2; ++j ) {
1554  vertsArray[j] = verts[vertexConnectivity[i][j]];
1555  }
1556  edges[i] = MemoryManager<SpatialDomains::SegGeom>::
1557  AllocateSharedPtr( i, three, vertsArray);
1558  }
1559 
1560  /////////////////////////////
1561  // Set up Hexahedron faces //
1562  /////////////////////////////
1563 
1564  const int nFaces = 6;
1565  const int edgeConnectivity[][4] = {
1566  {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1567  {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1568  };
1569  const bool isEdgeFlipped[][4] = {
1570  {0,0,0,1}, {0,0,1,1}, {0,0,1,1},
1571  {0,0,1,1}, {0,0,1,1}, {0,0,0,1}
1572  };
1573 
1574  // Populate the list of faces
1575  SpatialDomains::QuadGeomSharedPtr faces[nFaces];
1576  for( int i = 0; i < nFaces; ++i ) {
1577  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1578  StdRegions::Orientation eorientArray[4];
1579  for( int j = 0; j < 4; ++j ) {
1580  edgeArray[j] = edges[edgeConnectivity[i][j]];
1581  eorientArray[j] = isEdgeFlipped[i][j] ?
1582  StdRegions::eBackwards : StdRegions::eForwards;
1583  }
1584  faces[i] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(i, edgeArray,
1585  eorientArray);
1586  }
1587 
1588  SpatialDomains::HexGeomSharedPtr geom =
1589  MemoryManager<SpatialDomains::HexGeom>::AllocateSharedPtr
1590  (faces);
1591 
1592  geom->SetOwnData();
1593 
1594  return geom;
1595  }
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
point
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
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 1334 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

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

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

Referenced by SetUpReferenceElements().

1425  {
1426  /////////////////////////////////
1427  // Set up Tetrahedron vertices //
1428  /////////////////////////////////
1429 
1430  int i,j;
1431  const int three=3;
1432  const int nVerts = 4;
1433  const double point[][3] = {
1434  {-1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1435  {1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1436  {0,2/sqrt(double(3)),-1/sqrt(double(6))},
1437  {0,0,3/sqrt(double(6))}};
1438 
1439  boost::shared_ptr<SpatialDomains::PointGeom> verts[4];
1440  for(i=0; i < nVerts; ++i)
1441  {
1442  verts[i] =
1443  MemoryManager<SpatialDomains::PointGeom>::
1444  AllocateSharedPtr
1445  ( three, i, point[i][0], point[i][1], point[i][2] );
1446  }
1447 
1448  //////////////////////////////
1449  // Set up Tetrahedron Edges //
1450  //////////////////////////////
1451 
1452  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1453  const int nEdges = 6;
1454  const int vertexConnectivity[][2] = {
1455  {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1456  };
1457 
1458  // Populate the list of edges
1459  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1460  for(i=0; i < nEdges; ++i)
1461  {
1462  boost::shared_ptr<SpatialDomains::PointGeom>
1463  vertsArray[2];
1464  for(j=0; j<2; ++j)
1465  {
1466  vertsArray[j] = verts[vertexConnectivity[i][j]];
1467  }
1468 
1469  edges[i] = MemoryManager<SpatialDomains::SegGeom>
1470  ::AllocateSharedPtr(i, three, vertsArray);
1471  }
1472 
1473  //////////////////////////////
1474  // Set up Tetrahedron faces //
1475  //////////////////////////////
1476 
1477  const int nFaces = 4;
1478  const int edgeConnectivity[][3] = {
1479  {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1480  };
1481  const bool isEdgeFlipped[][3] = {
1482  {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
1483  };
1484 
1485  // Populate the list of faces
1486  SpatialDomains::TriGeomSharedPtr faces[nFaces];
1487  for(i=0; i < nFaces; ++i)
1488  {
1489  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1490  StdRegions::Orientation eorientArray[3];
1491  for(j=0; j < 3; ++j)
1492  {
1493  edgeArray[j] = edges[edgeConnectivity[i][j]];
1494  eorientArray[j] = isEdgeFlipped[i][j] ?
1495  StdRegions::eBackwards : StdRegions::eForwards;
1496  }
1497 
1498 
1499  faces[i] = MemoryManager<SpatialDomains::TriGeom>
1500  ::AllocateSharedPtr(i, edgeArray, eorientArray);
1501  }
1502 
1503  SpatialDomains::TetGeomSharedPtr geom =
1504  MemoryManager<SpatialDomains::TetGeom>::AllocateSharedPtr
1505  (faces);
1506 
1507  geom->SetOwnData();
1508 
1509  return geom;
1510  }
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
point
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
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 1873 of file PreconditionerLowEnergy.cpp.

Referenced by SetUpReferenceElements().

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

919  {
920  boost::shared_ptr<MultiRegions::ExpList>
921  expList=((m_linsys.lock())->GetLocMat()).lock();
922  StdRegions::StdExpansionSharedPtr locExpansion;
923 
924  int n, nel;
925 
926  const Array<OneD,const unsigned int>& nbdry_size
927  = m_locToGloMap->GetNumLocalBndCoeffsPerPatch();
928 
929  int n_exp=expList->GetNumElmts();
930 
931  //maps for different element types
932  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
933  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
934  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransmatrixmap;
935  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransposedtransmatrixmap;
936 
937  //Transformation matrix map
938  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
939  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
940  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
941 
942  //Transposed transformation matrix map
943  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
944  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
945  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
946 
947  //Inverse transfomation map
948  invtransmatrixmap[LibUtilities::eTetrahedron]=m_Rinvtet;
949  invtransmatrixmap[LibUtilities::ePrism]=m_Rinvprism;
950  invtransmatrixmap[LibUtilities::eHexahedron]=m_Rinvhex;
951 
952  //Inverse transposed transformation map
953  invtransposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTinvtet;
954  invtransposedtransmatrixmap[LibUtilities::ePrism]=m_RTinvprism;
955  invtransposedtransmatrixmap[LibUtilities::eHexahedron]=m_RTinvhex;
956 
957  MatrixStorage blkmatStorage = eDIAGONAL;
958 
959  //Variants of R matrices required for low energy preconditioning
960  m_RBlk = MemoryManager<DNekScalBlkMat>
961  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
962  m_RTBlk = MemoryManager<DNekScalBlkMat>
963  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
964  m_InvRBlk = MemoryManager<DNekScalBlkMat>
965  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
966  m_InvRTBlk = MemoryManager<DNekScalBlkMat>
967  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
968 
969  for(n=0; n < n_exp; ++n)
970  {
971  nel = n;
972 
973  locExpansion = expList->GetExp(nel);
974  LibUtilities::ShapeType eType=locExpansion->DetShapeType();
975 
976  //Block R matrix
977  m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
978 
979  //Block RT matrix
980  m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
981 
982  //Block inverse R matrix
983  m_InvRBlk->SetBlock(n,n, invtransmatrixmap[eType]);
984 
985  //Block inverse RT matrix
986  m_InvRTBlk->SetBlock(n,n, invtransposedtransmatrixmap[eType]);
987  }
988  }
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:1878
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 1605 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().

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

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

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

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

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

1181  {
1182  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1183  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1184  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1185  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1186 
1187  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1188  "Input array is greater than the nGlobHomBndDofs");
1189  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1190  "Output array is greater than the nGlobHomBndDofs");
1191 
1192  //vectors of length number of non-dirichlet boundary dofs
1193  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1194  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1195  eWrapper);
1196  //Block inverse transformation matrix
1197  DNekScalBlkMat &invRT = *m_InvRTBlk;
1198 
1199  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1200  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1201  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1202 
1203  m_locToGloMap->GlobalToLocalBnd(pInput,pLocal, nDirBndDofs);
1204 
1205  //Multiply by the block transposed transformation matrix
1206  F_LocBnd=invRT*F_LocBnd;
1207 
1208  m_locToGloMap->AssembleBnd(pLocal,pOutput, nDirBndDofs);
1209 
1210  Vmath::Vmul(nGlobHomBndDofs,pOutput,1,m_multiplicity,1,pOutput,1);
1211  }
Nektar::DNekScalBlkMat
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:66
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
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:183
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 897 of file PreconditionerLowEnergy.cpp.

References Nektar::eWrapper, and m_locToGloMap.

900  {
901  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
902  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
903  int nNonDir = nGlobal-nDir;
904  DNekBlkMat &M = (*m_BlkMat);
905 
906  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
907  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
908 
909  z = M * r;
910  }
Nektar::DNekBlkMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:60
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 1091 of file PreconditionerLowEnergy.cpp.

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

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

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

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

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 75 of file PreconditionerLowEnergy.cpp.

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

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

1225  {
1226  boost::shared_ptr<MultiRegions::ExpList>
1227  expList=((m_linsys.lock())->GetLocMat()).lock();
1228 
1229  StdRegions::StdExpansionSharedPtr locExpansion;
1230  locExpansion = expList->GetExp(offset);
1231  unsigned int nbnd=locExpansion->NumBndryCoeffs();
1232 
1233  //This is the SC elemental matrix in the orginal basis (S1)
1234  DNekScalMatSharedPtr pS1=loc_mat;
1235 
1236  //Transformation matrices
1237  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
1238  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
1239  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
1240  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
1241  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
1242  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
1243  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
1244  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
1245 
1246  DNekScalMat &S1 = (*pS1);
1247 
1248  MatrixStorage storage = eFULL;
1249 
1250  DNekMatSharedPtr pS2 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1251  DNekMatSharedPtr pRS1 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1252 
1253  LibUtilities::ShapeType eType=
1254  (expList->GetExp(offset))->DetShapeType();
1255 
1256  //transformation matrices
1257  DNekScalMat &R = (*(transmatrixmap[eType]));
1258  DNekScalMat &RT = (*(transposedtransmatrixmap[eType]));
1259 
1260  //create low energy matrix
1261  DNekMat &RS1 = (*pRS1);
1262  DNekMat &S2 = (*pS2);
1263 
1264  //setup S2
1265  RS1=R*S1;
1266  S2=RS1*RT;
1267 
1268  DNekScalMatSharedPtr tmp_mat;
1269  tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
1270 
1271  return tmp_mat;
1272  }
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::DNekMat
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
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:1878
Nektar::DNekScalMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:72

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 Sat Nov 25 2017 16:19:02 for Nektar++ by   doxygen 1.8.8