Nektar++
GlobalLinSysDirectStaticCond.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: GlobalLinSysDirectStaticCond.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: GlobalLinSysDirectStaticCond definition
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
38{
39/**
40 * @class GlobalLinSysDirect
41 *
42 * Solves a linear system using single- or multi-level static
43 * condensation.
44 */
45
46/**
47 * Registers the class with the Factory.
48 */
51 "DirectStaticCond", GlobalLinSysDirectStaticCond::create,
52 "Direct static condensation.");
53
56 "DirectMultiLevelStaticCond", GlobalLinSysDirectStaticCond::create,
57 "Direct multi-level static condensation.");
58
59/**
60 * For a matrix system of the form @f[
61 * \left[ \begin{array}{cc}
62 * \boldsymbol{A} & \boldsymbol{B}\\
63 * \boldsymbol{C} & \boldsymbol{D}
64 * \end{array} \right]
65 * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2}
66 * \end{array}\right]
67 * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2}
68 * \end{array}\right],
69 * @f]
70 * where @f$\boldsymbol{D}@f$ and
71 * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble
72 * a static condensation system, according to a given local to global
73 * mapping. #m_linSys is constructed by AssembleSchurComplement().
74 * @param mKey Associated matrix key.
75 * @param pLocMatSys LocalMatrixSystem
76 * @param locToGloMap Local to global mapping.
77 */
79 const GlobalLinSysKey &pKey, const std::weak_ptr<ExpList> &pExpList,
80 const std::shared_ptr<AssemblyMap> &pLocToGloMap)
81 : GlobalLinSys(pKey, pExpList, pLocToGloMap),
82 GlobalLinSysDirect(pKey, pExpList, pLocToGloMap),
83 GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
84{
87 "This constructor is only valid when using static "
88 "condensation");
90 pLocToGloMap->GetGlobalSysSolnType(),
91 "The local to global map is not set up for the requested "
92 "solution type");
93}
94
95/**
96 *
97 */
99 const GlobalLinSysKey &pKey, const std::weak_ptr<ExpList> &pExpList,
100 const DNekScalBlkMatSharedPtr pSchurCompl,
102 const DNekScalBlkMatSharedPtr pInvD,
103 const std::shared_ptr<AssemblyMap> &pLocToGloMap)
104 : GlobalLinSys(pKey, pExpList, pLocToGloMap),
105 GlobalLinSysDirect(pKey, pExpList, pLocToGloMap),
106 GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
107{
108 m_schurCompl = pSchurCompl;
109 m_BinvD = pBinvD;
110 m_C = pC;
111 m_invD = pInvD;
112}
113
114/**
115 *
116 */
118{
119}
120
122 const AssemblyMapSharedPtr &pLocToGloMap)
123{
124 int nBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
125 int NumDirBCs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
126 unsigned int rows = nBndDofs - NumDirBCs;
127 int bwidth = pLocToGloMap->GetBndSystemBandWidth();
128
129 MatrixStorage matStorage;
130
131 switch (m_linSysKey.GetMatrixType())
132 {
133 // case for all symmetric matices
138 {
139 if ((2 * (bwidth + 1)) < rows)
140 {
142 }
143 else
144 {
145 matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
146 }
147 }
148 break;
152 default:
153 {
154 // Current inversion techniques do not seem to
155 // allow banded matrices to be used as a linear
156 // system
157 matStorage = eFULL;
158 }
159 break;
160 }
161
162 return matStorage;
163}
164
165/**
166 * Assemble the schur complement matrix from the block matrices stored
167 * in #m_blkMatrices and the given local to global mapping information.
168 * @param locToGloMap Local to global mapping information.
169 */
171 const AssemblyMapSharedPtr pLocToGloMap)
172{
173 int i, j, n, cnt, gid1, gid2;
174 NekDouble sign1, sign2, value;
175
176 int nBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
177 int NumDirBCs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
178
183
184 unsigned int rows = nBndDofs - NumDirBCs;
185 unsigned int cols = nBndDofs - NumDirBCs;
186
187 DNekMatSharedPtr Gmat;
188 int bwidth = pLocToGloMap->GetBndSystemBandWidth();
189
190 MatrixStorage matStorage = DetermineMatrixStorage(pLocToGloMap);
191
192 switch (matStorage)
193 {
195 {
196 try
197 {
199 rows, cols, 0.0, matStorage, bwidth, bwidth);
200 }
201 catch (...)
202 {
204 "Insufficient memory for GlobalLinSys.");
205 }
206 break;
207 }
208
210 case eFULL:
211 {
212 Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols, 0.0,
213 matStorage);
214 break;
215 }
216
217 default:
218 {
220 "Unknown matrix storage type of type not set up");
221 }
222 }
223
224 // fill global matrix
225 DNekScalMatSharedPtr loc_mat;
226 int loc_lda;
227 for (n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
228 {
229 loc_mat = SchurCompl->GetBlock(n, n);
230 loc_lda = loc_mat->GetRows();
231
232 // Set up Matrix;
233 for (i = 0; i < loc_lda; ++i)
234 {
235 gid1 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + i) - NumDirBCs;
236 sign1 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + i);
237
238 if (gid1 >= 0)
239 {
240 for (j = 0; j < loc_lda; ++j)
241 {
242 gid2 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + j) -
243 NumDirBCs;
244 sign2 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + j);
245
246 if (gid2 >= 0)
247 {
248 // As the global matrix should be symmetric,
249 // only add the value for the upper triangular
250 // part in order to avoid entries to be entered
251 // twice
252 if ((matStorage == eFULL) || (gid2 >= gid1))
253 {
254 value = Gmat->GetValue(gid1, gid2) +
255 sign1 * sign2 * (*loc_mat)(i, j);
256 Gmat->SetValue(gid1, gid2, value);
257 }
258 }
259 }
260 }
261 }
262 cnt += loc_lda;
263 }
264
265 if (rows)
266 {
269 }
270}
271
273 const GlobalLinSysKey &mkey, const std::weak_ptr<ExpList> &pExpList,
274 const DNekScalBlkMatSharedPtr pSchurCompl,
276 const DNekScalBlkMatSharedPtr pInvD,
277 const std::shared_ptr<AssemblyMap> &l2gMap)
278{
281 mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap);
282 sys->Initialise(l2gMap);
283 return sys;
284}
285
286/// Solve the linear system for given input and output vectors.
288 const int pNumRows, const Array<OneD, const NekDouble> &pInput,
289 Array<OneD, NekDouble> &pOutput, const AssemblyMapSharedPtr &pLocToGloMap,
290 const int pNumDir)
291{
292 Array<OneD, NekDouble> tmp(pNumRows);
293 Array<OneD, NekDouble> global(pNumRows, 0.0);
294
295 pLocToGloMap->AssembleBnd(pInput, tmp);
296
297 const int nHomDofs = pNumRows - pNumDir;
298 DNekVec Vin(nHomDofs, tmp + pNumDir);
299
300 Array<OneD, NekDouble> tmp1 = global + pNumDir;
301 DNekVec Vout(nHomDofs, tmp1, eWrapper);
302
303 m_linSys->Solve(Vin, Vout);
304
305 pLocToGloMap->GlobalToLocalBnd(global, pOutput);
306}
307
308} // namespace Nektar::MultiRegions
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
DNekLinSysSharedPtr m_linSys
Basic linear system object.
static GlobalLinSysSharedPtr create(const GlobalLinSysKey &pLinSysKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
GlobalLinSysDirectStaticCond(const GlobalLinSysKey &mkey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &locToGloMap)
Constructor for full direct matrix solve.
MatrixStorage DetermineMatrixStorage(const std::shared_ptr< AssemblyMap > &locToGloMap)
Matrix Storage type for known matrices.
GlobalLinSysStaticCondSharedPtr v_Recurse(const GlobalLinSysKey &mkey, const std::weak_ptr< ExpList > &pExpList, const DNekScalBlkMatSharedPtr pSchurCompl, const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, const DNekScalBlkMatSharedPtr pInvD, const std::shared_ptr< AssemblyMap > &l2gMap) override
void v_SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir) override
Solve the linear system for given input and output vectors.
void v_AssembleSchurComplement(std::shared_ptr< AssemblyMap > pLocToGloMap) override
A global linear system.
Definition: GlobalLinSys.h:70
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:120
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
DNekScalBlkMatSharedPtr m_schurCompl
Block Schur complement matrix.
DNekScalBlkMatSharedPtr m_BinvD
Block matrix.
DNekScalBlkMatSharedPtr m_C
Block matrix.
DNekScalBlkMatSharedPtr m_invD
Block matrix.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
std::shared_ptr< GlobalLinSysDirectStaticCond > GlobalLinSysDirectStaticCondSharedPtr
std::shared_ptr< GlobalLinSysStaticCond > GlobalLinSysStaticCondSharedPtr
GlobalLinSysFactory & GetGlobalLinSysFactory()
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:50
std::vector< double > w(NPUPPER)
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...