Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GlobalLinSysDirectStaticCond.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File GlobalLinSys.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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: GlobalLinSys definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40  namespace MultiRegions
41  {
42  /**
43  * @class GlobalLinSysDirect
44  *
45  * Solves a linear system using single- or multi-level static
46  * condensation.
47  */
48 
49  /**
50  * Registers the class with the Factory.
51  */
54  "DirectStaticCond",
56  "Direct static condensation.");
57 
60  "DirectMultiLevelStaticCond",
62  "Direct multi-level static condensation.");
63 
64  /**
65  * For a matrix system of the form @f[
66  * \left[ \begin{array}{cc}
67  * \boldsymbol{A} & \boldsymbol{B}\\
68  * \boldsymbol{C} & \boldsymbol{D}
69  * \end{array} \right]
70  * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2}
71  * \end{array}\right]
72  * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2}
73  * \end{array}\right],
74  * @f]
75  * where @f$\boldsymbol{D}@f$ and
76  * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble
77  * a static condensation system, according to a given local to global
78  * mapping. #m_linSys is constructed by AssembleSchurComplement().
79  * @param mKey Associated matrix key.
80  * @param pLocMatSys LocalMatrixSystem
81  * @param locToGloMap Local to global mapping.
82  */
84  const GlobalLinSysKey &pKey,
85  const boost::weak_ptr<ExpList> &pExpList,
86  const boost::shared_ptr<AssemblyMap>
87  &pLocToGloMap)
88  : GlobalLinSys (pKey, pExpList, pLocToGloMap),
89  GlobalLinSysDirect (pKey, pExpList, pLocToGloMap),
90  GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
91  {
94  "This constructor is only valid when using static "
95  "condensation");
97  == pLocToGloMap->GetGlobalSysSolnType(),
98  "The local to global map is not set up for the requested "
99  "solution type");
100  }
101 
102  /**
103  *
104  */
106  const GlobalLinSysKey &pKey,
107  const boost::weak_ptr<ExpList> &pExpList,
108  const DNekScalBlkMatSharedPtr pSchurCompl,
109  const DNekScalBlkMatSharedPtr pBinvD,
110  const DNekScalBlkMatSharedPtr pC,
111  const DNekScalBlkMatSharedPtr pInvD,
112  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
113  : GlobalLinSys (pKey, pExpList, pLocToGloMap),
114  GlobalLinSysDirect (pKey, pExpList, pLocToGloMap),
115  GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
116  {
117  m_schurCompl = pSchurCompl;
118  m_BinvD = pBinvD;
119  m_C = pC;
120  m_invD = pInvD;
121  }
122 
123 
124  /**
125  *
126  */
128  {
129 
130  }
131 
132 
134  const AssemblyMapSharedPtr &pLocToGloMap)
135  {
136  int nBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
137  int NumDirBCs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
138  unsigned int rows = nBndDofs - NumDirBCs;
139  int bwidth = pLocToGloMap->GetBndSystemBandWidth();
140 
141  MatrixStorage matStorage;
142 
143  switch(m_linSysKey.GetMatrixType())
144  {
145  // case for all symmetric matices
146  case StdRegions::eMass:
150  {
151  if( (2*(bwidth+1)) < rows)
152  {
154  }
155  else
156  {
157  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
158  }
159  }
160  break;
163  default:
164  {
165  // Current inversion techniques do not seem to
166  // allow banded matrices to be used as a linear
167  // system
168  matStorage = eFULL;
169  }
170  break;
171  }
172 
173  return matStorage;
174  }
175 
176  /**
177  * Assemble the schur complement matrix from the block matrices stored
178  * in #m_blkMatrices and the given local to global mapping information.
179  * @param locToGloMap Local to global mapping information.
180  */
182  const AssemblyMapSharedPtr pLocToGloMap)
183  {
184  int i, j, n, cnt, gid1, gid2;
185  NekDouble sign1, sign2, value;
186 
187  int nBndDofs = pLocToGloMap->GetNumGlobalBndCoeffs();
188  int NumDirBCs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
189 
194 
195  unsigned int rows = nBndDofs - NumDirBCs;
196  unsigned int cols = nBndDofs - NumDirBCs;
197 
198  DNekMatSharedPtr Gmat;
199  int bwidth = pLocToGloMap->GetBndSystemBandWidth();
200 
201  MatrixStorage matStorage = DetermineMatrixStorage(pLocToGloMap);
202 
203  switch(matStorage)
204  {
206  {
207  try {
209  ::AllocateSharedPtr(rows, cols, 0.0,
210  matStorage,
211  bwidth, bwidth);
212  }
213  catch (...) {
215  "Insufficient memory for GlobalLinSys.");
216  }
217  break;
218  }
219 
221  case eFULL:
222  {
224  ::AllocateSharedPtr(rows, cols, 0.0, matStorage);
225  break;
226  }
227 
228  default:
229  {
231  "Unknown matrix storage type of type not set up");
232  }
233  }
234 
235  // fill global matrix
236  DNekScalMatSharedPtr loc_mat;
237  int loc_lda;
238  for(n = cnt = 0; n < SchurCompl->GetNumberOfBlockRows(); ++n)
239  {
240  loc_mat = SchurCompl->GetBlock(n,n);
241  loc_lda = loc_mat->GetRows();
242 
243  // Set up Matrix;
244  for(i = 0; i < loc_lda; ++i)
245  {
246  gid1 = pLocToGloMap->GetLocalToGlobalBndMap (cnt + i)
247  - NumDirBCs;
248  sign1 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + i);
249 
250  if(gid1 >= 0)
251  {
252  for(j = 0; j < loc_lda; ++j)
253  {
254  gid2 = pLocToGloMap->GetLocalToGlobalBndMap(cnt+j)
255  - NumDirBCs;
256  sign2 = pLocToGloMap->GetLocalToGlobalBndSign(cnt+j);
257 
258  if(gid2 >= 0)
259  {
260  // As the global matrix should be symmetric,
261  // only add the value for the upper triangular
262  // part in order to avoid entries to be entered
263  // twice
264  if((matStorage == eFULL)||(gid2 >= gid1))
265  {
266  value = Gmat->GetValue(gid1,gid2)
267  + sign1*sign2*(*loc_mat)(i,j);
268  Gmat->SetValue(gid1,gid2,value);
269  }
270  }
271  }
272  }
273  }
274  cnt += loc_lda;
275  }
276 
277  if(rows)
278  {
281  }
282  }
283 
285  const GlobalLinSysKey &mkey,
286  const boost::weak_ptr<ExpList> &pExpList,
287  const DNekScalBlkMatSharedPtr pSchurCompl,
288  const DNekScalBlkMatSharedPtr pBinvD,
289  const DNekScalBlkMatSharedPtr pC,
290  const DNekScalBlkMatSharedPtr pInvD,
291  const boost::shared_ptr<AssemblyMap> &l2gMap)
292  {
294  GlobalLinSysDirectStaticCond>::AllocateSharedPtr(
295  mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap);
296  sys->Initialise(l2gMap);
297  return sys;
298  }
299  }
300 }
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:185
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
MatrixStorage DetermineMatrixStorage(const boost::shared_ptr< AssemblyMap > &locToGloMap)
Matrix Storage type for known matrices.
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
DNekScalBlkMatSharedPtr m_schurCompl
Block Schur complement matrix.
DNekLinSysSharedPtr m_linSys
Basic linear system object.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
DNekScalBlkMatSharedPtr m_C
Block matrix.
boost::shared_ptr< GlobalLinSysDirectStaticCond > GlobalLinSysDirectStaticCondSharedPtr
DNekScalBlkMatSharedPtr m_invD
Block matrix.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
GlobalLinSysDirectStaticCond(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &locToGloMap)
Constructor for full direct matrix solve.
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
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 > &l2gMap)
double NekDouble
Describe a linear system.
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
static GlobalLinSysSharedPtr create(const GlobalLinSysKey &pLinSysKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:127
A global linear system.
Definition: GlobalLinSys.h:74
virtual void v_AssembleSchurComplement(boost::shared_ptr< AssemblyMap > pLocToGloMap)
GlobalLinSysFactory & GetGlobalLinSysFactory()
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
DNekScalBlkMatSharedPtr m_BinvD
Block matrix.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
boost::shared_ptr< GlobalLinSysStaticCond > GlobalLinSysStaticCondSharedPtr