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;
175 V_LocBnd = BinvD*F_Int;
178 pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
180 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 F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
217 nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
224 &pert[nDirBndDofs],1,&out[nDirBndDofs],1);
230 pLocToGloMap->GetNextLevelLocalToGlobalMap());
239 if(nGlobHomBndDofs || nDirBndDofs)
243 if(dirForcCalculated && nDirBndDofs)
245 pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd,
250 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
252 F_Int = F_Int - C*V_LocBnd;
268 const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
274 if (pLocToGloMap->AtLastLevel())
281 pLocToGloMap->GetNextLevelLocalToGlobalMap());
298 const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
301 int n_exp =
m_expList.lock()->GetNumElmts();
304 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
306 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
319 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];
416 for(i = 0; i < nPatches; i++)
419 tmparray = storageA+cntA;
425 tmparray = storageB+cntB;
431 tmparray = storageC+cntC;
437 tmparray = storageD+cntD;
443 cntA += nBndDofsPerPatch[i] * nBndDofsPerPatch[i];
444 cntB += nBndDofsPerPatch[i] * nIntDofsPerPatch[i];
445 cntC += nIntDofsPerPatch[i] * nBndDofsPerPatch[i];
446 cntD += nIntDofsPerPatch[i] * nIntDofsPerPatch[i];
453 int schurComplSubMatnRows;
459 for(n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
461 schurComplSubMat = SchurCompl->GetBlock(n,n);
462 schurComplSubMatnRows = schurComplSubMat->GetRows();
464 scale = SchurCompl->GetBlock(n,n)->Scale();
466 = SchurCompl->GetBlock(n,n)->GetOwnedMatrix()->GetPtr();
468 patchId = pLocToGloMap->GetPatchMapFromPrevLevel()
469 ->GetPatchId() + cnt;
470 dofId = pLocToGloMap->GetPatchMapFromPrevLevel()
472 isBndDof = pLocToGloMap->GetPatchMapFromPrevLevel()
474 sign = pLocToGloMap->GetPatchMapFromPrevLevel()
478 for(i = 0; i < schurComplSubMatnRows; ++i)
480 int pId = patchId[i];
482 = substructuredMat[0][pId]->GetPtr();
484 = substructuredMat[1][patchId[i]]->GetPtr();
486 = substructuredMat[2][patchId[i]]->GetPtr();
488 = substructuredMat[3][patchId[i]]->GetPtr();
489 int subMat0rows = substructuredMat[0][pId]->GetRows();
490 int subMat1rows = substructuredMat[1][pId]->GetRows();
491 int subMat2rows = substructuredMat[2][pId]->GetRows();
492 int subMat3rows = substructuredMat[3][pId]->GetRows();
496 for(j = 0; j < schurComplSubMatnRows; ++j)
499 "These values should be equal");
503 subMat0[dofId[i]+dofId[j]*subMat0rows] +=
506 i+j*schurComplSubMatnRows]);
510 subMat1[dofId[i]+dofId[j]*subMat1rows] +=
513 i+j*schurComplSubMatnRows]);
519 for(j = 0; j < schurComplSubMatnRows; ++j)
522 "These values should be equal");
526 subMat2[dofId[i]+dofId[j]*subMat2rows] +=
529 i+j*schurComplSubMatnRows]);
533 subMat3[dofId[i]+dofId[j]*subMat3rows] +=
536 i+j*schurComplSubMatnRows]);
541 cnt += schurComplSubMatnRows;
546 for(i = 0; i < nPatches; i++)
548 if(nIntDofsPerPatch[i])
551 = substructuredMat[0][i]->GetPtr();
553 = substructuredMat[1][i]->GetPtr();
555 = substructuredMat[2][i]->GetPtr();
556 int subMat0rows = substructuredMat[0][i]->GetRows();
557 int subMat1rows = substructuredMat[1][i]->GetRows();
558 int subMat2rows = substructuredMat[2][i]->GetRows();
559 int subMat2cols = substructuredMat[2][i]->GetColumns();
562 substructuredMat[3][i]->Invert();
564 (*substructuredMat[1][i]) = (*substructuredMat[1][i])*
565 (*substructuredMat[3][i]);
570 Blas::Dgemm(
'N',
'N', subMat1rows, subMat2cols,
571 subMat2rows, -1.0, &subMat1[0], subMat1rows,
572 &subMat2[0], subMat2rows, 1.0, &subMat0[0],
581 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
583 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
596 for(i = 0; i < nPatches; i++)
598 for(j = 0; j < 4; j++)
602 blkMatrices[j]->SetBlock(i,i,tmpscalmat);
615 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
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< DNekScalMat > DNekScalMatSharedPtr
virtual void v_AssembleSchurComplement(boost::shared_ptr< AssemblyMap > pLoctoGloMap)
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.
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.