Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PreconditionerLowEnergy.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Preconditioner.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: Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
41 #include <LocalRegions/MatrixKey.h>
42 #include <math.h>
43 
44 using namespace std;
45 
46 namespace Nektar
47 {
48  namespace MultiRegions
49  {
50  /**
51  * Registers the class with the Factory.
52  */
53  string PreconditionerLowEnergy::className
55  "LowEnergyBlock",
56  PreconditionerLowEnergy::create,
57  "LowEnergy Preconditioning");
58 
59  /**
60  * @class PreconditionerLowEnergy
61  *
62  * This class implements low energy preconditioning for the conjugate
63  * gradient matrix solver.
64  */
65 
66  PreconditionerLowEnergy::PreconditionerLowEnergy(
67  const boost::shared_ptr<GlobalLinSys> &plinsys,
68  const AssemblyMapSharedPtr &pLocToGloMap)
69  : Preconditioner(plinsys, pLocToGloMap),
70  m_linsys(plinsys),
71  m_locToGloMap(pLocToGloMap)
72  {
73  }
74 
76  {
77  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
78  ASSERTL0(solvertype == eIterativeStaticCond ||
79  solvertype == ePETScStaticCond, "Solver type not valid");
80 
81  boost::shared_ptr<MultiRegions::ExpList>
82  expList=((m_linsys.lock())->GetLocMat()).lock();
83 
85 
86  locExpansion = expList->GetExp(0);
87 
88  int nDim = locExpansion->GetShapeDimension();
89 
90  ASSERTL0(nDim==3,
91  "Preconditioner type only valid in 3D");
92 
93  //Sets up reference element and builds transformation matrix
95 
96  //Set up block transformation matrix
98 
99  //Sets up multiplicity map for transformation from global to local
101  }
102 
103 
104  /**
105  * \brief Construct the low energy preconditioner from
106  * \f$\mathbf{S}_{2}\f$
107  *
108  * \f[\mathbf{M}^{-1}=\left[\begin{array}{ccc}
109  * Diag[(\mathbf{S_{2}})_{vv}] & & \\ & (\mathbf{S}_{2})_{eb} & \\ & &
110  * (\mathbf{S}_{2})_{fb} \end{array}\right] \f]
111  *
112  * where \f$\mathbf{R}\f$ is the transformation matrix and
113  * \f$\mathbf{S}_{2}\f$ the Schur complement of the modified basis,
114  * given by
115  *
116  * \f[\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\f]
117  *
118  * where \f$\mathbf{S}_{1}\f$ is the local schur complement matrix for
119  * each element.
120  */
122  {
123  boost::shared_ptr<MultiRegions::ExpList>
124  expList=((m_linsys.lock())->GetLocMat()).lock();
126  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
127  StdRegions::VarCoeffMap vVarCoeffMap;
128  int i, j, k;
129  int nVerts, nEdges,nFaces;
130  int eid, fid, n, cnt, nedgemodes, nfacemodes;
131  NekDouble zero = 0.0;
132 
133  int vMap1, vMap2, sign1, sign2;
134  int m, v, eMap1, eMap2, fMap1, fMap2;
135  int offset, globalrow, globalcol, nCoeffs;
136 
137  // Periodic information
138  PeriodicMap periodicVerts;
139  PeriodicMap periodicEdges;
140  PeriodicMap periodicFaces;
141  expList->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
142 
143  //matrix storage
144  MatrixStorage storage = eFULL;
145  MatrixStorage vertstorage = eDIAGONAL;
146  MatrixStorage blkmatStorage = eDIAGONAL;
147 
148  //local element static condensed matrices
149  DNekScalBlkMatSharedPtr loc_mat;
150  DNekScalMatSharedPtr bnd_mat;
151 
152  DNekMatSharedPtr pRS;
153  DNekMatSharedPtr pRSRT;
154 
155  //Transformation matrices
156  DNekMat R;
157  DNekMat RT;
158  DNekMat RS;
159  DNekMat RSRT;
160 
161  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
162  int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
163 
164  //Vertex, edge and face preconditioner matrices
166  AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
167 
168  Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
169  Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
170 
171  //maps for different element types
172  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
173  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
174 
175  //Transformation matrix
176  transmatrixmap[LibUtilities::eTetrahedron]= m_Rtet;
177  transmatrixmap[LibUtilities::ePrism] = m_Rprism;
178  transmatrixmap[LibUtilities::eHexahedron] = m_Rhex;
179 
180  //Transposed transformation matrix
181  transposedtransmatrixmap[LibUtilities::eTetrahedron]= m_RTtet;
182  transposedtransmatrixmap[LibUtilities::ePrism] = m_RTprism;
183  transposedtransmatrixmap[LibUtilities::eHexahedron] = m_RThex;
184 
185  int n_exp = expList->GetNumElmts();
186  int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
187  int nNonDirFaceIDs=m_locToGloMap->GetNumNonDirFaces();
188 
189  //set the number of blocks in the matrix
190  int numBlks = 1+nNonDirEdgeIDs+nNonDirFaceIDs;
191  Array<OneD,unsigned int> n_blks(numBlks);
192  for(i = 0; i < numBlks; ++i)
193  {
194  n_blks[i] = 0;
195  }
196  n_blks[0]=nNonDirVerts;
197 
198  set<int> edgeDirMap;
199  set<int> faceDirMap;
200  map<int,int> uniqueEdgeMap;
201  map<int,int> uniqueFaceMap;
202 
203  //this should be of size total number of local edges
204  Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
205  Array<OneD, int> facemodeoffset(nNonDirFaceIDs,0);
206 
207  Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
208  Array<OneD, int> faceglobaloffset(nNonDirFaceIDs,0);
209 
210  const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
212  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>& bndConditions = expList->GetBndConditions();
213 
214  int meshVertId;
215  int meshEdgeId;
216  int meshFaceId;
217 
218  const Array<OneD, const int> &extradiredges
219  = m_locToGloMap->GetExtraDirEdges();
220  for(i=0; i<extradiredges.num_elements(); ++i)
221  {
222  meshEdgeId=extradiredges[i];
223  edgeDirMap.insert(meshEdgeId);
224  }
225 
226  //Determine which boundary edges and faces have dirichlet values
227  for(i = 0; i < bndCondExp.num_elements(); i++)
228  {
229  cnt = 0;
230  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
231  {
232  bndCondFaceExp = boost::dynamic_pointer_cast<
233  StdRegions::StdExpansion2D>(bndCondExp[i]->GetExp(j));
234  if (bndConditions[i]->GetBoundaryConditionType() ==
236  {
237  for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
238  {
239  meshEdgeId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(k);
240  if(edgeDirMap.count(meshEdgeId) == 0)
241  {
242  edgeDirMap.insert(meshEdgeId);
243  }
244  }
245  meshFaceId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetFid();
246  faceDirMap.insert(meshFaceId);
247  }
248  }
249  }
250 
251  int dof=0;
252  int maxFaceDof=0;
253  int maxEdgeDof=0;
254  int nlocalNonDirEdges=0;
255  int nlocalNonDirFaces=0;
256 
257  int edgematrixlocation=0;
258  int ntotaledgeentries=0;
259 
260  map<int,int> EdgeSize;
261  map<int,int> FaceSize;
262 
263  /// - Count edges, face and add up edges and face sizes
264  for(n = 0; n < n_exp; ++n)
265  {
266  eid = expList->GetOffset_Elmt_Id(n);
267  locExpansion = expList->GetExp(eid);
268 
269  nEdges = locExpansion->GetNedges();
270  for(j = 0; j < nEdges; ++j)
271  {
272  int nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
273  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
274  EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
275  }
276 
277  nFaces = locExpansion->GetNfaces();
278  for(j = 0; j < nFaces; ++j)
279  {
280  int nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
281  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
282  FaceSize[meshFaceId] = nFaceInteriorCoeffs;
283  }
284  }
285 
286  m_comm = expList->GetComm();
287 
288  // Loop over all the elements in the domain and compute max edge
289  // DOF and set up unique ordering.
290 
291  // First do periodic edges
292  PeriodicMap::const_iterator pIt;
293  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
294  {
295  meshEdgeId = pIt->first;
296 
297  if(edgeDirMap.count(meshEdgeId)==0)
298  {
299  dof = EdgeSize[meshEdgeId];
300  if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
301  {
302  bool SetUpNewEdge = true;
303 
304 
305  for (i = 0; i < pIt->second.size(); ++i)
306  {
307  if (!pIt->second[i].isLocal)
308  {
309  continue;
310  }
311 
312  int meshEdgeId2 = pIt->second[i].id;
313 
314  if(edgeDirMap.count(meshEdgeId2)==0)
315  {
316  if(uniqueEdgeMap.count(meshEdgeId2)!=0)
317  {
318  // set unique map to same location
319  uniqueEdgeMap[meshEdgeId] =
320  uniqueEdgeMap[meshEdgeId2];
321  SetUpNewEdge = false;
322  }
323  }
324  else
325  {
326  edgeDirMap.insert(meshEdgeId);
327  SetUpNewEdge = false;
328  }
329  }
330 
331  if(SetUpNewEdge)
332  {
333  uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
334 
335  edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
336 
337  edgemodeoffset[edgematrixlocation]=dof*dof;
338 
339  ntotaledgeentries+=dof*dof;
340 
341  n_blks[1+edgematrixlocation++]=dof;
342 
343  }
344  }
345  }
346  }
347 
348 
349  for(cnt=n=0; n < n_exp; ++n)
350  {
351  eid = expList->GetOffset_Elmt_Id(n);
352  locExpansion = expList->GetExp(eid);
353 
354  for (j = 0; j < locExpansion->GetNedges(); ++j)
355  {
356  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
357  dof = EdgeSize[meshEdgeId];
358  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
359 
360  if(edgeDirMap.count(meshEdgeId)==0)
361  {
362  if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
363 
364  {
365  uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
366 
367  edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
368 
369  edgemodeoffset[edgematrixlocation]=dof*dof;
370 
371  ntotaledgeentries+=dof*dof;
372 
373  n_blks[1+edgematrixlocation++]=dof;
374 
375  }
376  nlocalNonDirEdges+=dof*dof;
377  }
378  }
379  }
380 
381  int facematrixlocation=0;
382  int ntotalfaceentries=0;
383 
384  // Loop over all the elements in the domain and compute max face
385  // DOF. Reduce across all processes to get universal maximum.
386  // - Periodic faces
387  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
388  {
389  meshFaceId = pIt->first;
390 
391  if(faceDirMap.count(meshFaceId)==0)
392  {
393  dof = FaceSize[meshFaceId];
394 
395  if(uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
396  {
397  bool SetUpNewFace = true;
398 
399  if(pIt->second[0].isLocal)
400  {
401  int meshFaceId2 = pIt->second[0].id;
402 
403  if(faceDirMap.count(meshFaceId2)==0)
404  {
405  if(uniqueFaceMap.count(meshFaceId2)!=0)
406  {
407  // set unique map to same location
408  uniqueFaceMap[meshFaceId] =
409  uniqueFaceMap[meshFaceId2];
410  SetUpNewFace = false;
411  }
412  }
413  else // set face to be a Dirichlet face
414  {
415  faceDirMap.insert(meshFaceId);
416  SetUpNewFace = false;
417  }
418  }
419 
420  if(SetUpNewFace)
421  {
422  uniqueFaceMap[meshFaceId]=facematrixlocation;
423 
424  facemodeoffset[facematrixlocation]=dof*dof;
425 
426  faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
427 
428  ntotalfaceentries+=dof*dof;
429 
430  n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
431  }
432  }
433  }
434  }
435 
436 
437  for(cnt=n=0; n < n_exp; ++n)
438  {
439  eid = expList->GetOffset_Elmt_Id(n);
440 
441  locExpansion = expList->GetExp(eid);
442 
443  for (j = 0; j < locExpansion->GetNfaces(); ++j)
444  {
445  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
446 
447  dof = FaceSize[meshFaceId];
448  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
449 
450 
451  if(faceDirMap.count(meshFaceId)==0)
452  {
453  if(uniqueFaceMap.count(meshFaceId)==0 && dof > 0)
454  {
455  uniqueFaceMap[meshFaceId]=facematrixlocation;
456 
457  facemodeoffset[facematrixlocation]=dof*dof;
458 
459  faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
460 
461  ntotalfaceentries+=dof*dof;
462 
463  n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
464 
465  }
466  nlocalNonDirFaces+=dof*dof;
467  }
468  }
469  }
470 
471  m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
472  m_comm->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
473 
474  //Allocate arrays for block to universal map (number of expansions * p^2)
475  Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
476  Array<OneD, long> FaceBlockToUniversalMap(ntotalfaceentries,-1);
477 
478  Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
479  Array<OneD, int> localFaceToGlobalMatrixMap(nlocalNonDirFaces,-1);
480 
481  //Allocate arrays to store matrices (number of expansions * p^2)
482  Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
483  Array<OneD, NekDouble> FaceBlockArray(nlocalNonDirFaces,-1);
484 
485  int edgematrixoffset=0;
486  int facematrixoffset=0;
487  int vGlobal;
488 
489  for(n=0; n < n_exp; ++n)
490  {
491  eid = expList->GetOffset_Elmt_Id(n);
492 
493  locExpansion = expList->GetExp(eid);
494 
495  //loop over the edges of the expansion
496  for(j = 0; j < locExpansion->GetNedges(); ++j)
497  {
498  //get mesh edge id
499  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
500 
501  nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
502 
503  if(edgeDirMap.count(meshEdgeId)==0)
504  {
505  // Determine the Global edge offset
506  int edgeOffset = edgeglobaloffset[uniqueEdgeMap[meshEdgeId]];
507 
508  // Determine a universal map offset
509  int uniOffset = meshEdgeId;
510  pIt = periodicEdges.find(meshEdgeId);
511  if (pIt != periodicEdges.end())
512  {
513  for (int l = 0; l < pIt->second.size(); ++l)
514  {
515  uniOffset = min(uniOffset, pIt->second[l].id);
516  }
517  }
518  uniOffset = uniOffset *maxEdgeDof*maxEdgeDof;
519 
520  for(k=0; k<nedgemodes*nedgemodes; ++k)
521  {
522  vGlobal=edgeOffset+k;
523  localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
524  EdgeBlockToUniversalMap[vGlobal] = uniOffset + k + 1;
525  }
526  edgematrixoffset+=nedgemodes*nedgemodes;
527  }
528  }
529 
530  Array<OneD, unsigned int> faceInteriorMap;
531  Array<OneD, int> faceInteriorSign;
532  //loop over the faces of the expansion
533  for(j = 0; j < locExpansion->GetNfaces(); ++j)
534  {
535  //get mesh face id
536  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
537 
538  nfacemodes = locExpansion->GetFaceIntNcoeffs(j);
539 
540  //Check if face has dirichlet values
541  if(faceDirMap.count(meshFaceId)==0)
542  {
543  // Determine the Global edge offset
544  int faceOffset = faceglobaloffset[uniqueFaceMap[meshFaceId]];
545 
546  // Determine a universal map offset
547  int uniOffset = meshFaceId;
548  // use minimum face edge when periodic
549  pIt = periodicFaces.find(meshFaceId);
550  if (pIt != periodicFaces.end())
551  {
552  uniOffset = min(uniOffset, pIt->second[0].id);
553  }
554  uniOffset = uniOffset * maxFaceDof * maxFaceDof;
555 
556  for(k=0; k<nfacemodes*nfacemodes; ++k)
557  {
558  vGlobal=faceOffset+k;
559 
560  localFaceToGlobalMatrixMap[facematrixoffset+k]
561  = vGlobal;
562 
563  FaceBlockToUniversalMap[vGlobal] = uniOffset + k + 1;
564  }
565  facematrixoffset+=nfacemodes*nfacemodes;
566  }
567  }
568  }
569 
570  edgematrixoffset=0;
571  facematrixoffset=0;
572 
574  ::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
575 
576  const Array<OneD,const unsigned int>& nbdry_size
577  = m_locToGloMap->GetNumLocalBndCoeffsPerPatch();
578 
579  //Variants of R matrices required for low energy preconditioning
581  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
583  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
584 
585  //Here we loop over the expansion and build the block low energy
586  //preconditioner as well as the block versions of the transformation
587  //matrices.
588  for(cnt=n=0; n < n_exp; ++n)
589  {
590  eid = expList->GetOffset_Elmt_Id(n);
591 
592  locExpansion = expList->GetExp(eid);
593  nCoeffs=locExpansion->NumBndryCoeffs();
594  LibUtilities::ShapeType eType=locExpansion->DetShapeType();
595 
596  //Get correct transformation matrix for element type
597  R=(*(transmatrixmap[eType]));
598  RT=(*(transposedtransmatrixmap[eType]));
599 
601  (nCoeffs, nCoeffs, zero, storage);
602  RS = (*pRS);
603 
605  (nCoeffs, nCoeffs, zero, storage);
606  RSRT = (*pRSRT);
607 
608  nVerts=locExpansion->GetGeom()->GetNumVerts();
609  nEdges=locExpansion->GetGeom()->GetNumEdges();
610  nFaces=locExpansion->GetGeom()->GetNumFaces();
611 
612  //Get statically condensed matrix
613  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
614 
615  //Extract boundary block (elemental S1)
616  bnd_mat=loc_mat->GetBlock(0,0);
617 
618  //offset by number of rows
619  offset = bnd_mat->GetRows();
620 
621  DNekScalMat &S=(*bnd_mat);
622 
623  //Calculate S*trans(R) (note R is already transposed)
624  RS=R*S;
625 
626  //Calculate R*S*trans(R)
627  RSRT=RS*RT;
628 
629  //loop over vertices of the element and return the vertex map
630  //for each vertex
631  for (v=0; v<nVerts; ++v)
632  {
633  vMap1=locExpansion->GetVertexMap(v);
634 
635  //Get vertex map
636  globalrow = m_locToGloMap->
637  GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
638 
639  if(globalrow >= 0)
640  {
641  for (m=0; m<nVerts; ++m)
642  {
643  vMap2=locExpansion->GetVertexMap(m);
644 
645  //global matrix location (without offset due to
646  //dirichlet values)
647  globalcol = m_locToGloMap->
648  GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
649 
650  //offset for dirichlet conditions
651  if (globalcol == globalrow)
652  {
653  //modal connectivity between elements
654  sign1 = m_locToGloMap->
655  GetLocalToGlobalBndSign(cnt + vMap1);
656  sign2 = m_locToGloMap->
657  GetLocalToGlobalBndSign(cnt + vMap2);
658 
659  vertArray[globalrow]
660  += sign1*sign2*RSRT(vMap1,vMap2);
661 
662 
663  meshVertId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetVid(v);
664 
665  pIt = periodicVerts.find(meshVertId);
666  if (pIt != periodicVerts.end())
667  {
668  for (k = 0; k < pIt->second.size(); ++k)
669  {
670  meshVertId = min(meshVertId, pIt->second[k].id);
671  }
672  }
673 
674  VertBlockToUniversalMap[globalrow]
675  = meshVertId + 1;
676  }
677  }
678  }
679  }
680 
681  //loop over edges of the element and return the edge map
682  for (eid=0; eid<nEdges; ++eid)
683  {
684  nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
685 
686  DNekMatSharedPtr m_locMat =
688  (nedgemodes,nedgemodes,zero,storage);
689 
690  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(eid);
691  Array<OneD, unsigned int> edgemodearray = locExpansion->GetEdgeInverseBoundaryMap(eid);
692 
693  if(edgeDirMap.count(meshEdgeId)==0)
694  {
695  for (v=0; v<nedgemodes; ++v)
696  {
697  eMap1=edgemodearray[v];
698 
699  for (m=0; m<nedgemodes; ++m)
700  {
701  eMap2=edgemodearray[m];
702 
703  //modal connectivity between elements
704  sign1 = m_locToGloMap->
705  GetLocalToGlobalBndSign(cnt + eMap1);
706  sign2 = m_locToGloMap->
707  GetLocalToGlobalBndSign(cnt + eMap2);
708 
709  NekDouble globalEdgeValue = sign1*sign2*RSRT(eMap1,eMap2);
710 
711  EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=globalEdgeValue;
712  }
713  }
714  edgematrixoffset+=nedgemodes*nedgemodes;
715  }
716  }
717 
718  //loop over faces of the element and return the face map
719  for (fid=0; fid<nFaces; ++fid)
720  {
721  nfacemodes = locExpansion->GetFaceIntNcoeffs(fid);
722 
723  DNekMatSharedPtr m_locMat =
725  (nfacemodes,nfacemodes,zero,storage);
726 
727  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(fid);
728 
729  if(faceDirMap.count(meshFaceId)==0)
730  {
731  Array<OneD, unsigned int> facemodearray;
732  StdRegions::Orientation faceOrient = locExpansion->GetForient(fid);
733 
734  pIt = periodicFaces.find(meshFaceId);
735  if (pIt != periodicFaces.end())
736  {
737  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
738  {
739  facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid,faceOrient);
740  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
741  }
742  }
743 
744  facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid,faceOrient);
745 
746 
747  for (v=0; v<nfacemodes; ++v)
748  {
749  fMap1=facemodearray[v];
750 
751  for (m=0; m<nfacemodes; ++m)
752  {
753  fMap2=facemodearray[m];
754 
755  //modal connectivity between elements
756  sign1 = m_locToGloMap->
757  GetLocalToGlobalBndSign(cnt + fMap1);
758  sign2 = m_locToGloMap->
759  GetLocalToGlobalBndSign(cnt + fMap2);
760 
761  // Get the face-face value from the low energy matrix (S2)
762  NekDouble globalFaceValue = sign1*sign2*RSRT(fMap1,fMap2);
763 
764  //local face value to global face value
765  FaceBlockArray[facematrixoffset+v*nfacemodes+m]=globalFaceValue;
766  }
767  }
768  facematrixoffset+=nfacemodes*nfacemodes;
769  }
770  }
771 
772  //offset for the expansion
773  cnt+=offset;
774 
775  //Here we build the block matrices for R and RT
776  m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
777  m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
778  }
779 
780  if(nNonDirVerts!=0)
781  {
782  //Exchange vertex data over different processes
783  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
784  Gs::Gather(vertArray, Gs::gs_add, tmp);
785 
786  }
787 
788  Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries,0.0);
789  if(ntotaledgeentries)
790  {
791  //Assemble edge matrices of each process
792  Vmath::Assmb(EdgeBlockArray.num_elements(),
793  EdgeBlockArray,
794  localEdgeToGlobalMatrixMap,
795  GlobalEdgeBlock);
796  }
797 
798  //Exchange edge data over different processes
799  Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
800  Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
801 
802  Array<OneD, NekDouble> GlobalFaceBlock(ntotalfaceentries,0.0);
803  if(ntotalfaceentries)
804  {
805  //Assemble face matrices of each process
806  Vmath::Assmb(FaceBlockArray.num_elements(),
807  FaceBlockArray,
808  localFaceToGlobalMatrixMap,
809  GlobalFaceBlock);
810  }
811 
812  //Exchange face data over different processes
813  Gs::gs_data *tmp2 = Gs::Init(FaceBlockToUniversalMap, m_comm);
814  Gs::Gather(GlobalFaceBlock, Gs::gs_add, tmp2);
815 
816  // Populate vertex block
817  for (int i = 0; i < nNonDirVerts; ++i)
818  {
819  VertBlk->SetValue(i,i,1.0/vertArray[i]);
820  }
821 
822  //Set the first block to be the diagonal of the vertex space
823  m_BlkMat->SetBlock(0,0, VertBlk);
824 
825  offset=0;
826  //Build the edge matrices from the vector
827  DNekMatSharedPtr gmat;
828  for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
829  {
830  nedgemodes = n_blks[1+loc];
832  (nedgemodes,nedgemodes,zero,storage);
833 
834  for (v=0; v<nedgemodes; ++v)
835  {
836  for (m=0; m<nedgemodes; ++m)
837  {
838  NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
839  gmat->SetValue(v,m,EdgeValue);
840 
841  }
842  }
843 
844  m_BlkMat->SetBlock(1+loc,1+loc, gmat);
845  offset+=edgemodeoffset[loc];
846  }
847 
848  offset=0;
849 
850  Array<OneD, int> globalToUniversalMap = m_locToGloMap->GetGlobalToUniversalBndMap();
851  //Build the face matrices from the vector
852  for(int loc=0; loc<nNonDirFaceIDs; ++loc)
853  {
854  nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
856  (nfacemodes,nfacemodes,zero,storage);
857 
858  for (v=0; v<nfacemodes; ++v)
859  {
860  for (m=0; m<nfacemodes; ++m)
861  {
862  NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
863  gmat->SetValue(v,m,FaceValue);
864 
865  }
866  }
867  m_BlkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
868  offset+=facemodeoffset[loc];
869  }
870 
871  int totblks=m_BlkMat->GetNumberOfBlockRows();
872  for (i=1; i< totblks; ++i)
873  {
874  unsigned int nmodes=m_BlkMat->GetNumberOfRowsInBlockRow(i);
875  if(nmodes)
876  {
877  DNekMatSharedPtr tmp_mat =
879  (nmodes,nmodes,zero,storage);
880 
881  tmp_mat=m_BlkMat->GetBlock(i,i);
882  tmp_mat->Invert();
883 
884  m_BlkMat->SetBlock(i,i,tmp_mat);
885  }
886  }
887  }
888 
889 
890  /**
891  * Apply the low energy preconditioner during the conjugate gradient
892  * routine
893  */
895  const Array<OneD, NekDouble>& pInput,
896  Array<OneD, NekDouble>& pOutput)
897  {
898  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
899  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
900  int nNonDir = nGlobal-nDir;
901  DNekBlkMat &M = (*m_BlkMat);
902 
903  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
904  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
905 
906  z = M * r;
907  }
908 
909 
910  /**
911  * Set a block transformation matrices for each element type. These are
912  * needed in routines that transform the schur complement matrix to and
913  * from the low energy basis.
914  */
916  {
917  boost::shared_ptr<MultiRegions::ExpList>
918  expList=((m_linsys.lock())->GetLocMat()).lock();
920 
921  int n, nel;
922 
923  const Array<OneD,const unsigned int>& nbdry_size
924  = m_locToGloMap->GetNumLocalBndCoeffsPerPatch();
925 
926  int n_exp=expList->GetNumElmts();
927 
928  //maps for different element types
929  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
930  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
931  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransmatrixmap;
932  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransposedtransmatrixmap;
933 
934  //Transformation matrix map
935  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
936  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
937  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
938 
939  //Transposed transformation matrix map
940  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
941  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
942  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
943 
944  //Inverse transfomation map
945  invtransmatrixmap[LibUtilities::eTetrahedron]=m_Rinvtet;
946  invtransmatrixmap[LibUtilities::ePrism]=m_Rinvprism;
947  invtransmatrixmap[LibUtilities::eHexahedron]=m_Rinvhex;
948 
949  //Inverse transposed transformation map
950  invtransposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTinvtet;
951  invtransposedtransmatrixmap[LibUtilities::ePrism]=m_RTinvprism;
952  invtransposedtransmatrixmap[LibUtilities::eHexahedron]=m_RTinvhex;
953 
954  MatrixStorage blkmatStorage = eDIAGONAL;
955 
956  //Variants of R matrices required for low energy preconditioning
958  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
960  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
962  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
964  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
965 
966  for(n=0; n < n_exp; ++n)
967  {
968  nel = expList->GetOffset_Elmt_Id(n);
969 
970  locExpansion = expList->GetExp(nel);
971  LibUtilities::ShapeType eType=locExpansion->DetShapeType();
972 
973  //Block R matrix
974  m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
975 
976  //Block RT matrix
977  m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
978 
979  //Block inverse R matrix
980  m_InvRBlk->SetBlock(n,n, invtransmatrixmap[eType]);
981 
982  //Block inverse RT matrix
983  m_InvRTBlk->SetBlock(n,n, invtransposedtransmatrixmap[eType]);
984  }
985  }
986 
987 
988 
989  /**
990  * \brief Transform the solution vector vector to low energy.
991  *
992  * As the conjugate gradient system is solved for the low energy basis,
993  * the solution vector \f$\mathbf{x}\f$ must be transformed to the low
994  * energy basis i.e. \f$\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\f$.
995  */
997  Array<OneD, NekDouble>& pInOut,
998  int offset)
999  {
1000  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1001  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1002  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1003  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1004 
1005  //Non-dirichlet boundary dofs
1006  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pInOut+offset,
1007  eWrapper);
1008 
1009  //Block transformation matrix
1010  DNekScalBlkMat &R = *m_RBlk;
1011 
1012  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1013  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1014  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1015 
1016  //Not actually needed but we should only work with the Global boundary dofs
1017  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1018  Vmath::Vcopy(nGlobBndDofs, pInOut.get(), 1, tmp.get(), 1);
1019 
1020  //Global boundary (with dirichlet values) to local boundary with multiplicity
1021  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1022 
1023  //Multiply by the block transformation matrix
1024  F_LocBnd=R*F_LocBnd;
1025 
1026  //Assemble local boundary to global non-dirichlet Dofs
1027  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd, nDirBndDofs);
1028  }
1029 
1030  /**
1031  * \brief Transform the solution vector vector to low energy.
1032  *
1033  * As the conjugate gradient system is solved for the low energy basis,
1034  * the solution vector \f$\mathbf{x}\f$ must be transformed to the low
1035  * energy basis i.e. \f$\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\f$.
1036  */
1038  const Array<OneD, NekDouble>& pInput,
1039  Array<OneD, NekDouble>& pOutput)
1040  {
1041  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1042  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1043  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1044  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1045 
1046  //Input/output vectors should be length nGlobHomBndDofs
1047  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1048  "Input array is greater than the nGlobHomBndDofs");
1049  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1050  "Output array is greater than the nGlobHomBndDofs");
1051 
1052  //vectors of length number of non-dirichlet boundary dofs
1053  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1054  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1055  eWrapper);
1056  //Block transformation matrix
1057  DNekScalBlkMat &R = *m_RBlk;
1058 
1059  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1060  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1061  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1062 
1063  // Allocated array of size number of global boundary dofs and copy
1064  // the input array to the tmp array offset by Dirichlet boundary
1065  // conditions.
1066  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1067  Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1068 
1069  //Global boundary dofs (with zeroed dirichlet values) to local boundary dofs
1070  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1071 
1072  //Multiply by the block transformation matrix
1073  F_LocBnd=R*F_LocBnd;
1074 
1075  //Assemble local boundary to global non-dirichlet boundary
1076  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
1077  }
1078 
1079  /**
1080  * \brief transform the solution vector from low energy back to the
1081  * original basis.
1082  *
1083  * After the conjugate gradient routine the output vector is in the low
1084  * energy basis and must be trasnformed back to the original basis in
1085  * order to get the correct solution out. the solution vector
1086  * i.e. \f$\mathbf{x}=\mathbf{R^{T}}\mathbf{\overline{x}}\f$.
1087  */
1089  Array<OneD, NekDouble>& pInOut)
1090  {
1091  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1092  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1093  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1094  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1095 
1096  ASSERTL1(pInOut.num_elements() >= nGlobBndDofs,
1097  "Output array is greater than the nGlobBndDofs");
1098 
1099  //Block transposed transformation matrix
1100  DNekScalBlkMat &RT = *m_RTBlk;
1101 
1102  NekVector<NekDouble> V_GlobHomBnd(nGlobHomBndDofs,pInOut+nDirBndDofs,
1103  eWrapper);
1104 
1105  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1106  NekVector<NekDouble> V_LocBnd(nLocBndDofs,pLocal,eWrapper);
1107  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1108  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1109 
1110  //Global boundary (less dirichlet) to local boundary
1111  m_locToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
1112 
1113  //Multiply by the block transposed transformation matrix
1114  V_LocBnd=RT*V_LocBnd;
1115 
1116 
1117  //Assemble local boundary to global boundary
1118  Vmath::Assmb(nLocBndDofs, m_locToGloSignMult.get(),pLocal.get(), m_map.get(), tmp.get());
1119 
1120  //Universal assemble across processors
1121  m_locToGloMap->UniversalAssembleBnd(tmp);
1122 
1123  //copy non-dirichlet boundary values
1124  Vmath::Vcopy(nGlobBndDofs-nDirBndDofs, tmp.get() + nDirBndDofs, 1, pInOut.get() + nDirBndDofs, 1);
1125  }
1126 
1127  /**
1128  * \brief Multiply by the block inverse transformation matrix
1129  */
1131  const Array<OneD, NekDouble>& pInput,
1132  Array<OneD, NekDouble>& pOutput)
1133  {
1134  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1135  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1136  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1137  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1138 
1139  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1140  "Input array is greater than the nGlobHomBndDofs");
1141  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1142  "Output array is greater than the nGlobHomBndDofs");
1143 
1144  //vectors of length number of non-dirichlet boundary dofs
1145  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1146  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1147  eWrapper);
1148  //Block inverse transformation matrix
1149  DNekScalBlkMat &invR = *m_InvRBlk;
1150 
1151  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1152  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1153  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1154 
1155  // Allocated array of size number of global boundary dofs and copy
1156  // the input array to the tmp array offset by Dirichlet boundary
1157  // conditions.
1158  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1159  Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1160 
1161  //Global boundary dofs (with zeroed dirichlet values) to local boundary dofs
1162  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1163 
1164  //Multiply by block inverse transformation matrix
1165  F_LocBnd=invR*F_LocBnd;
1166 
1167  //Assemble local boundary to global non-dirichlet boundary
1168  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
1169 
1170  }
1171 
1172  /**
1173  * \brief Multiply by the block tranposed inverse transformation matrix
1174  */
1176  const Array<OneD, NekDouble>& pInput,
1177  Array<OneD, NekDouble>& pOutput)
1178  {
1179  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1180  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1181  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1182  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1183 
1184  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1185  "Input array is greater than the nGlobHomBndDofs");
1186  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1187  "Output array is greater than the nGlobHomBndDofs");
1188 
1189  //vectors of length number of non-dirichlet boundary dofs
1190  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1191  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1192  eWrapper);
1193  //Block inverse transformation matrix
1194  DNekScalBlkMat &invRT = *m_InvRTBlk;
1195 
1196  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1197  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1198  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1199 
1200  m_locToGloMap->GlobalToLocalBnd(pInput,pLocal, nDirBndDofs);
1201 
1202  //Multiply by the block transposed transformation matrix
1203  F_LocBnd=invRT*F_LocBnd;
1204 
1205  m_locToGloMap->AssembleBnd(pLocal,pOutput, nDirBndDofs);
1206 
1207  Vmath::Vmul(nGlobHomBndDofs,pOutput,1,m_multiplicity,1,pOutput,1);
1208  }
1209 
1210 
1211  /**
1212  * \brief Set up the transformed block matrix system
1213  *
1214  * Sets up a block elemental matrix in which each of the block matrix is
1215  * the low energy equivalent
1216  * i.e. \f$\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\f$
1217  */
1220  int offset,
1221  const boost::shared_ptr<DNekScalMat > &loc_mat)
1222  {
1223  boost::shared_ptr<MultiRegions::ExpList>
1224  expList=((m_linsys.lock())->GetLocMat()).lock();
1225 
1226  StdRegions::StdExpansionSharedPtr locExpansion;
1227  locExpansion = expList->GetExp(offset);
1228  unsigned int nbnd=locExpansion->NumBndryCoeffs();
1229 
1230  //This is the SC elemental matrix in the orginal basis (S1)
1231  DNekScalMatSharedPtr pS1=loc_mat;
1232 
1233  //Transformation matrices
1234  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
1235  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
1236  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
1237  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
1238  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
1239  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
1240  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
1241  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
1242 
1243  DNekScalMat &S1 = (*pS1);
1244 
1245  MatrixStorage storage = eFULL;
1246 
1247  DNekMatSharedPtr pS2 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1248  DNekMatSharedPtr pRS1 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1249 
1251  (expList->GetExp(offset))->DetShapeType();
1252 
1253  //transformation matrices
1254  DNekScalMat &R = (*(transmatrixmap[eType]));
1255  DNekScalMat &RT = (*(transposedtransmatrixmap[eType]));
1256 
1257  //create low energy matrix
1258  DNekMat &RS1 = (*pRS1);
1259  DNekMat &S2 = (*pS2);
1260 
1261  //setup S2
1262  RS1=R*S1;
1263  S2=RS1*RT;
1264 
1265  DNekScalMatSharedPtr tmp_mat;
1267 
1268  return tmp_mat;
1269  }
1270 
1271  /**
1272  * Create the inverse multiplicity map.
1273  */
1275  {
1276  unsigned int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
1277  unsigned int nEntries = m_locToGloMap->GetNumLocalBndCoeffs();
1278  unsigned int i;
1279 
1280  const Array<OneD, const int> &vMap
1281  = m_locToGloMap->GetLocalToGlobalBndMap();
1282 
1284  = m_locToGloMap->GetLocalToGlobalBndSign();
1285 
1286  bool m_signChange=m_locToGloMap->GetSignChange();
1287 
1288  // Count the multiplicity of each global DOF on this process
1289  Array<OneD, NekDouble> vCounts(nGlobalBnd, 0.0);
1290  for (i = 0; i < nEntries; ++i)
1291  {
1292  vCounts[vMap[i]] += 1.0;
1293  }
1294 
1295  // Get universal multiplicity by globally assembling counts
1296  m_locToGloMap->UniversalAssembleBnd(vCounts);
1297 
1298  // Construct a map of 1/multiplicity
1300  for (i = 0; i < nEntries; ++i)
1301  {
1302  if(m_signChange)
1303  {
1304  m_locToGloSignMult[i] = sign[i]*1.0/vCounts[vMap[i]];
1305  }
1306  else
1307  {
1308  m_locToGloSignMult[i] = 1.0/vCounts[vMap[i]];
1309  }
1310  }
1311 
1312  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1313  int nGlobHomBnd = nGlobalBnd - nDirBnd;
1314  int nLocBnd = m_locToGloMap->GetNumLocalBndCoeffs();
1315 
1316  //Set up multiplicity array for inverse transposed transformation matrix
1317  Array<OneD,NekDouble> tmp(nGlobHomBnd,1.0);
1318  m_multiplicity = Array<OneD,NekDouble>(nGlobHomBnd,1.0);
1319  Array<OneD,NekDouble> loc(nLocBnd,1.0);
1320 
1321  m_locToGloMap->GlobalToLocalBnd(tmp,loc, nDirBnd);
1322  m_locToGloMap->AssembleBnd(loc,m_multiplicity, nDirBnd);
1323  Vmath::Sdiv(nGlobHomBnd,1.0,m_multiplicity,1,m_multiplicity,1);
1324 
1325  }
1326 
1327  /**
1328  *\brief Sets up the reference prismatic element needed to construct
1329  *a low energy basis
1330  */
1332  {
1333  //////////////////////////
1334  // Set up Prism element //
1335  //////////////////////////
1336 
1337  const int three=3;
1338  const int nVerts = 6;
1339  const double point[][3] = {
1340  {-1,-1,0}, {1,-1,0}, {1,1,0},
1341  {-1,1,0}, {0,-1,sqrt(double(3))}, {0,1,sqrt(double(3))},
1342  };
1343 
1344  //boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
1346  for(int i=0; i < nVerts; ++i)
1347  {
1349  ( three, i, point[i][0], point[i][1], point[i][2] );
1350  }
1351  const int nEdges = 9;
1352  const int vertexConnectivity[][2] = {
1353  {0,1}, {1,2}, {3,2}, {0,3}, {0,4},
1354  {1,4}, {2,5}, {3,5}, {4,5}
1355  };
1356 
1357  // Populate the list of edges
1358  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1359  for(int i=0; i < nEdges; ++i){
1360  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1361  for(int j=0; j<2; ++j)
1362  {
1363  vertsArray[j] = verts[vertexConnectivity[i][j]];
1364  }
1365  edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(i, three, vertsArray);
1366  }
1367 
1368  ////////////////////////
1369  // Set up Prism faces //
1370  ////////////////////////
1371 
1372  const int nFaces = 5;
1373  //quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
1374  const int quadEdgeConnectivity[][4] = { {0,1,2,3}, {1,6,8,5}, {3,7,8,4} };
1375  const bool isQuadEdgeFlipped[][4] = { {0,0,1,1}, {0,0,1,1}, {0,0,1,1} };
1376  // QuadId ordered as 0, 1, 2, otherwise return false
1377  const int quadId[] = { 0,-1,1,-1,2 };
1378 
1379  //triangle-edge connectivity side-triface-1, side triface-3
1380  const int triEdgeConnectivity[][3] = { {0,5,4}, {2,6,7} };
1381  const bool isTriEdgeFlipped[][3] = { {0,0,1}, {0,0,1} };
1382  // TriId ordered as 0, 1, otherwise return false
1383  const int triId[] = { -1,0,-1,1,-1 };
1384 
1385  // Populate the list of faces
1386  SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1387  for(int f = 0; f < nFaces; ++f){
1388  if(f == 1 || f == 3) {
1389  int i = triId[f];
1390  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1391  StdRegions::Orientation eorientArray[3];
1392  for(int j = 0; j < 3; ++j){
1393  edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1394  eorientArray[j] = isTriEdgeFlipped[i][j] ? StdRegions::eBackwards : StdRegions::eForwards;
1395  }
1396  faces[f] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(f, edgeArray, eorientArray);
1397  }
1398  else {
1399  int i = quadId[f];
1400  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1401  StdRegions::Orientation eorientArray[4];
1402  for(int j=0; j < 4; ++j){
1403  edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1404  eorientArray[j] = isQuadEdgeFlipped[i][j] ? StdRegions::eBackwards : StdRegions::eForwards;
1405  }
1406  faces[f] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(f, edgeArray, eorientArray);
1407  }
1408  }
1409 
1411 
1412  geom->SetOwnData();
1413 
1414  return geom;
1415  }
1416 
1417  /**
1418  *\brief Sets up the reference tretrahedral element needed to construct
1419  *a low energy basis
1420  */
1422  {
1423  /////////////////////////////////
1424  // Set up Tetrahedron vertices //
1425  /////////////////////////////////
1426 
1427  int i,j;
1428  const int three=3;
1429  const int nVerts = 4;
1430  const double point[][3] = {
1431  {-1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1432  {1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1433  {0,2/sqrt(double(3)),-1/sqrt(double(6))},
1434  {0,0,3/sqrt(double(6))}};
1435 
1436  boost::shared_ptr<SpatialDomains::PointGeom> verts[4];
1437  for(i=0; i < nVerts; ++i)
1438  {
1439  verts[i] =
1442  ( three, i, point[i][0], point[i][1], point[i][2] );
1443  }
1444 
1445  //////////////////////////////
1446  // Set up Tetrahedron Edges //
1447  //////////////////////////////
1448 
1449  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1450  const int nEdges = 6;
1451  const int vertexConnectivity[][2] = {
1452  {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1453  };
1454 
1455  // Populate the list of edges
1456  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1457  for(i=0; i < nEdges; ++i)
1458  {
1459  boost::shared_ptr<SpatialDomains::PointGeom>
1460  vertsArray[2];
1461  for(j=0; j<2; ++j)
1462  {
1463  vertsArray[j] = verts[vertexConnectivity[i][j]];
1464  }
1465 
1467  ::AllocateSharedPtr(i, three, vertsArray);
1468  }
1469 
1470  //////////////////////////////
1471  // Set up Tetrahedron faces //
1472  //////////////////////////////
1473 
1474  const int nFaces = 4;
1475  const int edgeConnectivity[][3] = {
1476  {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1477  };
1478  const bool isEdgeFlipped[][3] = {
1479  {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
1480  };
1481 
1482  // Populate the list of faces
1483  SpatialDomains::TriGeomSharedPtr faces[nFaces];
1484  for(i=0; i < nFaces; ++i)
1485  {
1486  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1487  StdRegions::Orientation eorientArray[3];
1488  for(j=0; j < 3; ++j)
1489  {
1490  edgeArray[j] = edges[edgeConnectivity[i][j]];
1491  eorientArray[j] = isEdgeFlipped[i][j] ?
1493  }
1494 
1495 
1497  ::AllocateSharedPtr(i, edgeArray, eorientArray);
1498  }
1499 
1502  (faces);
1503 
1504  geom->SetOwnData();
1505 
1506  return geom;
1507  }
1508 
1509  /**
1510  *\brief Sets up the reference hexahedral element needed to construct
1511  *a low energy basis
1512  */
1514  {
1515  ////////////////////////////////
1516  // Set up Hexahedron vertices //
1517  ////////////////////////////////
1518 
1519  const int three=3;
1520 
1521  const int nVerts = 8;
1522  const double point[][3] = {
1523  {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1524  {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1525  };
1526 
1527  // Populate the list of verts
1529  for( int i = 0; i < nVerts; ++i ) {
1531  ::AllocateSharedPtr(three, i, point[i][0],
1532  point[i][1], point[i][2]);
1533  }
1534 
1535  /////////////////////////////
1536  // Set up Hexahedron Edges //
1537  /////////////////////////////
1538 
1539  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1540  const int nEdges = 12;
1541  const int vertexConnectivity[][2] = {
1542  {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1543  {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1544  };
1545 
1546  // Populate the list of edges
1547  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1548  for( int i = 0; i < nEdges; ++i ) {
1549  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1550  for( int j = 0; j < 2; ++j ) {
1551  vertsArray[j] = verts[vertexConnectivity[i][j]];
1552  }
1554  AllocateSharedPtr( i, three, vertsArray);
1555  }
1556 
1557  /////////////////////////////
1558  // Set up Hexahedron faces //
1559  /////////////////////////////
1560 
1561  const int nFaces = 6;
1562  const int edgeConnectivity[][4] = {
1563  {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1564  {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1565  };
1566  const bool isEdgeFlipped[][4] = {
1567  {0,0,0,1}, {0,0,1,1}, {0,0,1,1},
1568  {0,0,1,1}, {0,0,1,1}, {0,0,0,1}
1569  };
1570 
1571  // Populate the list of faces
1572  SpatialDomains::QuadGeomSharedPtr faces[nFaces];
1573  for( int i = 0; i < nFaces; ++i ) {
1574  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1575  StdRegions::Orientation eorientArray[4];
1576  for( int j = 0; j < 4; ++j ) {
1577  edgeArray[j] = edges[edgeConnectivity[i][j]];
1578  eorientArray[j] = isEdgeFlipped[i][j] ?
1580  }
1582  eorientArray);
1583  }
1584 
1587  (faces);
1588 
1589  geom->SetOwnData();
1590 
1591  return geom;
1592  }
1593 
1594 
1595  /**
1596  * \brief Sets up the reference elements needed by the preconditioner
1597  *
1598  * Sets up reference elements which are used to preconditioning the
1599  * corresponding matrices. Currently we support tetrahedral, prismatic
1600  * and hexahedral elements
1601  */
1603  {
1604  int cnt,i,j;
1605  boost::shared_ptr<MultiRegions::ExpList>
1606  expList=((m_linsys.lock())->GetLocMat()).lock();
1607  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
1608  StdRegions::VarCoeffMap vVarCoeffMap;
1609  StdRegions::StdExpansionSharedPtr locExpansion;
1610  locExpansion = expList->GetExp(0);
1611 
1612  DNekScalBlkMatSharedPtr RtetBlk, RprismBlk;
1613  DNekScalBlkMatSharedPtr RTtetBlk, RTprismBlk;
1614 
1615  DNekScalMatSharedPtr Rprismoriginal;
1616  DNekScalMatSharedPtr RTprismoriginal;
1617  DNekMatSharedPtr Rtettmp, RTtettmp, Rhextmp, RThextmp, Rprismtmp, RTprismtmp ;
1618 
1619  /*
1620  * Set up a Tetrahral & prismatic element which comprises
1621  * equilateral triangles as all faces for the tet and the end faces
1622  * for the prism. Using these elements a new expansion is created
1623  * (which is the same as the expansion specified in the input
1624  * file).
1625  */
1629 
1630  //Expansion as specified in the input file - here we need to alter
1631  //this so we can read in different exapansions for different element
1632  //types
1633  int nummodes=locExpansion->GetBasisNumModes(0);
1634 
1635  //Bases for Tetrahedral element
1636  const LibUtilities::BasisKey TetBa(
1637  LibUtilities::eModified_A, nummodes,
1639  const LibUtilities::BasisKey TetBb(
1640  LibUtilities::eModified_B, nummodes,
1642  const LibUtilities::BasisKey TetBc(
1643  LibUtilities::eModified_C, nummodes,
1645 
1646  //Create reference tetrahedral expansion
1648 
1650  ::AllocateSharedPtr(TetBa,TetBb,TetBc,
1651  tetgeom);
1652 
1653  //Bases for prismatic element
1654  const LibUtilities::BasisKey PrismBa(
1655  LibUtilities::eModified_A, nummodes,
1657  const LibUtilities::BasisKey PrismBb(
1658  LibUtilities::eModified_A, nummodes,
1660  const LibUtilities::BasisKey PrismBc(
1661  LibUtilities::eModified_B, nummodes,
1663 
1664  //Create reference prismatic expansion
1666 
1668  ::AllocateSharedPtr(PrismBa,PrismBb,PrismBc,
1669  prismgeom);
1670 
1671  //Bases for prismatic element
1672  const LibUtilities::BasisKey HexBa(
1673  LibUtilities::eModified_A, nummodes,
1675  const LibUtilities::BasisKey HexBb(
1676  LibUtilities::eModified_A, nummodes,
1678  const LibUtilities::BasisKey HexBc(
1679  LibUtilities::eModified_A, nummodes,
1681 
1682  //Create reference prismatic expansion
1684 
1686  ::AllocateSharedPtr(HexBa,HexBb,HexBc,
1687  hexgeom);
1688 
1689 
1690  // retrieve variable coefficient
1691  if(m_linSysKey.GetNVarCoeffs() > 0)
1692  {
1693  StdRegions::VarCoeffMap::const_iterator x;
1694  cnt = expList->GetPhys_Offset(0);
1695  for (x = m_linSysKey.GetVarCoeffs().begin();
1696  x != m_linSysKey.GetVarCoeffs().end(); ++x)
1697  {
1698  vVarCoeffMap[x->first] = x->second + cnt;
1699  }
1700  }
1701 
1702  StdRegions::MatrixType PreconR,PreconRT;
1703 
1704  if(m_linSysKey.GetMatrixType() == StdRegions::eMass)
1705  {
1706  PreconR = StdRegions::ePreconRMass;
1707  PreconRT = StdRegions::ePreconRTMass;
1708  }
1709  else
1710  {
1711  PreconR = StdRegions::ePreconR;
1712  PreconRT = StdRegions::ePreconRT;
1713  }
1714 
1715 
1716 
1717  /*
1718  * Matrix keys - for each element type there are two matrix keys
1719  * corresponding to the transformation matrix R and its transpose
1720  */
1721 
1722 
1723  //Matrix keys for tetrahedral element transformation matrix
1725  (PreconR, LibUtilities::eTetrahedron,
1726  *TetExp, m_linSysKey.GetConstFactors(),
1727  vVarCoeffMap);
1728 
1729  //Matrix keys for tetrahedral transposed transformation matrix
1731  (PreconRT, LibUtilities::eTetrahedron,
1732  *TetExp, m_linSysKey.GetConstFactors(),
1733  vVarCoeffMap);
1734 
1735  //Matrix keys for prismatic element transformation matrix
1737  (PreconR, LibUtilities::ePrism,
1738  *PrismExp, m_linSysKey.GetConstFactors(),
1739  vVarCoeffMap);
1740 
1741  //Matrix keys for prismatic element transposed transformation matrix
1742  LocalRegions::MatrixKey PrismRT
1743  (PreconRT, LibUtilities::ePrism,
1744  *PrismExp, m_linSysKey.GetConstFactors(),
1745  vVarCoeffMap);
1746 
1747  //Matrix keys for hexahedral element transformation matrix
1749  (PreconR, LibUtilities::eHexahedron,
1750  *HexExp, m_linSysKey.GetConstFactors(),
1751  vVarCoeffMap);
1752 
1753  //Matrix keys for hexahedral element transposed transformation
1754  //matrix
1756  (PreconRT, LibUtilities::eHexahedron,
1757  *HexExp, m_linSysKey.GetConstFactors(),
1758  vVarCoeffMap);
1759 
1760  /*
1761  * Create transformation matrices for the tetrahedral element
1762  */
1763 
1764  //Get tetrahedral transformation matrix
1765  m_Rtet = TetExp->GetLocMatrix(TetR);
1766 
1767  //Get tetrahedral transposed transformation matrix
1768  m_RTtet = TetExp->GetLocMatrix(TetRT);
1769 
1770  // Using the transformation matrix and the inverse transformation
1771  // matrix create the inverse matrices
1772  Rtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1773 
1774  //Inverse transposed transformation matrix
1775  RTtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1776  RTtettmp->Transpose();
1777 
1779  ::AllocateSharedPtr(1.0,Rtettmp);
1781  ::AllocateSharedPtr(1.0,RTtettmp);
1782 
1783  /*
1784  * Create transformation matrices for the hexahedral element
1785  */
1786 
1787  //Get hexahedral transformation matrix
1788  m_Rhex = HexExp->GetLocMatrix(HexR);
1789  //Get hexahedral transposed transformation matrix
1790  m_RThex = HexExp->GetLocMatrix(HexRT);
1791 
1792  // Using the transformation matrix and the inverse transformation
1793  // matrix create the inverse matrices
1794  Rhextmp=HexExp->BuildInverseTransformationMatrix(m_Rhex);
1795  //Inverse transposed transformation matrix
1796  RThextmp=HexExp->BuildInverseTransformationMatrix(m_Rhex);
1797  RThextmp->Transpose();
1798 
1800  ::AllocateSharedPtr(1.0,Rhextmp);
1802  ::AllocateSharedPtr(1.0,RThextmp);
1803 
1804  /*
1805  * Create transformation matrices for the prismatic element
1806  */
1807 
1808  //Get prism transformation matrix
1809  Rprismoriginal = PrismExp->GetLocMatrix(PrismR);
1810  //Get prism transposed transformation matrix
1811  RTprismoriginal = PrismExp->GetLocMatrix(PrismRT);
1812 
1813  unsigned int nRows=Rprismoriginal->GetRows();
1814  NekDouble zero=0.0;
1816  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1818  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1819  NekDouble Rvalue, RTvalue;
1820 
1821  //Copy values from the prism transformation matrix
1822  for(i=0; i<nRows; ++i)
1823  {
1824  for(j=0; j<nRows; ++j)
1825  {
1826  Rvalue=(*Rprismoriginal)(i,j);
1827  RTvalue=(*RTprismoriginal)(i,j);
1828  Rtmpprism->SetValue(i,j,Rvalue);
1829  RTtmpprism->SetValue(i,j,RTvalue);
1830  }
1831  }
1832 
1833  //Replace triangular faces and edges of the prims transformation
1834  //matrix with the corresponding values of the tetrahedral
1835  //transformation matrix.
1836  ModifyPrismTransformationMatrix(TetExp,PrismExp,Rtmpprism,RTtmpprism);
1837 
1839  ::AllocateSharedPtr(1.0,Rtmpprism);
1840 
1842  ::AllocateSharedPtr(1.0,RTtmpprism);
1843 
1844  //Inverse transformation matrix
1845  Rprismtmp=PrismExp->BuildInverseTransformationMatrix(m_Rprism);
1846 
1847  //Inverse transposed transformation matrix
1848  RTprismtmp=PrismExp->BuildInverseTransformationMatrix(m_Rprism);
1849  RTprismtmp->Transpose();
1850 
1852  ::AllocateSharedPtr(1.0,Rprismtmp);
1853 
1855  ::AllocateSharedPtr(1.0,RTprismtmp);
1856  }
1857 
1858  /**
1859  * \brief Modify the prism transformation matrix to align with the
1860  * tetrahedral modes.
1861  *
1862  * This routine replaces the edge and triangular face components of the
1863  * prismatic vertex transformation matrices \f$\mathbf{R}_{ve}\f$ and
1864  * \f$\mathbf{R}_{vf}\f$ with the corresponding components from the
1865  * tetrahedral transformation matrices. Additionally, triangular face
1866  * components in the prismatic edge transformation matrix
1867  * \f$\mathbf{R}_{ef}\f$ with the corresponding component from the
1868  * tetrahedral transformation matrix.
1869  */
1873  DNekMatSharedPtr Rmodprism,
1874  DNekMatSharedPtr RTmodprism)
1875  {
1876  NekDouble Rvalue, RTvalue;
1877  int i, j;
1878 
1879  //For a tet element the bottom face is made up of the following:
1880  //vertices: 0, 1 and 2 edges: 0, 1 and 2 face: 0. We first need to
1881  //determine the mode locations of these vertices, edges and face so
1882  //we can extract the correct values from the tetrahedral R matrix.
1883 
1884  //These are the vertex mode locations of R which need to be replaced
1885  //in the prism element
1886  int TetVertex0=TetExp->GetVertexMap(0);
1887  int TetVertex1=TetExp->GetVertexMap(1);
1888  int TetVertex2=TetExp->GetVertexMap(2);
1889  int TetVertex3=TetExp->GetVertexMap(3);
1890 
1891 
1892  //These are the edge mode locations of R which need to be replaced
1893  //in the prism element
1894  Array<OneD, unsigned int> TetEdge0=TetExp->GetEdgeInverseBoundaryMap(0);
1895  Array<OneD, unsigned int> TetEdge1=TetExp->GetEdgeInverseBoundaryMap(1);
1896  Array<OneD, unsigned int> TetEdge2=TetExp->GetEdgeInverseBoundaryMap(2);
1897  Array<OneD, unsigned int> TetEdge3=TetExp->GetEdgeInverseBoundaryMap(3);
1898  Array<OneD, unsigned int> TetEdge4=TetExp->GetEdgeInverseBoundaryMap(4);
1899  Array<OneD, unsigned int> TetEdge5=TetExp->GetEdgeInverseBoundaryMap(5);
1900 
1901  //These are the face mode locations of R which need to be replaced
1902  //in the prism element
1903  Array<OneD, unsigned int> TetFace=TetExp->GetFaceInverseBoundaryMap(1);
1904 
1905  //Prism vertex modes
1906  int PrismVertex0=PrismExp->GetVertexMap(0);
1907  int PrismVertex1=PrismExp->GetVertexMap(1);
1908  int PrismVertex2=PrismExp->GetVertexMap(2);
1909  int PrismVertex3=PrismExp->GetVertexMap(3);
1910  int PrismVertex4=PrismExp->GetVertexMap(4);
1911  int PrismVertex5=PrismExp->GetVertexMap(5);
1912 
1913  //Prism edge modes
1914  Array<OneD, unsigned int> PrismEdge0=
1915  PrismExp->GetEdgeInverseBoundaryMap(0);
1916  Array<OneD, unsigned int> PrismEdge1=
1917  PrismExp->GetEdgeInverseBoundaryMap(1);
1918  Array<OneD, unsigned int> PrismEdge2=
1919  PrismExp->GetEdgeInverseBoundaryMap(2);
1920  Array<OneD, unsigned int> PrismEdge3=
1921  PrismExp->GetEdgeInverseBoundaryMap(3);
1922  Array<OneD, unsigned int> PrismEdge4=
1923  PrismExp->GetEdgeInverseBoundaryMap(4);
1924  Array<OneD, unsigned int> PrismEdge5=
1925  PrismExp->GetEdgeInverseBoundaryMap(5);
1926  Array<OneD, unsigned int> PrismEdge6=
1927  PrismExp->GetEdgeInverseBoundaryMap(6);
1928  Array<OneD, unsigned int> PrismEdge7=
1929  PrismExp->GetEdgeInverseBoundaryMap(7);
1930  Array<OneD, unsigned int> PrismEdge8=
1931  PrismExp->GetEdgeInverseBoundaryMap(8);
1932 
1933  //Prism face 1 & 3 face modes
1934  Array<OneD, unsigned int> PrismFace1=
1935  PrismExp->GetFaceInverseBoundaryMap(1);
1936  Array<OneD, unsigned int> PrismFace3=
1937  PrismExp->GetFaceInverseBoundaryMap(3);
1938  Array<OneD, unsigned int> PrismFace0=
1939  PrismExp->GetFaceInverseBoundaryMap(0);
1940  Array<OneD, unsigned int> PrismFace2=
1941  PrismExp->GetFaceInverseBoundaryMap(2);
1942  Array<OneD, unsigned int> PrismFace4=
1943  PrismExp->GetFaceInverseBoundaryMap(4);
1944 
1945  //vertex 0 edge 0 3 & 4
1946  for(i=0; i< PrismEdge0.num_elements(); ++i)
1947  {
1948  Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
1949  Rmodprism->SetValue(PrismVertex0,PrismEdge0[i],Rvalue);
1950  Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
1951  Rmodprism->SetValue(PrismVertex0,PrismEdge3[i],Rvalue);
1952  Rvalue=(*m_Rtet)(TetVertex0,TetEdge3[i]);
1953  Rmodprism->SetValue(PrismVertex0,PrismEdge4[i],Rvalue);
1954 
1955  //transposed values
1956  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
1957  RTmodprism->SetValue(PrismEdge0[i],PrismVertex0,RTvalue);
1958  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
1959  RTmodprism->SetValue(PrismEdge3[i],PrismVertex0,RTvalue);
1960  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex0);
1961  RTmodprism->SetValue(PrismEdge4[i],PrismVertex0,RTvalue);
1962  }
1963 
1964  //vertex 1 edge 0 1 & 5
1965  for(i=0; i< PrismEdge1.num_elements(); ++i)
1966  {
1967  Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1968  Rmodprism->SetValue(PrismVertex1,PrismEdge0[i],Rvalue);
1969  Rvalue=(*m_Rtet)(TetVertex1,TetEdge1[i]);
1970  Rmodprism->SetValue(PrismVertex1,PrismEdge1[i],Rvalue);
1971  Rvalue=(*m_Rtet)(TetVertex1,TetEdge4[i]);
1972  Rmodprism->SetValue(PrismVertex1,PrismEdge5[i],Rvalue);
1973 
1974  //transposed values
1975  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1976  RTmodprism->SetValue(PrismEdge0[i],PrismVertex1,RTvalue);
1977  RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex1);
1978  RTmodprism->SetValue(PrismEdge1[i],PrismVertex1,RTvalue);
1979  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex1);
1980  RTmodprism->SetValue(PrismEdge5[i],PrismVertex1,RTvalue);
1981  }
1982 
1983  //vertex 2 edge 1 2 & 6
1984  for(i=0; i< PrismEdge2.num_elements(); ++i)
1985  {
1986  Rvalue=(*m_Rtet)(TetVertex2,TetEdge1[i]);
1987  Rmodprism->SetValue(PrismVertex2,PrismEdge1[i],Rvalue);
1988  Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1989  Rmodprism->SetValue(PrismVertex2,PrismEdge2[i],Rvalue);
1990  Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
1991  Rmodprism->SetValue(PrismVertex2,PrismEdge6[i],Rvalue);
1992 
1993  //transposed values
1994  RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex2);
1995  RTmodprism->SetValue(PrismEdge1[i],PrismVertex2,RTvalue);
1996  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1997  RTmodprism->SetValue(PrismEdge2[i],PrismVertex2,RTvalue);
1998  RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
1999  RTmodprism->SetValue(PrismEdge6[i],PrismVertex2,RTvalue);
2000  }
2001 
2002  //vertex 3 edge 3 2 & 7
2003  for(i=0; i< PrismEdge3.num_elements(); ++i)
2004  {
2005  Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2006  Rmodprism->SetValue(PrismVertex3,PrismEdge3[i],Rvalue);
2007  Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
2008  Rmodprism->SetValue(PrismVertex3,PrismEdge2[i],Rvalue);
2009  Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
2010  Rmodprism->SetValue(PrismVertex3,PrismEdge7[i],Rvalue);
2011 
2012  //transposed values
2013  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2014  RTmodprism->SetValue(PrismEdge3[i],PrismVertex3,RTvalue);
2015  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
2016  RTmodprism->SetValue(PrismEdge2[i],PrismVertex3,RTvalue);
2017  RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2018  RTmodprism->SetValue(PrismEdge7[i],PrismVertex3,RTvalue);
2019  }
2020 
2021  //vertex 4 edge 4 5 & 8
2022  for(i=0; i< PrismEdge4.num_elements(); ++i)
2023  {
2024  Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2025  Rmodprism->SetValue(PrismVertex4,PrismEdge4[i],Rvalue);
2026  Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2027  Rmodprism->SetValue(PrismVertex4,PrismEdge5[i],Rvalue);
2028  Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
2029  Rmodprism->SetValue(PrismVertex4,PrismEdge8[i],Rvalue);
2030 
2031  //transposed values
2032  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2033  RTmodprism->SetValue(PrismEdge4[i],PrismVertex4,RTvalue);
2034  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2035  RTmodprism->SetValue(PrismEdge5[i],PrismVertex4,RTvalue);
2036  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
2037  RTmodprism->SetValue(PrismEdge8[i],PrismVertex4,RTvalue);
2038  }
2039 
2040  //vertex 5 edge 6 7 & 8
2041  for(i=0; i< PrismEdge5.num_elements(); ++i)
2042  {
2043  Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2044  Rmodprism->SetValue(PrismVertex5,PrismEdge6[i],Rvalue);
2045  Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2046  Rmodprism->SetValue(PrismVertex5,PrismEdge7[i],Rvalue);
2047  Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2048  Rmodprism->SetValue(PrismVertex5,PrismEdge8[i],Rvalue);
2049 
2050  //transposed values
2051  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2052  RTmodprism->SetValue(PrismEdge6[i],PrismVertex5,RTvalue);
2053  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2054  RTmodprism->SetValue(PrismEdge7[i],PrismVertex5,RTvalue);
2055  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2056  RTmodprism->SetValue(PrismEdge8[i],PrismVertex5,RTvalue);
2057  }
2058 
2059  // face 1 vertices 0 1 4
2060  for(i=0; i< PrismFace1.num_elements(); ++i)
2061  {
2062  Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2063  Rmodprism->SetValue(PrismVertex0,PrismFace1[i],Rvalue);
2064  Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2065  Rmodprism->SetValue(PrismVertex1,PrismFace1[i],Rvalue);
2066  Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2067  Rmodprism->SetValue(PrismVertex4,PrismFace1[i],Rvalue);
2068 
2069  //transposed values
2070  RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2071  RTmodprism->SetValue(PrismFace1[i],PrismVertex0,RTvalue);
2072  RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2073  RTmodprism->SetValue(PrismFace1[i],PrismVertex1,RTvalue);
2074  RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2075  RTmodprism->SetValue(PrismFace1[i],PrismVertex4,RTvalue);
2076  }
2077 
2078  // face 3 vertices 2, 3 & 5
2079  for(i=0; i< PrismFace3.num_elements(); ++i)
2080  {
2081  Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2082  Rmodprism->SetValue(PrismVertex2,PrismFace3[i],Rvalue);
2083  Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2084  Rmodprism->SetValue(PrismVertex3,PrismFace3[i],Rvalue);
2085  Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2086  Rmodprism->SetValue(PrismVertex5,PrismFace3[i],Rvalue);
2087 
2088  //transposed values
2089  RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2090  RTmodprism->SetValue(PrismFace3[i],PrismVertex2,RTvalue);
2091  RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2092  RTmodprism->SetValue(PrismFace3[i],PrismVertex3,RTvalue);
2093  RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2094  RTmodprism->SetValue(PrismFace3[i],PrismVertex5,RTvalue);
2095  }
2096 
2097  // Face 1 edge 0 4 5
2098  for(i=0; i< PrismFace1.num_elements(); ++i)
2099  {
2100  for(j=0; j<PrismEdge0.num_elements(); ++j)
2101  {
2102  Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2103  Rmodprism->SetValue(PrismEdge0[j],PrismFace1[i],Rvalue);
2104  Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2105  Rmodprism->SetValue(PrismEdge4[j],PrismFace1[i],Rvalue);
2106  Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2107  Rmodprism->SetValue(PrismEdge5[j],PrismFace1[i],Rvalue);
2108 
2109  //transposed values
2110  RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2111  RTmodprism->SetValue(PrismFace1[i],PrismEdge0[j],RTvalue);
2112  RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2113  RTmodprism->SetValue(PrismFace1[i],PrismEdge4[j],RTvalue);
2114  RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2115  RTmodprism->SetValue(PrismFace1[i],PrismEdge5[j],RTvalue);
2116  }
2117  }
2118 
2119  // Face 3 edge 2 6 7
2120  for(i=0; i< PrismFace3.num_elements(); ++i)
2121  {
2122  for(j=0; j<PrismEdge2.num_elements(); ++j)
2123  {
2124  Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2125  Rmodprism->SetValue(PrismEdge2[j],PrismFace3[i],Rvalue);
2126  Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2127  Rmodprism->SetValue(PrismEdge6[j],PrismFace3[i],Rvalue);
2128  Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2129  Rmodprism->SetValue(PrismEdge7[j],PrismFace3[i],Rvalue);
2130 
2131  RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2132  RTmodprism->SetValue(PrismFace3[i],PrismEdge2[j],RTvalue);
2133  RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2134  RTmodprism->SetValue(PrismFace3[i],PrismEdge6[j],RTvalue);
2135  RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2136  RTmodprism->SetValue(PrismFace3[i],PrismEdge7[j],RTvalue);
2137  }
2138  }
2139  }
2140 
2141  }
2142 }
2143 
2144 
2145 
2146 
2147 
2148 
const StdRegions::VarCoeffMap & GetVarCoeffs() const
LibUtilities::CommSharedPtr m_comm
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Principle Modified Functions .
Definition: BasisType.h:51
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
virtual void v_DoMultiplybyInverseTransposedTransformationMatrix(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Multiply by the block tranposed inverse transformation matrix.
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
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
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
Principle Modified Functions .
Definition: BasisType.h:49
PreconFactory & GetPreconFactory()
STL namespace.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:257
boost::shared_ptr< HexExp > HexExpSharedPtr
Definition: HexExp.h:57
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion.h:287
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< TetExp > TetExpSharedPtr
Definition: TetExp.h:223
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
virtual DNekScalMatSharedPtr v_TransformedSchurCompl(int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
Set up the transformed block matrix system.
virtual void v_BuildPreconditioner()
Construct the low energy preconditioner from .
int GetFaceIntNcoeffs(const int i) const
Definition: StdExpansion.h:359
void SetUpReferenceElements(void)
Sets up the reference elements needed by the preconditioner.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:226
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Principle Modified Functions .
Definition: BasisType.h:50
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
virtual void v_DoTransformToLowEnergy(Array< OneD, NekDouble > &pInOut, int offset)
Transform the solution vector vector to low energy.
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:695
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< PrismExp > PrismExpSharedPtr
Definition: PrismExp.h:217
boost::shared_ptr< T > as()
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_DoMultiplybyInverseTransformationMatrix(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Multiply by the block inverse transformation matrix.
Describe a linear system.
virtual void v_DoTransformFromLowEnergy(Array< OneD, NekDouble > &pInOut)
transform the solution vector from low energy back to the original basis.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom(void)
Sets up the reference tretrahedral element needed to construct a low energy basis.
void ModifyPrismTransformationMatrix(LocalRegions::TetExpSharedPtr TetExp, LocalRegions::PrismExpSharedPtr PrismExp, DNekMatSharedPtr Rmodprism, DNekMatSharedPtr RTmodprism)
Modify the prism transformation matrix to align with the tetrahedral modes.
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis.
const boost::weak_ptr< GlobalLinSys > m_linsys
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Describes the specification for a Basis.
Definition: Basis.h:50
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215