Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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  nel = expList->GetOffset_Elmt_Id(n);
166  exp = expList->GetExp(nel);
167 
168  // Grab reference to local Schur complement matrix.
169  DNekScalMatSharedPtr schurMat =
170  m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
171 
172  // Process vertices to extract relevant portion of the Schur
173  // complement matrix.
174  for (i = 0; i < exp->GetNverts(); ++i)
175  {
176  meshVertId = exp->GetGeom()->GetVid(i);
177  int locId = exp->GetVertexMap(i);
178 
179  // Get the global ID of this vertex.
180  gId = m_locToGloMap->GetLocalToGlobalMap(
181  cnt + locId) - nDirBnd;
182 
183  // Ignore all Dirichlet vertices.
184  if (gId < 0)
185  {
186  continue;
187  }
188 
189  gidDofs[0][gId] = 1;
190 
191  // Extract vertex value from Schur complement matrix.
192  NekDouble vertVal = (*schurMat)(locId,locId);
193 
194  // See if we have processed this vertex from another
195  // element.
196  gIt = idMats[0].find(gId);
197 
198  if (gIt == idMats[0].end())
199  {
200  // If not then put our 'matrix' inside idMats.
201  idMats[0][gId] = vector<NekDouble>(1, vertVal);
202  }
203  else
204  {
205  // Otherwise combine with the value that is already
206  // there (i.e. do assembly on this degree of freedom).
207  gIt->second[0] += vertVal;
208  }
209 
210  // Now check to see if the vertex is periodic. If it is,
211  // then we change meshVertId to be the minimum of all the
212  // other periodic vertices, so that we don't end up
213  // duplicating the matrix in our final block matrix.
214  pIt = periodicVerts.find(meshVertId);
215  if (pIt != periodicVerts.end())
216  {
217  for (j = 0; j < pIt->second.size(); ++j)
218  {
219  meshVertId = min(meshVertId, pIt->second[j].id);
220  }
221  }
222 
223  // Finally record the other information we need into the
224  // other matrices.
225  gidMeshIds[0][gId] = meshVertId;
226  maxVertIds[0] = max(maxVertIds[0], meshVertId);
227  maxVertIds[1] = 1;
228  }
229 
230  // Process edges. This logic is mostly the same as the previous
231  // block.
232  for (i = 0; i < exp->GetNedges(); ++i)
233  {
234  meshEdgeId = exp->GetGeom()->GetEid(i);
235 
236  Array<OneD, unsigned int> bmap, bmap2;
238  StdRegions::Orientation edgeOrient =
239  exp->GetGeom()->GetEorient(i);
240 
241  // Check if this edge is periodic. We may need to flip
242  // orientation if it is.
243  pIt = periodicEdges.find(meshEdgeId);
244  if (pIt != periodicEdges.end())
245  {
246  pair<int, StdRegions::Orientation> idOrient =
248  meshEdgeId, edgeOrient, pIt->second);
249  meshEdgeId = idOrient.first;
250  edgeOrient = idOrient.second;
251  }
252 
253  // Grab edge interior map, and the edge inverse boundary
254  // map, so that we can extract this edge from the Schur
255  // complement matrix.
256  exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
257  bmap2 = exp->GetEdgeInverseBoundaryMap(i);
258 
259  // Allocate temporary storage for the extracted edge matrix.
260  const int nEdgeCoeffs = bmap.num_elements();
261  vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
262 
263  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
264 
265  for (j = 0; j < nEdgeCoeffs; ++j)
266  {
267  // We record the minimum ID from the edge for our
268  // maps. This follows the logic that the assembly map
269  // ordering will always give us a contiguous ordering of
270  // global degrees of freedom for edge interior
271  // coefficients.
272  gId = min(gId,
273  m_locToGloMap->GetLocalToGlobalMap(
274  cnt + bmap[j])
275  - nDirBnd);
276 
277  // Ignore Dirichlet edges.
278  if (gId < 0)
279  {
280  continue;
281  }
282 
283  const NekDouble sign1 = sign[j];
284 
285  // Extract this edge, along with sign array for assembly
286  // later.
287  for (k = 0; k < nEdgeCoeffs; ++k)
288  {
289  tmpStore[k+j*nEdgeCoeffs] =
290  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
291  }
292  }
293 
294  if (gId < 0)
295  {
296  continue;
297  }
298 
299  gidDofs[1][gId] = nEdgeCoeffs;
300 
301  // Assemble this edge matrix with another one, if it exists.
302  gIt = idMats[1].find(gId);
303  if (gIt == idMats[1].end())
304  {
305  idMats[1][gId] = tmpStore;
306  }
307  else
308  {
309  ASSERTL1(tmpStore.size() == gIt->second.size(),
310  "Number of modes mismatch");
311  Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
312  &tmpStore[0], 1, &gIt->second[0], 1);
313  }
314 
315  gidMeshIds[1][gId] = meshEdgeId;
316  maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
317  maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
318  }
319 
320  // Process faces. This logic is mostly the same as the previous
321  // block.
322  for (i = 0; i < exp->GetNfaces(); ++i)
323  {
324  meshFaceId = exp->GetGeom()->GetFid(i);
325 
326  Array<OneD, unsigned int> bmap, bmap2;
328  StdRegions::Orientation faceOrient =
329  exp->GetGeom()->GetForient(i);
330 
331  // Check if this face is periodic. We may need to flip
332  // orientation if it is.
333  pIt = periodicFaces.find(meshFaceId);
334  if (pIt != periodicFaces.end())
335  {
336  meshFaceId = min(meshFaceId, pIt->second[0].id);
337  faceOrient = DeterminePeriodicFaceOrient(
338  faceOrient, pIt->second[0].orient);
339  }
340 
341  exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
342  bmap2 = exp->GetFaceInverseBoundaryMap(i);
343 
344  // Allocate temporary storage for the extracted face matrix.
345  const int nFaceCoeffs = bmap.num_elements();
346  vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
347 
348  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
349 
350  for (j = 0; j < nFaceCoeffs; ++j)
351  {
352  gId = min(gId,
353  m_locToGloMap->GetLocalToGlobalMap(
354  cnt + bmap[j])
355  - nDirBnd);
356 
357  // Ignore Dirichlet faces.
358  if (gId < 0)
359  {
360  continue;
361  }
362 
363  const NekDouble sign1 = sign[j];
364 
365  // Extract this face, along with sign array for assembly
366  // later.
367  for (k = 0; k < nFaceCoeffs; ++k)
368  {
369  tmpStore[k+j*nFaceCoeffs] =
370  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
371  }
372  }
373 
374  if (gId < 0)
375  {
376  continue;
377  }
378 
379  gidDofs[2][gId] = nFaceCoeffs;
380 
381  // Assemble this face matrix with another one, if it exists.
382  gIt = idMats[2].find(gId);
383  if (gIt == idMats[2].end())
384  {
385  idMats[2][gId] = tmpStore;
386  }
387  else
388  {
389  ASSERTL1(tmpStore.size() == gIt->second.size(),
390  "Number of modes mismatch");
391  Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
392  &tmpStore[0], 1, &gIt->second[0], 1);
393  }
394 
395  gidMeshIds[2][gId] = meshFaceId;
396  maxVertIds[4] = max(maxVertIds[4], meshFaceId);
397  maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
398  }
399 
400  cnt += exp->GetNcoeffs();
401  }
402 
403  // Perform a reduction to find maximum vertex, edge and face
404  // geometry IDs.
405  m_comm = expList->GetSession()->GetComm()->GetRowComm();
406  m_comm->AllReduce(maxVertIds, LibUtilities::ReduceMax);
407 
408  // Concatenate all matrices into contiguous storage and figure out
409  // universal ID numbering.
410  vector<NekDouble> storageBuf;
411  vector<long> globalToUniversal;
412 
413  for (i = 0, cnt = 1; i < 3; ++i)
414  {
415  const int maxDofs = maxVertIds[2*i+1];
416 
417  // Note that iterating over the map uses the property that keys
418  // are accessed in order of ascending order, putting everything
419  // in the right order for the global system.
420  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
421  {
422  // Copy matrix into storage.
423  storageBuf.insert(storageBuf.end(),
424  gIt->second.begin(), gIt->second.end());
425 
426  // Get mesh ID from global ID number.
427  ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
428  "Unable to find global ID " +
429  boost::lexical_cast<string>(gIt->first) +
430  " inside map");
431  meshVertId = gidMeshIds[i][gIt->first];
432 
433  for (j = 0; j < gIt->second.size(); ++j)
434  {
435  globalToUniversal.push_back(
436  cnt + meshVertId*maxDofs*maxDofs + j);
437  }
438 
439  // Free up the temporary storage.
440  gIt->second.clear();
441  }
442 
443  cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
444  }
445 
446  ASSERTL1(storageBuf.size() == globalToUniversal.size(),
447  "Storage buffer and global to universal map size does "
448  "not match");
449 
450  Array<OneD, NekDouble> storageData(
451  storageBuf.size(), &storageBuf[0]);
452  Array<OneD, long> globalToUniversalMap(
453  globalToUniversal.size(), &globalToUniversal[0]);
454 
455  // Use GS to assemble data between processors.
456  Gs::gs_data *tmpGs = Gs::Init(globalToUniversalMap, m_comm);
457  Gs::Gather(storageData, Gs::gs_add, tmpGs);
458 
459  // Figure out what storage we need in the block matrix.
461  1 + idMats[1].size() + idMats[2].size());
462 
463  // Vertex block is a diagonal matrix.
464  n_blks[0] = idMats[0].size();
465 
466  // Now extract number of rows in each edge and face block from the
467  // gidDofs map.
468  cnt = 1;
469  for (i = 1; i < 3; ++i)
470  {
471  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
472  {
473  ASSERTL1(gidDofs[i].count(gIt->first) > 0,
474  "Unable to find number of degrees of freedom for "
475  "global ID " + boost::lexical_cast<string>(
476  gIt->first));
477 
478  n_blks[cnt++] = gidDofs[i][gIt->first];
479  }
480  }
481 
482  // Allocate storage for the block matrix.
484  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
485 
486  // We deal with the vertex matrix separately since all vertices can
487  // be condensed into a single, block-diagonal matrix.
489  ::AllocateSharedPtr(n_blks[0], n_blks[0], 0.0, eDIAGONAL);
490 
491  // Fill the vertex matrix with the inverse of each vertex value.
492  cnt = 0;
493  for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
494  {
495  (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
496  }
497 
498  // Put the vertex matrix in the block matrix.
499  m_blkMat->SetBlock(0,0,vertMat);
500 
501  // Finally, grab the matrices from the block storage, invert them
502  // and place them in the correct position inside the block matrix.
503  int cnt2 = 1;
504  for (i = 1; i < 3; ++i)
505  {
506  for (gIt = idMats[i].begin(); gIt != idMats[i].end();
507  ++gIt, ++cnt2)
508  {
509  int nDofs = gidDofs[i][gIt->first];
510 
512  ::AllocateSharedPtr(nDofs, nDofs);
513 
514  for (j = 0; j < nDofs; ++j)
515  {
516  for (k = 0; k < nDofs; ++k)
517  {
518  (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
519  }
520  }
521 
522  cnt += nDofs*nDofs;
523 
524  tmp->Invert();
525  m_blkMat->SetBlock(cnt2, cnt2, tmp);
526  }
527  }
528  }
529 
530  /**
531  * @brief Construct a block preconditioner for the hybridized
532  * discontinuous Galerkin system.
533  *
534  * This system is built in a similar fashion to its continuous variant
535  * found in PreconditionerBlock::BlockPreconditionerCG. In this setting
536  * however, the matrix is constructed as:
537  *
538  * \f[ M^{-1} = \mathrm{Diag}[ (\mathbf{S_{1}})_{f}^{-1} ] \f]
539  *
540  * where each matrix is the Schur complement system restricted to a
541  * single face of the trace system.
542  */
544  {
545  boost::shared_ptr<MultiRegions::ExpList>
546  expList = ((m_linsys.lock())->GetLocMat()).lock();
547  boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
549  DNekScalBlkMatSharedPtr loc_mat;
550  DNekScalMatSharedPtr bnd_mat;
551 
552  AssemblyMapDGSharedPtr asmMap = boost::dynamic_pointer_cast<
554 
555  int i, j, k, n, cnt, cnt2;
556 
557  // Figure out number of Dirichlet trace elements
558  int nTrace = expList->GetTrace()->GetExpSize();
559  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
560 
561  for (cnt = n = 0; n < nTrace; ++n)
562  {
563  if (cnt >= nDir)
564  {
565  break;
566  }
567 
568  cnt += trace->GetExp(n)->GetNcoeffs();
569  }
570 
571  nDir = n;
572 
573  // Allocate storage for block matrix. Need number of unique faces in
574  // trace space.
575  int maxTraceSize = -1;
576  map<int, int> arrayOffsets;
577 
578  for (cnt = 0, n = nDir; n < nTrace; ++n)
579  {
580  int nCoeffs = trace->GetExp(n)->GetNcoeffs();
581  int nCoeffs2 = nCoeffs * nCoeffs;
582 
583  arrayOffsets[n] = cnt;
584  cnt += nCoeffs2;
585 
586  if (maxTraceSize < nCoeffs2)
587  {
588  maxTraceSize = nCoeffs2;
589  }
590  }
591 
592  // Find maximum trace size.
593  m_comm = expList->GetSession()->GetComm()->GetRowComm();
594  m_comm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
595 
596  // Zero matrix storage.
597  Array<OneD, NekDouble> tmpStore(cnt, 0.0);
598 
599  // Assemble block matrices for each trace element.
600  for (cnt = n = 0; n < expList->GetExpSize(); ++n)
601  {
602  int elmt = expList->GetOffset_Elmt_Id(n);
603  locExpansion = expList->GetExp(elmt);
604 
606  asmMap->GetElmtToTrace()[elmt];
607 
608  // Block matrix (lambda)
609  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
610  bnd_mat = loc_mat->GetBlock(0,0);
611 
612  int nFacets = locExpansion->GetNumBases() == 2 ?
613  locExpansion->GetNedges() : locExpansion->GetNfaces();
614 
615  for (cnt2 = i = 0; i < nFacets; ++i)
616  {
617  int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
618  int elmtId = elmtToTraceMap[i]->GetElmtId ();
619 
620  // Ignore Dirichlet edges.
621  if (elmtId < nDir)
622  {
623  cnt += nCoeffs;
624  cnt2 += nCoeffs;
625  continue;
626  }
627 
628  NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
629 
630  for (j = 0; j < nCoeffs; ++j)
631  {
632  NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
633  cnt + j);
634  for (k = 0; k < nCoeffs; ++k)
635  {
636  NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
637  cnt + k);
638  off[j*nCoeffs + k] +=
639  (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
640  }
641  }
642 
643  cnt += nCoeffs;
644  cnt2 += nCoeffs;
645  }
646  }
647 
648  // Set up IDs for universal numbering.
649  Array<OneD, long> uniIds(tmpStore.num_elements());
650  for (cnt = 0, n = nDir; n < nTrace; ++n)
651  {
652  LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
653  int nCoeffs = traceExp->GetNcoeffs();
654  int geomId = traceExp->GetGeom()->GetGlobalID();
655 
656  for (i = 0; i < nCoeffs*nCoeffs; ++i)
657  {
658  uniIds[cnt++] = geomId * maxTraceSize + i + 1;
659  }
660  }
661 
662  // Assemble matrices across partitions.
663  Gs::gs_data *gsh = Gs::Init(uniIds, m_comm);
664  Gs::Gather(tmpStore, Gs::gs_add, gsh);
665 
666  // Set up diagonal block matrix
667  Array<OneD, unsigned int> n_blks(nTrace - nDir);
668  for (n = 0; n < nTrace - nDir; ++n)
669  {
670  n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
671  }
672 
674  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
675 
676  for (n = 0; n < nTrace - nDir; ++n)
677  {
678  int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
680  ::AllocateSharedPtr(nCoeffs, nCoeffs);
681  NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
682 
683  for (i = 0; i < nCoeffs; ++i)
684  {
685  for (j = 0; j < nCoeffs; ++j)
686  {
687  (*tmp)(i,j) = off[i*nCoeffs + j];
688  }
689  }
690 
691  tmp->Invert();
692  m_blkMat->SetBlock(n, n, tmp);
693  }
694  }
695 
696  /**
697  * @brief Apply preconditioner to \p pInput and store the result in \p
698  * pOutput.
699  */
701  const Array<OneD, NekDouble>& pInput,
702  Array<OneD, NekDouble>& pOutput)
703  {
704  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
705  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
706  int nNonDir = nGlobal-nDir;
707  DNekBlkMat &M = (*m_blkMat);
708  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
709  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
710  z = M * r;
711  }
712  }
713 }
boost::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:49
LibUtilities::CommSharedPtr m_comm
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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:224
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:22
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
PreconFactory & GetPreconFactory()
STL namespace.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
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:191
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:285
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215