48 namespace MultiRegions
76 GlobalLinSysStaticCond::GlobalLinSysStaticCond(
78 const boost::weak_ptr<ExpList> &pExpList,
79 const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
81 m_locToGloMap (pLocToGloMap)
112 bool dirForcCalculated = (bool) dirForcing.num_elements();
113 bool atLastLevel = pLocToGloMap->AtLastLevel();
114 int scLevel = pLocToGloMap->GetStaticCondLevel();
116 int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
117 int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
118 int nDirBndDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
119 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
120 int nLocBndDofs = pLocToGloMap->GetNumLocalBndCoeffs();
121 int nIntDofs = pLocToGloMap->GetNumGlobalCoeffs()
126 if(nDirBndDofs && dirForcCalculated)
128 Vmath::Vsub(nGlobDofs,in.get(),1,dirForcing.get(),1,F.get(),1);
156 if( nIntDofs && ((!dirForcCalculated) && (atLastLevel)) )
162 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
163 V_LocBnd = BinvD*F_Int + SchurCompl*V_LocBnd;
165 else if((!dirForcCalculated) && (atLastLevel))
169 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
170 V_LocBnd = SchurCompl*V_LocBnd;
178 pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
180 Subtract(F_HomBnd, F_HomBnd, V_GlobHomBndTmp);
191 int lcLevel = pLocToGloMap->GetLowestStaticCondLevel();
192 if(atLastLevel && scLevel < lcLevel)
198 for (
int i = scLevel; i < lcLevel; ++i)
201 pLocToGloMap->UniversalAssembleBnd(tmp);
203 tmp.get()+nDirBndDofs, 1,
204 V_GlobHomBndTmp.
GetPtr().get(), 1);
205 Subtract( F_HomBnd, F_HomBnd, V_GlobHomBndTmp);
216 nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
223 &pert[nDirBndDofs],1,&out[nDirBndDofs],1);
229 pLocToGloMap->GetNextLevelLocalToGlobalMap());
238 if(nGlobHomBndDofs || nDirBndDofs)
242 if(dirForcCalculated && nDirBndDofs)
244 pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd,
249 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
251 F_Int = F_Int - C*V_LocBnd;
266 const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
270 int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
271 int nDirBndDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
272 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
274 (2*nLocalBnd + nGlobal + nGlobHomBndDofs, 0.0);
276 if (pLocToGloMap->AtLastLevel())
283 pLocToGloMap->GetNextLevelLocalToGlobalMap());
300 const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
303 int n_exp =
m_expList.lock()->GetNumElmts();
306 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
308 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
321 for(n = 0; n < n_exp; ++n)
335 m_schurCompl->SetBlock(n, n, t = loc_schur->GetBlock(0,0));
336 m_BinvD ->SetBlock(n, n, t = loc_schur->GetBlock(0,1));
337 m_C ->SetBlock(n, n, t = loc_schur->GetBlock(1,0));
338 m_invD ->SetBlock(n, n, t = loc_schur->GetBlock(1,1));
361 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
363 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
381 int nPatches = pLocToGloMap->GetNumPatches();
382 int nEntriesA = 0;
int nEntriesB = 0;
383 int nEntriesC = 0;
int nEntriesD = 0;
385 for(i = 0; i < nPatches; i++)
387 nEntriesA += nBndDofsPerPatch[i]*nBndDofsPerPatch[i];
388 nEntriesB += nBndDofsPerPatch[i]*nIntDofsPerPatch[i];
389 nEntriesC += nIntDofsPerPatch[i]*nBndDofsPerPatch[i];
390 nEntriesD += nIntDofsPerPatch[i]*nIntDofsPerPatch[i];
425 for(i = 0; i < nPatches; i++)
428 tmparray = storageA+cntA;
434 tmparray = storageB+cntB;
440 tmparray = storageC+cntC;
446 tmparray = storageD+cntD;
453 cntA += nBndDofsPerPatch[i] * nBndDofsPerPatch[i];
454 cntB += nBndDofsPerPatch[i] * nIntDofsPerPatch[i];
455 cntC += nIntDofsPerPatch[i] * nBndDofsPerPatch[i];
456 cntD += nIntDofsPerPatch[i] * nIntDofsPerPatch[i];
463 int schurComplSubMatnRows;
468 for(n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
470 schurComplSubMat = SchurCompl->GetBlock(n,n);
471 schurComplSubMatnRows = schurComplSubMat->GetRows();
473 patchId = pLocToGloMap->GetPatchMapFromPrevLevel()
474 ->GetPatchId() + cnt;
475 dofId = pLocToGloMap->GetPatchMapFromPrevLevel()
477 isBndDof = pLocToGloMap->GetPatchMapFromPrevLevel()
479 sign = pLocToGloMap->GetPatchMapFromPrevLevel()
483 for(i = 0; i < schurComplSubMatnRows; ++i)
485 int pId = patchId[i];
487 = substructuredMat[0][pId]->GetPtr();
489 = substructuredMat[1][patchId[i]]->GetPtr();
491 = substructuredMat[2][patchId[i]]->GetPtr();
493 = substructuredMat[3][patchId[i]];
494 int subMat0rows = substructuredMat[0][pId]->GetRows();
495 int subMat1rows = substructuredMat[1][pId]->GetRows();
496 int subMat2rows = substructuredMat[2][pId]->GetRows();
500 for(j = 0; j < schurComplSubMatnRows; ++j)
503 "These values should be equal");
507 subMat0[dofId[i]+dofId[j]*subMat0rows] +=
509 (*schurComplSubMat)(i,j);
513 subMat1[dofId[i]+dofId[j]*subMat1rows] +=
515 (*schurComplSubMat)(i,j);
521 for(j = 0; j < schurComplSubMatnRows; ++j)
524 "These values should be equal");
528 subMat2[dofId[i]+dofId[j]*subMat2rows] +=
530 (*schurComplSubMat)(i,j);
536 if (dofId[i] <= dofId[j])
538 (*subMat3)(dofId[i],dofId[j]) +=
540 (*schurComplSubMat)(i,j);
545 (*subMat3)(dofId[i],dofId[j]) +=
547 (*schurComplSubMat)(i,j);
553 cnt += schurComplSubMatnRows;
558 for(i = 0; i < nPatches; i++)
560 if(nIntDofsPerPatch[i])
563 = substructuredMat[0][i]->GetPtr();
565 = substructuredMat[1][i]->GetPtr();
567 = substructuredMat[2][i]->GetPtr();
568 int subMat0rows = substructuredMat[0][i]->GetRows();
569 int subMat1rows = substructuredMat[1][i]->GetRows();
570 int subMat2rows = substructuredMat[2][i]->GetRows();
571 int subMat2cols = substructuredMat[2][i]->GetColumns();
574 substructuredMat[3][i]->Invert();
576 (*substructuredMat[1][i]) = (*substructuredMat[1][i])*
577 (*substructuredMat[3][i]);
582 Blas::Dgemm(
'N',
'N', subMat1rows, subMat2cols,
583 subMat2rows, -1.0, &subMat1[0], subMat1rows,
584 &subMat2[0], subMat2rows, 1.0, &subMat0[0],
593 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
595 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
608 for(i = 0; i < nPatches; i++)
610 for(j = 0; j < 4; j++)
614 blkMatrices[j]->SetBlock(i,i,tmpscalmat);
627 blkMatrices[2], blkMatrices[3], pLocToGloMap);
Array< OneD, NekDouble > m_wsp
Workspace array for matrix multiplication.
#define ASSERTL0(condition, msg)
virtual void v_Solve(const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Solve the linear system for given input and output vectors using a specified local to global map...
virtual GlobalLinSysStaticCondSharedPtr v_Recurse(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const DNekScalBlkMatSharedPtr pSchurCompl, const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, const DNekScalBlkMatSharedPtr pInvD, const boost::shared_ptr< AssemblyMap > &locToGloMap)=0
void DiagonalBlockFullScalMatrixMultiply(NekVector< double > &result, const NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > &lhs, const NekVector< double > &rhs)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
virtual void v_InitObject()
DNekScalBlkMatSharedPtr m_schurCompl
Block Schur complement matrix.
GlobalLinSysStaticCondSharedPtr m_recursiveSchurCompl
Schur complement for Direct Static Condensation.
DNekScalBlkMatSharedPtr m_C
Block matrix.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir=0)
Solve the linear system for given input and output vectors.
virtual void v_BasisTransform(Array< OneD, NekDouble > &pInOut, int offset)
virtual void v_Initialise(const boost::shared_ptr< AssemblyMap > &locToGloMap)
Initialise this object.
void Initialise(const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n)
Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided ...
DNekScalBlkMatSharedPtr m_invD
Block matrix.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
virtual void v_AssembleSchurComplement(boost::shared_ptr< AssemblyMap > pLoctoGloMap)
void Multiply(NekMatrix< ResultDataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const ResultDataType &rhs)
virtual DNekScalMatSharedPtr v_GetBlock(unsigned int n)
Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey...
virtual DNekScalBlkMatSharedPtr v_PreSolve(int scLevel, NekVector< NekDouble > &F_GlobBnd)
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
void SetupTopLevel(const boost::shared_ptr< AssemblyMap > &locToGloMap)
Set up the storage for the Schur complement or the top level of the multi-level Schur complement...
Describe a linear system.
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
virtual int v_GetNumBlocks()
Get the number of blocks in this system.
virtual ~GlobalLinSysStaticCond()
boost::shared_ptr< AssemblyMap > m_locToGloMap
Local to global map.
RhsMatrixType void Subtract(NekMatrix< DataType, StandardMatrixTag > &result, const NekMatrix< LhsDataType, LhsMatrixType > &lhs, const NekMatrix< RhsDataType, RhsMatrixType > &rhs)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
virtual void v_BasisInvTransform(Array< OneD, NekDouble > &pInOut)
DNekScalBlkMatSharedPtr m_BinvD
Block matrix.
Array< OneD, DataType > & GetPtr()
void ConstructNextLevelCondensedSystem(const boost::shared_ptr< AssemblyMap > &locToGloMap)
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.