46 namespace MultiRegions
76 const boost::weak_ptr<ExpList> &pExpList,
77 const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
79 m_locToGloMap (pLocToGloMap)
110 bool dirForcCalculated = (bool) dirForcing.num_elements();
111 bool atLastLevel = pLocToGloMap->AtLastLevel();
112 int scLevel = pLocToGloMap->GetStaticCondLevel();
114 int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
115 int nGlobBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
116 int nDirBndDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
117 int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
118 int nLocBndDofs = pLocToGloMap->GetNumLocalBndCoeffs();
119 int nIntDofs = pLocToGloMap->GetNumGlobalCoeffs()
124 if(nDirBndDofs && dirForcCalculated)
126 Vmath::Vsub(nGlobDofs,in.get(),1,dirForcing.get(),1,F.get(),1);
154 if( nIntDofs && ((!dirForcCalculated) && (atLastLevel)) )
160 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
161 V_LocBnd = BinvD*F_Int + SchurCompl*V_LocBnd;
163 else if((!dirForcCalculated) && (atLastLevel))
167 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
168 V_LocBnd = SchurCompl*V_LocBnd;
173 V_LocBnd = BinvD*F_Int;
176 pLocToGloMap->AssembleBnd(V_LocBnd,V_GlobHomBndTmp,
178 F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
189 int lcLevel = pLocToGloMap->GetLowestStaticCondLevel();
190 if(atLastLevel && scLevel < lcLevel)
196 for (
int i = scLevel; i < lcLevel; ++i)
199 pLocToGloMap->UniversalAssembleBnd(tmp);
201 tmp.get()+nDirBndDofs, 1,
202 V_GlobHomBndTmp.GetPtr().get(), 1);
203 F_HomBnd = F_HomBnd - V_GlobHomBndTmp;
215 nGlobBndDofs, F, pert, pLocToGloMap, nDirBndDofs);
222 &pert[nDirBndDofs],1,&out[nDirBndDofs],1);
228 pLocToGloMap->GetNextLevelLocalToGlobalMap());
237 if(nGlobHomBndDofs || nDirBndDofs)
241 if(dirForcCalculated && nDirBndDofs)
243 pLocToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd,
248 pLocToGloMap->GlobalToLocalBnd(V_GlobBnd,V_LocBnd);
250 F_Int = F_Int - C*V_LocBnd;
266 const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
272 if (pLocToGloMap->AtLastLevel())
279 pLocToGloMap->GetNextLevelLocalToGlobalMap());
296 const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
299 int n_exp =
m_expList.lock()->GetNumElmts();
302 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
304 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
317 for(n = 0; n < n_exp; ++n)
333 m_schurCompl->SetBlock(n, n, t = loc_schur->GetBlock(0,0));
334 m_BinvD ->SetBlock(n, n, t = loc_schur->GetBlock(0,1));
335 m_C ->SetBlock(n, n, t = loc_schur->GetBlock(1,0));
336 m_invD ->SetBlock(n, n, t = loc_schur->GetBlock(1,1));
359 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
361 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
379 int nPatches = pLocToGloMap->GetNumPatches();
380 int nEntriesA = 0;
int nEntriesB = 0;
381 int nEntriesC = 0;
int nEntriesD = 0;
383 for(i = 0; i < nPatches; i++)
385 nEntriesA += nBndDofsPerPatch[i]*nBndDofsPerPatch[i];
386 nEntriesB += nBndDofsPerPatch[i]*nIntDofsPerPatch[i];
387 nEntriesC += nIntDofsPerPatch[i]*nBndDofsPerPatch[i];
388 nEntriesD += nIntDofsPerPatch[i]*nIntDofsPerPatch[i];
414 for(i = 0; i < nPatches; i++)
417 tmparray = storageA+cntA;
423 tmparray = storageB+cntB;
429 tmparray = storageC+cntC;
435 tmparray = storageD+cntD;
441 cntA += nBndDofsPerPatch[i] * nBndDofsPerPatch[i];
442 cntB += nBndDofsPerPatch[i] * nIntDofsPerPatch[i];
443 cntC += nIntDofsPerPatch[i] * nBndDofsPerPatch[i];
444 cntD += nIntDofsPerPatch[i] * nIntDofsPerPatch[i];
451 int schurComplSubMatnRows;
457 for(n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
459 schurComplSubMat = SchurCompl->GetBlock(n,n);
460 schurComplSubMatnRows = schurComplSubMat->GetRows();
462 scale = SchurCompl->GetBlock(n,n)->Scale();
464 = SchurCompl->GetBlock(n,n)->GetOwnedMatrix()->GetPtr();
466 patchId = pLocToGloMap->GetPatchMapFromPrevLevel()
467 ->GetPatchId() + cnt;
468 dofId = pLocToGloMap->GetPatchMapFromPrevLevel()
470 isBndDof = pLocToGloMap->GetPatchMapFromPrevLevel()
472 sign = pLocToGloMap->GetPatchMapFromPrevLevel()
476 for(i = 0; i < schurComplSubMatnRows; ++i)
478 int pId = patchId[i];
480 = substructuredMat[0][pId]->GetPtr();
482 = substructuredMat[1][patchId[i]]->GetPtr();
484 = substructuredMat[2][patchId[i]]->GetPtr();
486 = substructuredMat[3][patchId[i]]->GetPtr();
487 int subMat0rows = substructuredMat[0][pId]->GetRows();
488 int subMat1rows = substructuredMat[1][pId]->GetRows();
489 int subMat2rows = substructuredMat[2][pId]->GetRows();
490 int subMat3rows = substructuredMat[3][pId]->GetRows();
494 for(j = 0; j < schurComplSubMatnRows; ++j)
497 "These values should be equal");
501 subMat0[dofId[i]+dofId[j]*subMat0rows] +=
504 i+j*schurComplSubMatnRows]);
508 subMat1[dofId[i]+dofId[j]*subMat1rows] +=
511 i+j*schurComplSubMatnRows]);
517 for(j = 0; j < schurComplSubMatnRows; ++j)
520 "These values should be equal");
524 subMat2[dofId[i]+dofId[j]*subMat2rows] +=
527 i+j*schurComplSubMatnRows]);
531 subMat3[dofId[i]+dofId[j]*subMat3rows] +=
534 i+j*schurComplSubMatnRows]);
539 cnt += schurComplSubMatnRows;
544 for(i = 0; i < nPatches; i++)
546 if(nIntDofsPerPatch[i])
549 = substructuredMat[0][i]->GetPtr();
551 = substructuredMat[1][i]->GetPtr();
553 = substructuredMat[2][i]->GetPtr();
554 int subMat0rows = substructuredMat[0][i]->GetRows();
555 int subMat1rows = substructuredMat[1][i]->GetRows();
556 int subMat2rows = substructuredMat[2][i]->GetRows();
557 int subMat2cols = substructuredMat[2][i]->GetColumns();
560 substructuredMat[3][i]->Invert();
562 (*substructuredMat[1][i]) = (*substructuredMat[1][i])*
563 (*substructuredMat[3][i]);
568 Blas::Dgemm(
'N',
'N', subMat1rows, subMat2cols,
569 subMat2rows, -1.0, &subMat1[0], subMat1rows,
570 &subMat2[0], subMat2rows, 1.0, &subMat0[0],
579 = pLocToGloMap->GetNumLocalBndCoeffsPerPatch();
581 = pLocToGloMap->GetNumLocalIntCoeffsPerPatch();
594 for(i = 0; i < nPatches; i++)
596 for(j = 0; j < 4; j++)
600 blkMatrices[j]->SetBlock(i,i,tmpscalmat);
613 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.
GlobalLinSysStaticCond(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &locToGloMap)
Constructor for full direct matrix solve.
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.