Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PreconditionerBlock.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PreconditionerBlock.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: Block Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
42 #include <LocalRegions/MatrixKey.h>
43 #include <LocalRegions/SegExp.h>
44 #include <math.h>
45 
46 using namespace std;
47 
48 namespace Nektar
49 {
50  namespace MultiRegions
51  {
52  /**
53  * Registers the class with the Factory.
54  */
55  string PreconditionerBlock::className
57  "Block",
58  PreconditionerBlock::create,
59  "Block Diagonal Preconditioning");
60  /**
61  * @class Block Preconditioner
62  *
63  * This class implements Block preconditioning for the conjugate
64  * gradient matrix solver.
65  */
66 
67  PreconditionerBlock::PreconditionerBlock(
68  const boost::shared_ptr<GlobalLinSys> &plinsys,
69  const AssemblyMapSharedPtr &pLocToGloMap)
70  : Preconditioner(plinsys, pLocToGloMap),
71  m_linsys(plinsys),
72  m_preconType(pLocToGloMap->GetPreconType()),
73  m_locToGloMap(pLocToGloMap)
74  {
75  }
76 
78  {
79  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
81  solvertype == MultiRegions::ePETScStaticCond,
82  "Solver type not valid");
83  }
84 
86  {
87  GlobalLinSysKey key = m_linsys.lock()->GetKey();
88 
89  // Different setup for HDG and CG.
91  {
93  }
94  else
95  {
97  }
98  }
99 
100  /**
101  * \brief Construct a block preconditioner from \f$\mathbf{S}_{1}\f$ for
102  * the continuous Galerkin system.
103  *
104  * The preconditioner is defined as:
105  *
106  * \f[\mathbf{M}^{-1}=\left[\begin{array}{ccc}
107  * \mathrm{Diag}[(\mathbf{S_{1}})_{vv}] & & \\
108  * & (\mathbf{S}_{1})_{eb} & \\
109  * & & (\mathbf{S}_{1})_{fb} \end{array}\right] \f]
110  *
111  * where \f$\mathbf{S}_{1}\f$ is the local Schur complement matrix for
112  * each element and the subscript denotes the portion of the Schur
113  * complement associated with the vertex, edge and face blocks
114  * respectively.
115  */
117  {
118  ExpListSharedPtr expList = m_linsys.lock()->GetLocMat().lock();
120  DNekScalBlkMatSharedPtr loc_mat;
121  DNekScalMatSharedPtr bnd_mat;
122 
123  int nel, i, j, k, n, cnt, gId;
124  int meshVertId, meshEdgeId, meshFaceId;
125 
126  const int nExp = expList->GetExpSize();
127  const int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
128 
129  // Grab periodic geometry information.
131  PeriodicMap periodicVerts, periodicEdges, periodicFaces;
132  expList->GetPeriodicEntities(
133  periodicVerts, periodicEdges, periodicFaces);
134 
135  // The vectors below are of size 3 to have separate storage for
136  // vertices, edges and faces.
137 
138  // Maps from geometry ID to the matrix representing the extracted
139  // portion of S_1. For example idMats[2] folds the S_1 face blocks.
140  vector<map<int, vector<NekDouble> > > idMats(3);
141 
142  // Maps from the global ID, as obtained from AssemblyMapCG's
143  // localToGlobalMap, to the geometry ID.
144  vector<map<int, int> > gidMeshIds(3);
145 
146  // Maps from the global ID to the number of degrees of freedom for
147  // this geometry object.
148  vector<map<int, int> > gidDofs(3);
149 
150  // Array containing maximum information needed for the universal
151  // numbering later. For i = 0,1,2 for each geometry dimension:
152  //
153  // maxVertIds[2*i] = maximum geometry ID at dimension i
154  // maxVertIds[2*i+1] = maximum number of degrees of freedom for all
155  // elements of dimension i.
156  Array<OneD, int> maxVertIds(6, -1);
157 
158  // Iterator for idMats members.
159  map<int, vector<NekDouble> >::iterator gIt;
160 
161  // Figure out mapping from each elemental contribution to offset in
162  // (vert,edge,face) triples.
163  for (cnt = n = 0; n < nExp; ++n)
164  {
165  exp = expList->GetExp(n);
166 
167  // Grab reference to local Schur complement matrix.
168  DNekScalMatSharedPtr schurMat =
169  m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
170 
171  // Process vertices to extract relevant portion of the Schur
172  // complement matrix.
173  for (i = 0; i < exp->GetNverts(); ++i)
174  {
175  meshVertId = exp->GetGeom()->GetVid(i);
176  int locId = exp->GetVertexMap(i);
177 
178  // Get the global ID of this vertex.
179  gId = m_locToGloMap->GetLocalToGlobalMap(
180  cnt + locId) - nDirBnd;
181 
182  // Ignore all Dirichlet vertices.
183  if (gId < 0)
184  {
185  continue;
186  }
187 
188  gidDofs[0][gId] = 1;
189 
190  // Extract vertex value from Schur complement matrix.
191  NekDouble vertVal = (*schurMat)(locId,locId);
192 
193  // See if we have processed this vertex from another
194  // element.
195  gIt = idMats[0].find(gId);
196 
197  if (gIt == idMats[0].end())
198  {
199  // If not then put our 'matrix' inside idMats.
200  idMats[0][gId] = vector<NekDouble>(1, vertVal);
201  }
202  else
203  {
204  // Otherwise combine with the value that is already
205  // there (i.e. do assembly on this degree of freedom).
206  gIt->second[0] += vertVal;
207  }
208 
209  // Now check to see if the vertex is periodic. If it is,
210  // then we change meshVertId to be the minimum of all the
211  // other periodic vertices, so that we don't end up
212  // duplicating the matrix in our final block matrix.
213  pIt = periodicVerts.find(meshVertId);
214  if (pIt != periodicVerts.end())
215  {
216  for (j = 0; j < pIt->second.size(); ++j)
217  {
218  meshVertId = min(meshVertId, pIt->second[j].id);
219  }
220  }
221 
222  // Finally record the other information we need into the
223  // other matrices.
224  gidMeshIds[0][gId] = meshVertId;
225  maxVertIds[0] = max(maxVertIds[0], meshVertId);
226  maxVertIds[1] = 1;
227  }
228 
229  // Process edges. This logic is mostly the same as the previous
230  // block.
231  for (i = 0; i < exp->GetNedges(); ++i)
232  {
233  meshEdgeId = exp->GetGeom()->GetEid(i);
234 
235  Array<OneD, unsigned int> bmap, bmap2;
237  StdRegions::Orientation edgeOrient =
238  exp->GetGeom()->GetEorient(i);
239 
240  // Check if this edge is periodic. We may need to flip
241  // orientation if it is.
242  pIt = periodicEdges.find(meshEdgeId);
243  if (pIt != periodicEdges.end())
244  {
245  pair<int, StdRegions::Orientation> idOrient =
247  meshEdgeId, edgeOrient, pIt->second);
248  meshEdgeId = idOrient.first;
249  edgeOrient = idOrient.second;
250  }
251 
252  // Grab edge interior map, and the edge inverse boundary
253  // map, so that we can extract this edge from the Schur
254  // complement matrix.
255  exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
256  bmap2 = exp->GetEdgeInverseBoundaryMap(i);
257 
258  // Allocate temporary storage for the extracted edge matrix.
259  const int nEdgeCoeffs = bmap.num_elements();
260  vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
261 
262  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
263 
264  for (j = 0; j < nEdgeCoeffs; ++j)
265  {
266  // We record the minimum ID from the edge for our
267  // maps. This follows the logic that the assembly map
268  // ordering will always give us a contiguous ordering of
269  // global degrees of freedom for edge interior
270  // coefficients.
271  gId = min(gId,
272  m_locToGloMap->GetLocalToGlobalMap(
273  cnt + bmap[j])
274  - nDirBnd);
275 
276  // Ignore Dirichlet edges.
277  if (gId < 0)
278  {
279  continue;
280  }
281 
282  const NekDouble sign1 = sign[j];
283 
284  // Extract this edge, along with sign array for assembly
285  // later.
286  for (k = 0; k < nEdgeCoeffs; ++k)
287  {
288  tmpStore[k+j*nEdgeCoeffs] =
289  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
290  }
291  }
292 
293  if (gId < 0)
294  {
295  continue;
296  }
297 
298  gidDofs[1][gId] = nEdgeCoeffs;
299 
300  // Assemble this edge matrix with another one, if it exists.
301  gIt = idMats[1].find(gId);
302  if (gIt == idMats[1].end())
303  {
304  idMats[1][gId] = tmpStore;
305  }
306  else
307  {
308  ASSERTL1(tmpStore.size() == gIt->second.size(),
309  "Number of modes mismatch");
310  Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
311  &tmpStore[0], 1, &gIt->second[0], 1);
312  }
313 
314  gidMeshIds[1][gId] = meshEdgeId;
315  maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
316  maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
317  }
318 
319  // Process faces. This logic is mostly the same as the previous
320  // block.
321  for (i = 0; i < exp->GetNfaces(); ++i)
322  {
323  meshFaceId = exp->GetGeom()->GetFid(i);
324 
325  Array<OneD, unsigned int> bmap, bmap2;
327  StdRegions::Orientation faceOrient =
328  exp->GetGeom()->GetForient(i);
329 
330  // Check if this face is periodic. We may need to flip
331  // orientation if it is.
332  pIt = periodicFaces.find(meshFaceId);
333  if (pIt != periodicFaces.end())
334  {
335  meshFaceId = min(meshFaceId, pIt->second[0].id);
336  faceOrient = DeterminePeriodicFaceOrient(
337  faceOrient, pIt->second[0].orient);
338  }
339 
340  exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
341  bmap2 = exp->GetFaceInverseBoundaryMap(i);
342 
343  // Allocate temporary storage for the extracted face matrix.
344  const int nFaceCoeffs = bmap.num_elements();
345  vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
346 
347  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
348 
349  for (j = 0; j < nFaceCoeffs; ++j)
350  {
351  gId = min(gId,
352  m_locToGloMap->GetLocalToGlobalMap(
353  cnt + bmap[j])
354  - nDirBnd);
355 
356  // Ignore Dirichlet faces.
357  if (gId < 0)
358  {
359  continue;
360  }
361 
362  const NekDouble sign1 = sign[j];
363 
364  // Extract this face, along with sign array for assembly
365  // later.
366  for (k = 0; k < nFaceCoeffs; ++k)
367  {
368  tmpStore[k+j*nFaceCoeffs] =
369  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
370  }
371  }
372 
373  if (gId < 0)
374  {
375  continue;
376  }
377 
378  gidDofs[2][gId] = nFaceCoeffs;
379 
380  // Assemble this face matrix with another one, if it exists.
381  gIt = idMats[2].find(gId);
382  if (gIt == idMats[2].end())
383  {
384  idMats[2][gId] = tmpStore;
385  }
386  else
387  {
388  ASSERTL1(tmpStore.size() == gIt->second.size(),
389  "Number of modes mismatch");
390  Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
391  &tmpStore[0], 1, &gIt->second[0], 1);
392  }
393 
394  gidMeshIds[2][gId] = meshFaceId;
395  maxVertIds[4] = max(maxVertIds[4], meshFaceId);
396  maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
397  }
398 
399  cnt += exp->GetNcoeffs();
400  }
401 
402  // Perform a reduction to find maximum vertex, edge and face
403  // geometry IDs.
404  m_comm = expList->GetSession()->GetComm()->GetRowComm();
405  m_comm->AllReduce(maxVertIds, LibUtilities::ReduceMax);
406 
407  // Concatenate all matrices into contiguous storage and figure out
408  // universal ID numbering.
409  vector<NekDouble> storageBuf;
410  vector<long> globalToUniversal;
411 
412  for (i = 0, cnt = 1; i < 3; ++i)
413  {
414  const int maxDofs = maxVertIds[2*i+1];
415 
416  // Note that iterating over the map uses the property that keys
417  // are accessed in order of ascending order, putting everything
418  // in the right order for the global system.
419  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
420  {
421  // Copy matrix into storage.
422  storageBuf.insert(storageBuf.end(),
423  gIt->second.begin(), gIt->second.end());
424 
425  // Get mesh ID from global ID number.
426  ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
427  "Unable to find global ID " +
428  boost::lexical_cast<string>(gIt->first) +
429  " inside map");
430  meshVertId = gidMeshIds[i][gIt->first];
431 
432  for (j = 0; j < gIt->second.size(); ++j)
433  {
434  globalToUniversal.push_back(
435  cnt + meshVertId*maxDofs*maxDofs + j);
436  }
437 
438  // Free up the temporary storage.
439  gIt->second.clear();
440  }
441 
442  cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
443  }
444 
445  ASSERTL1(storageBuf.size() == globalToUniversal.size(),
446  "Storage buffer and global to universal map size does "
447  "not match");
448 
449  Array<OneD, NekDouble> storageData(
450  storageBuf.size(), &storageBuf[0]);
451  Array<OneD, long> globalToUniversalMap(
452  globalToUniversal.size(), &globalToUniversal[0]);
453 
454  // Use GS to assemble data between processors.
455  Gs::gs_data *tmpGs = Gs::Init(
456  globalToUniversalMap, m_comm,
457  expList->GetSession()->DefinesCmdLineArgument("verbose"));
458  Gs::Gather(storageData, Gs::gs_add, tmpGs);
459 
460  // Figure out what storage we need in the block matrix.
462  1 + idMats[1].size() + idMats[2].size());
463 
464  // Vertex block is a diagonal matrix.
465  n_blks[0] = idMats[0].size();
466 
467  // Now extract number of rows in each edge and face block from the
468  // gidDofs map.
469  cnt = 1;
470  for (i = 1; i < 3; ++i)
471  {
472  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
473  {
474  ASSERTL1(gidDofs[i].count(gIt->first) > 0,
475  "Unable to find number of degrees of freedom for "
476  "global ID " + boost::lexical_cast<string>(
477  gIt->first));
478 
479  n_blks[cnt++] = gidDofs[i][gIt->first];
480  }
481  }
482 
483  // Allocate storage for the block matrix.
485  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
486 
487  // We deal with the vertex matrix separately since all vertices can
488  // be condensed into a single, block-diagonal matrix.
490  ::AllocateSharedPtr(n_blks[0], n_blks[0], 0.0, eDIAGONAL);
491 
492  // Fill the vertex matrix with the inverse of each vertex value.
493  cnt = 0;
494  for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
495  {
496  (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
497  }
498 
499  // Put the vertex matrix in the block matrix.
500  m_blkMat->SetBlock(0,0,vertMat);
501 
502  // Finally, grab the matrices from the block storage, invert them
503  // and place them in the correct position inside the block matrix.
504  int cnt2 = 1;
505  for (i = 1; i < 3; ++i)
506  {
507  for (gIt = idMats[i].begin(); gIt != idMats[i].end();
508  ++gIt, ++cnt2)
509  {
510  int nDofs = gidDofs[i][gIt->first];
511 
513  ::AllocateSharedPtr(nDofs, nDofs);
514 
515  for (j = 0; j < nDofs; ++j)
516  {
517  for (k = 0; k < nDofs; ++k)
518  {
519  (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
520  }
521  }
522 
523  cnt += nDofs*nDofs;
524 
525  tmp->Invert();
526  m_blkMat->SetBlock(cnt2, cnt2, tmp);
527  }
528  }
529  }
530 
531  /**
532  * @brief Construct a block preconditioner for the hybridized
533  * discontinuous Galerkin system.
534  *
535  * This system is built in a similar fashion to its continuous variant
536  * found in PreconditionerBlock::BlockPreconditionerCG. In this setting
537  * however, the matrix is constructed as:
538  *
539  * \f[ M^{-1} = \mathrm{Diag}[ (\mathbf{S_{1}})_{f}^{-1} ] \f]
540  *
541  * where each matrix is the Schur complement system restricted to a
542  * single face of the trace system.
543  */
545  {
546  boost::shared_ptr<MultiRegions::ExpList>
547  expList = ((m_linsys.lock())->GetLocMat()).lock();
548  boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
550  DNekScalBlkMatSharedPtr loc_mat;
551  DNekScalMatSharedPtr bnd_mat;
552 
553  AssemblyMapDGSharedPtr asmMap = boost::dynamic_pointer_cast<
555 
556  int i, j, k, n, cnt, cnt2;
557 
558  // Figure out number of Dirichlet trace elements
559  int nTrace = expList->GetTrace()->GetExpSize();
560  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
561 
562  for (cnt = n = 0; n < nTrace; ++n)
563  {
564  if (cnt >= nDir)
565  {
566  break;
567  }
568 
569  cnt += trace->GetExp(n)->GetNcoeffs();
570  }
571 
572  nDir = n;
573 
574  // Allocate storage for block matrix. Need number of unique faces in
575  // trace space.
576  int maxTraceSize = -1;
577  map<int, int> arrayOffsets;
578 
579  for (cnt = 0, n = nDir; n < nTrace; ++n)
580  {
581  int nCoeffs = trace->GetExp(n)->GetNcoeffs();
582  int nCoeffs2 = nCoeffs * nCoeffs;
583 
584  arrayOffsets[n] = cnt;
585  cnt += nCoeffs2;
586 
587  if (maxTraceSize < nCoeffs2)
588  {
589  maxTraceSize = nCoeffs2;
590  }
591  }
592 
593  // Find maximum trace size.
594  m_comm = expList->GetSession()->GetComm()->GetRowComm();
595  m_comm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
596 
597  // Zero matrix storage.
598  Array<OneD, NekDouble> tmpStore(cnt, 0.0);
599 
600  // Assemble block matrices for each trace element.
601  for (cnt = n = 0; n < expList->GetExpSize(); ++n)
602  {
603  int elmt = n;
604  locExpansion = expList->GetExp(elmt);
605 
607  asmMap->GetElmtToTrace()[elmt];
608 
609  // Block matrix (lambda)
610  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
611  bnd_mat = loc_mat->GetBlock(0,0);
612 
613  int nFacets = locExpansion->GetNumBases() == 2 ?
614  locExpansion->GetNedges() : locExpansion->GetNfaces();
615 
616  for (cnt2 = i = 0; i < nFacets; ++i)
617  {
618  int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
619  int elmtId = elmtToTraceMap[i]->GetElmtId ();
620 
621  // Ignore Dirichlet edges.
622  if (elmtId < nDir)
623  {
624  cnt += nCoeffs;
625  cnt2 += nCoeffs;
626  continue;
627  }
628 
629  NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
630 
631  for (j = 0; j < nCoeffs; ++j)
632  {
633  NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
634  cnt + j);
635  for (k = 0; k < nCoeffs; ++k)
636  {
637  NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
638  cnt + k);
639  off[j*nCoeffs + k] +=
640  (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
641  }
642  }
643 
644  cnt += nCoeffs;
645  cnt2 += nCoeffs;
646  }
647  }
648 
649  // Set up IDs for universal numbering.
650  Array<OneD, long> uniIds(tmpStore.num_elements());
651  for (cnt = 0, n = nDir; n < nTrace; ++n)
652  {
653  LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
654  int nCoeffs = traceExp->GetNcoeffs();
655  int geomId = traceExp->GetGeom()->GetGlobalID();
656 
657  for (i = 0; i < nCoeffs*nCoeffs; ++i)
658  {
659  uniIds[cnt++] = geomId * maxTraceSize + i + 1;
660  }
661  }
662 
663  // Assemble matrices across partitions.
664  Gs::gs_data *gsh = Gs::Init(
665  uniIds, m_comm,
666  expList->GetSession()->DefinesCmdLineArgument("verbose"));
667  Gs::Gather(tmpStore, Gs::gs_add, gsh);
668 
669  // Set up diagonal block matrix
670  Array<OneD, unsigned int> n_blks(nTrace - nDir);
671  for (n = 0; n < nTrace - nDir; ++n)
672  {
673  n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
674  }
675 
677  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
678 
679  for (n = 0; n < nTrace - nDir; ++n)
680  {
681  int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
683  ::AllocateSharedPtr(nCoeffs, nCoeffs);
684  NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
685 
686  for (i = 0; i < nCoeffs; ++i)
687  {
688  for (j = 0; j < nCoeffs; ++j)
689  {
690  (*tmp)(i,j) = off[i*nCoeffs + j];
691  }
692  }
693 
694  tmp->Invert();
695  m_blkMat->SetBlock(n, n, tmp);
696  }
697  }
698 
699  /**
700  * @brief Apply preconditioner to \p pInput and store the result in \p
701  * pOutput.
702  */
704  const Array<OneD, NekDouble>& pInput,
705  Array<OneD, NekDouble>& pOutput)
706  {
707  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
708  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
709  int nNonDir = nGlobal-nDir;
710  DNekBlkMat &M = (*m_blkMat);
711  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
712  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
713  z = M * r;
714  }
715  }
716 }
boost::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:49
LibUtilities::CommSharedPtr m_comm
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Apply preconditioner to pInput and store the result in pOutput.
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:239
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
PreconFactory & GetPreconFactory()
STL namespace.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Describe a linear system.
boost::shared_ptr< AssemblyMap > m_locToGloMap
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
const boost::weak_ptr< GlobalLinSys > m_linsys
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
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.
Definition: Vmath.cpp:299
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215