Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Attributes | Private Member Functions | List of all members
Nektar::MultiRegions::PreconditionerLowEnergy Class Reference

#include <PreconditionerLowEnergy.h>

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

Public Member Functions

 PreconditionerLowEnergy (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~PreconditionerLowEnergy ()
 
- Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
 Preconditioner (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~Preconditioner ()
 
void DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)
 
void DoTransformToLowEnergy (Array< OneD, NekDouble > &pInOut, int offset)
 
void DoTransformToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoTransformFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 
void DoMultiplybyInverseTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoMultiplybyInverseTransposedTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void BuildPreconditioner ()
 
void InitObject ()
 
Array< OneD, NekDouble > AssembleStaticCondGlobalDiagonals ()
 Performs global assembly of diagonal entries to global Schur complement matrix. More...
 
const DNekScalBlkMatSharedPtr & GetBlockTransformedSchurCompl () const
 
const DNekScalBlkMatSharedPtr & GetBlockCMatrix () const
 
const DNekScalBlkMatSharedPtr & GetBlockInvDMatrix () const
 
const DNekScalBlkMatSharedPtr & GetBlockSchurCompl () const
 
const DNekScalBlkMatSharedPtr & GetBlockTransformationMatrix () const
 
const DNekScalBlkMatSharedPtr & GetBlockTransposedTransformationMatrix () const
 
DNekScalMatSharedPtr TransformedSchurCompl (int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
 

Static Public Member Functions

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

Static Public Attributes

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

Protected Attributes

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

Private Member Functions

void SetUpReferenceElements (void)
 Sets up the reference elements needed by the preconditioner. More...
 
void CreateMultiplicityMap (void)
 
void SetupBlockTransformationMatrix (void)
 
void ModifyPrismTransformationMatrix (LocalRegions::TetExpSharedPtr TetExp, LocalRegions::PrismExpSharedPtr PrismExp, DNekMatSharedPtr Rmodprism, DNekMatSharedPtr RTmodprism)
 Modify the prism transformation matrix to align with the tetrahedral modes. More...
 
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom (void)
 Sets up the reference tretrahedral element needed to construct a low energy basis. More...
 
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom (void)
 Sets up the reference prismatic element needed to construct a low energy basis. More...
 
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom (void)
 Sets up the reference hexahedral element needed to construct a low energy basis. More...
 
virtual void v_InitObject ()
 
virtual void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
virtual void v_DoTransformToLowEnergy (Array< OneD, NekDouble > &pInOut, int offset)
 Transform the solution vector vector to low energy. More...
 
virtual void v_DoTransformToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Transform the solution vector vector to low energy. More...
 
virtual void v_DoTransformFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 transform the solution vector from low energy back to the original basis. More...
 
virtual void v_DoMultiplybyInverseTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block inverse transformation matrix. More...
 
virtual void v_DoMultiplybyInverseTransposedTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block tranposed inverse transformation matrix. More...
 
virtual void v_BuildPreconditioner ()
 Construct the low energy preconditioner from $\mathbf{S}_{2}$. More...
 
virtual DNekScalMatSharedPtr v_TransformedSchurCompl (int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
 Set up the transformed block matrix system. More...
 

Additional Inherited Members

Detailed Description

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

Definition at line 53 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

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

Definition at line 64 of file PreconditionerLowEnergy.cpp.

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

Definition at line 75 of file PreconditionerLowEnergy.h.

75 {}

Member Function Documentation

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

Creates an instance of this class.

Definition at line 57 of file PreconditionerLowEnergy.h.

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

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

Create the inverse multiplicity map.

Definition at line 1271 of file PreconditionerLowEnergy.cpp.

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

Referenced by v_InitObject().

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

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

Definition at line 1510 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

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

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

Definition at line 1328 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

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

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

Definition at line 1418 of file PreconditionerLowEnergy.cpp.

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

Referenced by SetUpReferenceElements().

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

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

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

Definition at line 1867 of file PreconditionerLowEnergy.cpp.

Referenced by SetUpReferenceElements().

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

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

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

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

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

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

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

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 118 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.

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

Multiply by the block inverse transformation matrix.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1127 of file PreconditionerLowEnergy.cpp.

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

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

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

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

Apply the low energy preconditioner during the conjugate gradient routine

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 891 of file PreconditionerLowEnergy.cpp.

References Nektar::eWrapper, and m_locToGloMap.

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

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

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

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

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

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 73 of file PreconditionerLowEnergy.cpp.

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

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

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

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 Tue Dec 8 2015 23:45:39 for Nektar++ by   doxygen 1.8.9.1