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

#include <PreconditionerLowEnergy.h>

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

Public Member Functions

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

Static Public Member Functions

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

Static Public Attributes

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

Protected Attributes

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

Private Member Functions

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

Additional Inherited Members

Detailed Description

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

Definition at line 53 of file PreconditionerLowEnergy.h.

Constructor & Destructor Documentation

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

Definition at line 66 of file PreconditionerLowEnergy.cpp.

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

Definition at line 75 of file PreconditionerLowEnergy.h.

75 {}

Member Function Documentation

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

Creates an instance of this class.

Definition at line 57 of file PreconditionerLowEnergy.h.

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

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

Create the inverse multiplicity map.

Definition at line 1274 of file PreconditionerLowEnergy.cpp.

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

Referenced by v_InitObject().

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

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

Referenced by SetUpReferenceElements().

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

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

Referenced by SetUpReferenceElements().

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

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

Referenced by SetUpReferenceElements().

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

Referenced by SetUpReferenceElements().

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

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

Definition at line 915 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().

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

Sets up the reference elements needed by the preconditioner.

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

Definition at line 1602 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().

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

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

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

Multiply by the block inverse transformation matrix.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1130 of file PreconditionerLowEnergy.cpp.

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

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

Multiply by the block tranposed inverse transformation matrix.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1175 of file PreconditionerLowEnergy.cpp.

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

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

References Nektar::eWrapper, and m_locToGloMap.

897  {
898  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
899  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
900  int nNonDir = nGlobal-nDir;
901  DNekBlkMat &M = (*m_BlkMat);
902 
903  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
904  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
905 
906  z = M * r;
907  }
Nektar::DNekBlkMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:60
void Nektar::MultiRegions::PreconditionerLowEnergy::v_DoTransformFromLowEnergy ( Array< OneD, NekDouble > &  pInOut)
privatevirtual

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1088 of file PreconditionerLowEnergy.cpp.

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

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

Transform the solution vector vector to low energy.

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 996 of file PreconditionerLowEnergy.cpp.

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

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

Transform the solution vector vector to low energy.

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 1037 of file PreconditionerLowEnergy.cpp.

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

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 75 of file PreconditionerLowEnergy.cpp.

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

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

Set up the transformed block matrix system.

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

Reimplemented from Nektar::MultiRegions::Preconditioner.

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

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

Member Data Documentation

string Nektar::MultiRegions::PreconditionerLowEnergy::className
static
Initial value:
= GetPreconFactory().RegisterCreatorFunction(
"LowEnergyBlock",
PreconditionerLowEnergy::create,
"LowEnergy Preconditioning")

Name of class.

Registers the class with the Factory.

Definition at line 68 of file PreconditionerLowEnergy.h.

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

Definition at line 82 of file PreconditionerLowEnergy.h.

Referenced by v_BuildPreconditioner().

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

Definition at line 83 of file PreconditionerLowEnergy.h.

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

Definition at line 87 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and v_DoMultiplybyInverseTransformationMatrix().

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

Definition at line 88 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and v_DoMultiplybyInverseTransposedTransformationMatrix().

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

Definition at line 79 of file PreconditionerLowEnergy.h.

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

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

Definition at line 80 of file PreconditionerLowEnergy.h.

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

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

Definition at line 105 of file PreconditionerLowEnergy.h.

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

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

Definition at line 107 of file PreconditionerLowEnergy.h.

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

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

Definition at line 106 of file PreconditionerLowEnergy.h.

Referenced by CreateMultiplicityMap(), and v_DoMultiplybyInverseTransposedTransformationMatrix().

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

Definition at line 85 of file PreconditionerLowEnergy.h.

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

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

Definition at line 95 of file PreconditionerLowEnergy.h.

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

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

Definition at line 97 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

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

Definition at line 102 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

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

Definition at line 92 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

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

Definition at line 100 of file PreconditionerLowEnergy.h.

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

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

Definition at line 86 of file PreconditionerLowEnergy.h.

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

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

Definition at line 90 of file PreconditionerLowEnergy.h.

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

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

Definition at line 96 of file PreconditionerLowEnergy.h.

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

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

Definition at line 98 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

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

Definition at line 103 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

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

Definition at line 93 of file PreconditionerLowEnergy.h.

Referenced by SetupBlockTransformationMatrix(), and SetUpReferenceElements().

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

Definition at line 101 of file PreconditionerLowEnergy.h.

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

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

Definition at line 91 of file PreconditionerLowEnergy.h.

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

Generated on Wed May 4 2016 17:54:34 for Nektar++ by   doxygen 1.8.8