Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 = 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 = 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 = 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 = 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 = 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  bool verbose =
781  expList->GetSession()->DefinesCmdLineArgument("verbose");
782 
783  if(nNonDirVerts!=0)
784  {
785  //Exchange vertex data over different processes
786  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm, verbose);
787  Gs::Gather(vertArray, Gs::gs_add, tmp);
788 
789  }
790 
791  Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries,0.0);
792  if(ntotaledgeentries)
793  {
794  //Assemble edge matrices of each process
795  Vmath::Assmb(EdgeBlockArray.num_elements(),
796  EdgeBlockArray,
797  localEdgeToGlobalMatrixMap,
798  GlobalEdgeBlock);
799  }
800 
801  //Exchange edge data over different processes
802  Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm, verbose);
803  Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
804 
805  Array<OneD, NekDouble> GlobalFaceBlock(ntotalfaceentries,0.0);
806  if(ntotalfaceentries)
807  {
808  //Assemble face matrices of each process
809  Vmath::Assmb(FaceBlockArray.num_elements(),
810  FaceBlockArray,
811  localFaceToGlobalMatrixMap,
812  GlobalFaceBlock);
813  }
814 
815  //Exchange face data over different processes
816  Gs::gs_data *tmp2 = Gs::Init(FaceBlockToUniversalMap, m_comm, verbose);
817  Gs::Gather(GlobalFaceBlock, Gs::gs_add, tmp2);
818 
819  // Populate vertex block
820  for (int i = 0; i < nNonDirVerts; ++i)
821  {
822  VertBlk->SetValue(i,i,1.0/vertArray[i]);
823  }
824 
825  //Set the first block to be the diagonal of the vertex space
826  m_BlkMat->SetBlock(0,0, VertBlk);
827 
828  offset=0;
829  //Build the edge matrices from the vector
830  DNekMatSharedPtr gmat;
831  for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
832  {
833  nedgemodes = n_blks[1+loc];
835  (nedgemodes,nedgemodes,zero,storage);
836 
837  for (v=0; v<nedgemodes; ++v)
838  {
839  for (m=0; m<nedgemodes; ++m)
840  {
841  NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
842  gmat->SetValue(v,m,EdgeValue);
843 
844  }
845  }
846 
847  m_BlkMat->SetBlock(1+loc,1+loc, gmat);
848  offset+=edgemodeoffset[loc];
849  }
850 
851  offset=0;
852 
853  Array<OneD, int> globalToUniversalMap = m_locToGloMap->GetGlobalToUniversalBndMap();
854  //Build the face matrices from the vector
855  for(int loc=0; loc<nNonDirFaceIDs; ++loc)
856  {
857  nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
859  (nfacemodes,nfacemodes,zero,storage);
860 
861  for (v=0; v<nfacemodes; ++v)
862  {
863  for (m=0; m<nfacemodes; ++m)
864  {
865  NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
866  gmat->SetValue(v,m,FaceValue);
867 
868  }
869  }
870  m_BlkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
871  offset+=facemodeoffset[loc];
872  }
873 
874  int totblks=m_BlkMat->GetNumberOfBlockRows();
875  for (i=1; i< totblks; ++i)
876  {
877  unsigned int nmodes=m_BlkMat->GetNumberOfRowsInBlockRow(i);
878  if(nmodes)
879  {
880  DNekMatSharedPtr tmp_mat =
882  (nmodes,nmodes,zero,storage);
883 
884  tmp_mat=m_BlkMat->GetBlock(i,i);
885  tmp_mat->Invert();
886 
887  m_BlkMat->SetBlock(i,i,tmp_mat);
888  }
889  }
890  }
891 
892 
893  /**
894  * Apply the low energy preconditioner during the conjugate gradient
895  * routine
896  */
898  const Array<OneD, NekDouble>& pInput,
899  Array<OneD, NekDouble>& pOutput)
900  {
901  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
902  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
903  int nNonDir = nGlobal-nDir;
904  DNekBlkMat &M = (*m_BlkMat);
905 
906  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
907  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
908 
909  z = M * r;
910  }
911 
912 
913  /**
914  * Set a block transformation matrices for each element type. These are
915  * needed in routines that transform the schur complement matrix to and
916  * from the low energy basis.
917  */
919  {
920  boost::shared_ptr<MultiRegions::ExpList>
921  expList=((m_linsys.lock())->GetLocMat()).lock();
923 
924  int n, nel;
925 
926  const Array<OneD,const unsigned int>& nbdry_size
927  = m_locToGloMap->GetNumLocalBndCoeffsPerPatch();
928 
929  int n_exp=expList->GetNumElmts();
930 
931  //maps for different element types
932  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
933  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
934  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransmatrixmap;
935  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> invtransposedtransmatrixmap;
936 
937  //Transformation matrix map
938  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
939  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
940  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
941 
942  //Transposed transformation matrix map
943  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
944  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
945  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
946 
947  //Inverse transfomation map
948  invtransmatrixmap[LibUtilities::eTetrahedron]=m_Rinvtet;
949  invtransmatrixmap[LibUtilities::ePrism]=m_Rinvprism;
950  invtransmatrixmap[LibUtilities::eHexahedron]=m_Rinvhex;
951 
952  //Inverse transposed transformation map
953  invtransposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTinvtet;
954  invtransposedtransmatrixmap[LibUtilities::ePrism]=m_RTinvprism;
955  invtransposedtransmatrixmap[LibUtilities::eHexahedron]=m_RTinvhex;
956 
957  MatrixStorage blkmatStorage = eDIAGONAL;
958 
959  //Variants of R matrices required for low energy preconditioning
961  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
963  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
965  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
967  ::AllocateSharedPtr(nbdry_size, nbdry_size , blkmatStorage);
968 
969  for(n=0; n < n_exp; ++n)
970  {
971  nel = n;
972 
973  locExpansion = expList->GetExp(nel);
974  LibUtilities::ShapeType eType=locExpansion->DetShapeType();
975 
976  //Block R matrix
977  m_RBlk->SetBlock(n,n, transmatrixmap[eType]);
978 
979  //Block RT matrix
980  m_RTBlk->SetBlock(n,n, transposedtransmatrixmap[eType]);
981 
982  //Block inverse R matrix
983  m_InvRBlk->SetBlock(n,n, invtransmatrixmap[eType]);
984 
985  //Block inverse RT matrix
986  m_InvRTBlk->SetBlock(n,n, invtransposedtransmatrixmap[eType]);
987  }
988  }
989 
990 
991 
992  /**
993  * \brief Transform the solution vector vector to low energy.
994  *
995  * As the conjugate gradient system is solved for the low energy basis,
996  * the solution vector \f$\mathbf{x}\f$ must be transformed to the low
997  * energy basis i.e. \f$\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\f$.
998  */
1000  Array<OneD, NekDouble>& pInOut,
1001  int offset)
1002  {
1003  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1004  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1005  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1006  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1007 
1008  //Non-dirichlet boundary dofs
1009  Array<OneD, NekDouble> tmpOffset = pInOut + offset;
1010  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs, tmpOffset, eWrapper);
1011 
1012  //Block transformation matrix
1013  DNekScalBlkMat &R = *m_RBlk;
1014 
1015  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1016  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1017  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1018 
1019  //Not actually needed but we should only work with the Global boundary dofs
1020  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1021  Vmath::Vcopy(nGlobBndDofs, pInOut.get(), 1, tmp.get(), 1);
1022 
1023  //Global boundary (with dirichlet values) to local boundary with multiplicity
1024  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1025 
1026  //Multiply by the block transformation matrix
1027  F_LocBnd=R*F_LocBnd;
1028 
1029  //Assemble local boundary to global non-dirichlet Dofs
1030  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd, nDirBndDofs);
1031  }
1032 
1033  /**
1034  * \brief Transform the solution vector vector to low energy.
1035  *
1036  * As the conjugate gradient system is solved for the low energy basis,
1037  * the solution vector \f$\mathbf{x}\f$ must be transformed to the low
1038  * energy basis i.e. \f$\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\f$.
1039  */
1041  const Array<OneD, NekDouble>& pInput,
1042  Array<OneD, NekDouble>& pOutput)
1043  {
1044  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1045  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1046  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1047  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1048 
1049  //Input/output vectors should be length nGlobHomBndDofs
1050  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1051  "Input array is greater than the nGlobHomBndDofs");
1052  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1053  "Output array is greater than the nGlobHomBndDofs");
1054 
1055  //vectors of length number of non-dirichlet boundary dofs
1056  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1057  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1058  eWrapper);
1059  //Block transformation matrix
1060  DNekScalBlkMat &R = *m_RBlk;
1061 
1062  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1063  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1064  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1065 
1066  // Allocated array of size number of global boundary dofs and copy
1067  // the input array to the tmp array offset by Dirichlet boundary
1068  // conditions.
1069  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1070  Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1071 
1072  //Global boundary dofs (with zeroed dirichlet values) to local boundary dofs
1073  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1074 
1075  //Multiply by the block transformation matrix
1076  F_LocBnd=R*F_LocBnd;
1077 
1078  //Assemble local boundary to global non-dirichlet boundary
1079  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
1080  }
1081 
1082  /**
1083  * \brief transform the solution vector from low energy back to the
1084  * original basis.
1085  *
1086  * After the conjugate gradient routine the output vector is in the low
1087  * energy basis and must be trasnformed back to the original basis in
1088  * order to get the correct solution out. the solution vector
1089  * i.e. \f$\mathbf{x}=\mathbf{R^{T}}\mathbf{\overline{x}}\f$.
1090  */
1092  Array<OneD, NekDouble>& pInOut)
1093  {
1094  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1095  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1096  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1097  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1098 
1099  ASSERTL1(pInOut.num_elements() >= nGlobBndDofs,
1100  "Output array is greater than the nGlobBndDofs");
1101 
1102  //Block transposed transformation matrix
1103  DNekScalBlkMat &RT = *m_RTBlk;
1104 
1105  Array<OneD, NekDouble> tmpOffset = pInOut + nDirBndDofs;
1106  NekVector<NekDouble> V_GlobHomBnd(nGlobHomBndDofs, tmpOffset, eWrapper);
1107 
1108  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1109  NekVector<NekDouble> V_LocBnd(nLocBndDofs,pLocal,eWrapper);
1110  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1111  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1112 
1113  //Global boundary (less dirichlet) to local boundary
1114  m_locToGloMap->GlobalToLocalBnd(V_GlobHomBnd,V_LocBnd, nDirBndDofs);
1115 
1116  //Multiply by the block transposed transformation matrix
1117  V_LocBnd=RT*V_LocBnd;
1118 
1119 
1120  //Assemble local boundary to global boundary
1121  Vmath::Assmb(nLocBndDofs, m_locToGloSignMult.get(),pLocal.get(), m_map.get(), tmp.get());
1122 
1123  //Universal assemble across processors
1124  m_locToGloMap->UniversalAssembleBnd(tmp);
1125 
1126  //copy non-dirichlet boundary values
1127  Vmath::Vcopy(nGlobBndDofs-nDirBndDofs, tmp.get() + nDirBndDofs, 1, pInOut.get() + nDirBndDofs, 1);
1128  }
1129 
1130  /**
1131  * \brief Multiply by the block inverse transformation matrix
1132  */
1134  const Array<OneD, NekDouble>& pInput,
1135  Array<OneD, NekDouble>& pOutput)
1136  {
1137  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1138  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1139  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1140  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1141 
1142  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1143  "Input array is greater than the nGlobHomBndDofs");
1144  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1145  "Output array is greater than the nGlobHomBndDofs");
1146 
1147  //vectors of length number of non-dirichlet boundary dofs
1148  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1149  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1150  eWrapper);
1151  //Block inverse transformation matrix
1152  DNekScalBlkMat &invR = *m_InvRBlk;
1153 
1154  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1155  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1156  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1157 
1158  // Allocated array of size number of global boundary dofs and copy
1159  // the input array to the tmp array offset by Dirichlet boundary
1160  // conditions.
1161  Array<OneD,NekDouble> tmp(nGlobBndDofs,0.0);
1162  Vmath::Vcopy(nGlobHomBndDofs, pInput.get(), 1, tmp.get() + nDirBndDofs, 1);
1163 
1164  //Global boundary dofs (with zeroed dirichlet values) to local boundary dofs
1165  Vmath::Gathr(m_map.num_elements(), m_locToGloSignMult.get(), tmp.get(), m_map.get(), pLocal.get());
1166 
1167  //Multiply by block inverse transformation matrix
1168  F_LocBnd=invR*F_LocBnd;
1169 
1170  //Assemble local boundary to global non-dirichlet boundary
1171  m_locToGloMap->AssembleBnd(F_LocBnd,F_HomBnd,nDirBndDofs);
1172 
1173  }
1174 
1175  /**
1176  * \brief Multiply by the block tranposed inverse transformation matrix
1177  */
1179  const Array<OneD, NekDouble>& pInput,
1180  Array<OneD, NekDouble>& pOutput)
1181  {
1182  int nGlobBndDofs = m_locToGloMap->GetNumGlobalBndCoeffs();
1183  int nDirBndDofs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1184  int nGlobHomBndDofs = nGlobBndDofs - nDirBndDofs;
1185  int nLocBndDofs = m_locToGloMap->GetNumLocalBndCoeffs();
1186 
1187  ASSERTL1(pInput.num_elements() >= nGlobHomBndDofs,
1188  "Input array is greater than the nGlobHomBndDofs");
1189  ASSERTL1(pOutput.num_elements() >= nGlobHomBndDofs,
1190  "Output array is greater than the nGlobHomBndDofs");
1191 
1192  //vectors of length number of non-dirichlet boundary dofs
1193  NekVector<NekDouble> F_GlobBnd(nGlobHomBndDofs,pInput,eWrapper);
1194  NekVector<NekDouble> F_HomBnd(nGlobHomBndDofs,pOutput,
1195  eWrapper);
1196  //Block inverse transformation matrix
1197  DNekScalBlkMat &invRT = *m_InvRTBlk;
1198 
1199  Array<OneD, NekDouble> pLocal(nLocBndDofs, 0.0);
1200  NekVector<NekDouble> F_LocBnd(nLocBndDofs,pLocal,eWrapper);
1201  m_map = m_locToGloMap->GetLocalToGlobalBndMap();
1202 
1203  m_locToGloMap->GlobalToLocalBnd(pInput,pLocal, nDirBndDofs);
1204 
1205  //Multiply by the block transposed transformation matrix
1206  F_LocBnd=invRT*F_LocBnd;
1207 
1208  m_locToGloMap->AssembleBnd(pLocal,pOutput, nDirBndDofs);
1209 
1210  Vmath::Vmul(nGlobHomBndDofs,pOutput,1,m_multiplicity,1,pOutput,1);
1211  }
1212 
1213 
1214  /**
1215  * \brief Set up the transformed block matrix system
1216  *
1217  * Sets up a block elemental matrix in which each of the block matrix is
1218  * the low energy equivalent
1219  * i.e. \f$\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\f$
1220  */
1223  int offset,
1224  const boost::shared_ptr<DNekScalMat > &loc_mat)
1225  {
1226  boost::shared_ptr<MultiRegions::ExpList>
1227  expList=((m_linsys.lock())->GetLocMat()).lock();
1228 
1229  StdRegions::StdExpansionSharedPtr locExpansion;
1230  locExpansion = expList->GetExp(offset);
1231  unsigned int nbnd=locExpansion->NumBndryCoeffs();
1232 
1233  //This is the SC elemental matrix in the orginal basis (S1)
1234  DNekScalMatSharedPtr pS1=loc_mat;
1235 
1236  //Transformation matrices
1237  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transmatrixmap;
1238  map<LibUtilities::ShapeType,DNekScalMatSharedPtr> transposedtransmatrixmap;
1239  transmatrixmap[LibUtilities::eTetrahedron]=m_Rtet;
1240  transmatrixmap[LibUtilities::ePrism]=m_Rprism;
1241  transmatrixmap[LibUtilities::eHexahedron]=m_Rhex;
1242  transposedtransmatrixmap[LibUtilities::eTetrahedron]=m_RTtet;
1243  transposedtransmatrixmap[LibUtilities::ePrism]=m_RTprism;
1244  transposedtransmatrixmap[LibUtilities::eHexahedron]=m_RThex;
1245 
1246  DNekScalMat &S1 = (*pS1);
1247 
1248  MatrixStorage storage = eFULL;
1249 
1250  DNekMatSharedPtr pS2 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1251  DNekMatSharedPtr pRS1 = MemoryManager<DNekMat>::AllocateSharedPtr(nbnd,nbnd,0.0,storage);
1252 
1254  (expList->GetExp(offset))->DetShapeType();
1255 
1256  //transformation matrices
1257  DNekScalMat &R = (*(transmatrixmap[eType]));
1258  DNekScalMat &RT = (*(transposedtransmatrixmap[eType]));
1259 
1260  //create low energy matrix
1261  DNekMat &RS1 = (*pRS1);
1262  DNekMat &S2 = (*pS2);
1263 
1264  //setup S2
1265  RS1=R*S1;
1266  S2=RS1*RT;
1267 
1268  DNekScalMatSharedPtr tmp_mat;
1270 
1271  return tmp_mat;
1272  }
1273 
1274  /**
1275  * Create the inverse multiplicity map.
1276  */
1278  {
1279  unsigned int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
1280  unsigned int nEntries = m_locToGloMap->GetNumLocalBndCoeffs();
1281  unsigned int i;
1282 
1283  const Array<OneD, const int> &vMap
1284  = m_locToGloMap->GetLocalToGlobalBndMap();
1285 
1287  = m_locToGloMap->GetLocalToGlobalBndSign();
1288 
1289  bool m_signChange=m_locToGloMap->GetSignChange();
1290 
1291  // Count the multiplicity of each global DOF on this process
1292  Array<OneD, NekDouble> vCounts(nGlobalBnd, 0.0);
1293  for (i = 0; i < nEntries; ++i)
1294  {
1295  vCounts[vMap[i]] += 1.0;
1296  }
1297 
1298  // Get universal multiplicity by globally assembling counts
1299  m_locToGloMap->UniversalAssembleBnd(vCounts);
1300 
1301  // Construct a map of 1/multiplicity
1303  for (i = 0; i < nEntries; ++i)
1304  {
1305  if(m_signChange)
1306  {
1307  m_locToGloSignMult[i] = sign[i]*1.0/vCounts[vMap[i]];
1308  }
1309  else
1310  {
1311  m_locToGloSignMult[i] = 1.0/vCounts[vMap[i]];
1312  }
1313  }
1314 
1315  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1316  int nGlobHomBnd = nGlobalBnd - nDirBnd;
1317  int nLocBnd = m_locToGloMap->GetNumLocalBndCoeffs();
1318 
1319  //Set up multiplicity array for inverse transposed transformation matrix
1320  Array<OneD,NekDouble> tmp(nGlobHomBnd,1.0);
1321  m_multiplicity = Array<OneD,NekDouble>(nGlobHomBnd,1.0);
1322  Array<OneD,NekDouble> loc(nLocBnd,1.0);
1323 
1324  m_locToGloMap->GlobalToLocalBnd(tmp,loc, nDirBnd);
1325  m_locToGloMap->AssembleBnd(loc,m_multiplicity, nDirBnd);
1326  Vmath::Sdiv(nGlobHomBnd,1.0,m_multiplicity,1,m_multiplicity,1);
1327 
1328  }
1329 
1330  /**
1331  *\brief Sets up the reference prismatic element needed to construct
1332  *a low energy basis
1333  */
1335  {
1336  //////////////////////////
1337  // Set up Prism element //
1338  //////////////////////////
1339 
1340  const int three=3;
1341  const int nVerts = 6;
1342  const double point[][3] = {
1343  {-1,-1,0}, {1,-1,0}, {1,1,0},
1344  {-1,1,0}, {0,-1,sqrt(double(3))}, {0,1,sqrt(double(3))},
1345  };
1346 
1347  //boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
1349  for(int i=0; i < nVerts; ++i)
1350  {
1352  ( three, i, point[i][0], point[i][1], point[i][2] );
1353  }
1354  const int nEdges = 9;
1355  const int vertexConnectivity[][2] = {
1356  {0,1}, {1,2}, {3,2}, {0,3}, {0,4},
1357  {1,4}, {2,5}, {3,5}, {4,5}
1358  };
1359 
1360  // Populate the list of edges
1361  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1362  for(int i=0; i < nEdges; ++i){
1363  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1364  for(int j=0; j<2; ++j)
1365  {
1366  vertsArray[j] = verts[vertexConnectivity[i][j]];
1367  }
1368  edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(i, three, vertsArray);
1369  }
1370 
1371  ////////////////////////
1372  // Set up Prism faces //
1373  ////////////////////////
1374 
1375  const int nFaces = 5;
1376  //quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
1377  const int quadEdgeConnectivity[][4] = { {0,1,2,3}, {1,6,8,5}, {3,7,8,4} };
1378  const bool isQuadEdgeFlipped[][4] = { {0,0,1,1}, {0,0,1,1}, {0,0,1,1} };
1379  // QuadId ordered as 0, 1, 2, otherwise return false
1380  const int quadId[] = { 0,-1,1,-1,2 };
1381 
1382  //triangle-edge connectivity side-triface-1, side triface-3
1383  const int triEdgeConnectivity[][3] = { {0,5,4}, {2,6,7} };
1384  const bool isTriEdgeFlipped[][3] = { {0,0,1}, {0,0,1} };
1385  // TriId ordered as 0, 1, otherwise return false
1386  const int triId[] = { -1,0,-1,1,-1 };
1387 
1388  // Populate the list of faces
1389  SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1390  for(int f = 0; f < nFaces; ++f){
1391  if(f == 1 || f == 3) {
1392  int i = triId[f];
1393  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1394  StdRegions::Orientation eorientArray[3];
1395  for(int j = 0; j < 3; ++j){
1396  edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1397  eorientArray[j] = isTriEdgeFlipped[i][j] ? StdRegions::eBackwards : StdRegions::eForwards;
1398  }
1399  faces[f] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(f, edgeArray, eorientArray);
1400  }
1401  else {
1402  int i = quadId[f];
1403  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1404  StdRegions::Orientation eorientArray[4];
1405  for(int j=0; j < 4; ++j){
1406  edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1407  eorientArray[j] = isQuadEdgeFlipped[i][j] ? StdRegions::eBackwards : StdRegions::eForwards;
1408  }
1409  faces[f] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(f, edgeArray, eorientArray);
1410  }
1411  }
1412 
1414 
1415  geom->SetOwnData();
1416 
1417  return geom;
1418  }
1419 
1420  /**
1421  *\brief Sets up the reference tretrahedral element needed to construct
1422  *a low energy basis
1423  */
1425  {
1426  /////////////////////////////////
1427  // Set up Tetrahedron vertices //
1428  /////////////////////////////////
1429 
1430  int i,j;
1431  const int three=3;
1432  const int nVerts = 4;
1433  const double point[][3] = {
1434  {-1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1435  {1,-1/sqrt(double(3)),-1/sqrt(double(6))},
1436  {0,2/sqrt(double(3)),-1/sqrt(double(6))},
1437  {0,0,3/sqrt(double(6))}};
1438 
1439  boost::shared_ptr<SpatialDomains::PointGeom> verts[4];
1440  for(i=0; i < nVerts; ++i)
1441  {
1442  verts[i] =
1445  ( three, i, point[i][0], point[i][1], point[i][2] );
1446  }
1447 
1448  //////////////////////////////
1449  // Set up Tetrahedron Edges //
1450  //////////////////////////////
1451 
1452  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1453  const int nEdges = 6;
1454  const int vertexConnectivity[][2] = {
1455  {0,1},{1,2},{0,2},{0,3},{1,3},{2,3}
1456  };
1457 
1458  // Populate the list of edges
1459  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1460  for(i=0; i < nEdges; ++i)
1461  {
1462  boost::shared_ptr<SpatialDomains::PointGeom>
1463  vertsArray[2];
1464  for(j=0; j<2; ++j)
1465  {
1466  vertsArray[j] = verts[vertexConnectivity[i][j]];
1467  }
1468 
1470  ::AllocateSharedPtr(i, three, vertsArray);
1471  }
1472 
1473  //////////////////////////////
1474  // Set up Tetrahedron faces //
1475  //////////////////////////////
1476 
1477  const int nFaces = 4;
1478  const int edgeConnectivity[][3] = {
1479  {0,1,2}, {0,4,3}, {1,5,4}, {2,5,3}
1480  };
1481  const bool isEdgeFlipped[][3] = {
1482  {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
1483  };
1484 
1485  // Populate the list of faces
1486  SpatialDomains::TriGeomSharedPtr faces[nFaces];
1487  for(i=0; i < nFaces; ++i)
1488  {
1489  SpatialDomains::SegGeomSharedPtr edgeArray[3];
1490  StdRegions::Orientation eorientArray[3];
1491  for(j=0; j < 3; ++j)
1492  {
1493  edgeArray[j] = edges[edgeConnectivity[i][j]];
1494  eorientArray[j] = isEdgeFlipped[i][j] ?
1496  }
1497 
1498 
1500  ::AllocateSharedPtr(i, edgeArray, eorientArray);
1501  }
1502 
1505  (faces);
1506 
1507  geom->SetOwnData();
1508 
1509  return geom;
1510  }
1511 
1512  /**
1513  *\brief Sets up the reference hexahedral element needed to construct
1514  *a low energy basis
1515  */
1517  {
1518  ////////////////////////////////
1519  // Set up Hexahedron vertices //
1520  ////////////////////////////////
1521 
1522  const int three=3;
1523 
1524  const int nVerts = 8;
1525  const double point[][3] = {
1526  {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0},
1527  {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}
1528  };
1529 
1530  // Populate the list of verts
1532  for( int i = 0; i < nVerts; ++i ) {
1534  ::AllocateSharedPtr(three, i, point[i][0],
1535  point[i][1], point[i][2]);
1536  }
1537 
1538  /////////////////////////////
1539  // Set up Hexahedron Edges //
1540  /////////////////////////////
1541 
1542  // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1543  const int nEdges = 12;
1544  const int vertexConnectivity[][2] = {
1545  {0,1}, {1,2}, {2,3}, {0,3}, {0,4}, {1,5},
1546  {2,6}, {3,7}, {4,5}, {5,6}, {6,7}, {4,7}
1547  };
1548 
1549  // Populate the list of edges
1550  SpatialDomains::SegGeomSharedPtr edges[nEdges];
1551  for( int i = 0; i < nEdges; ++i ) {
1552  SpatialDomains::PointGeomSharedPtr vertsArray[2];
1553  for( int j = 0; j < 2; ++j ) {
1554  vertsArray[j] = verts[vertexConnectivity[i][j]];
1555  }
1557  AllocateSharedPtr( i, three, vertsArray);
1558  }
1559 
1560  /////////////////////////////
1561  // Set up Hexahedron faces //
1562  /////////////////////////////
1563 
1564  const int nFaces = 6;
1565  const int edgeConnectivity[][4] = {
1566  {0,1,2,3}, {0,5,8,4}, {1,6,9,5},
1567  {2,7,10,6}, {3,7,11,4}, {8,9,10,11}
1568  };
1569  const bool isEdgeFlipped[][4] = {
1570  {0,0,0,1}, {0,0,1,1}, {0,0,1,1},
1571  {0,0,1,1}, {0,0,1,1}, {0,0,0,1}
1572  };
1573 
1574  // Populate the list of faces
1575  SpatialDomains::QuadGeomSharedPtr faces[nFaces];
1576  for( int i = 0; i < nFaces; ++i ) {
1577  SpatialDomains::SegGeomSharedPtr edgeArray[4];
1578  StdRegions::Orientation eorientArray[4];
1579  for( int j = 0; j < 4; ++j ) {
1580  edgeArray[j] = edges[edgeConnectivity[i][j]];
1581  eorientArray[j] = isEdgeFlipped[i][j] ?
1583  }
1585  eorientArray);
1586  }
1587 
1590  (faces);
1591 
1592  geom->SetOwnData();
1593 
1594  return geom;
1595  }
1596 
1597 
1598  /**
1599  * \brief Sets up the reference elements needed by the preconditioner
1600  *
1601  * Sets up reference elements which are used to preconditioning the
1602  * corresponding matrices. Currently we support tetrahedral, prismatic
1603  * and hexahedral elements
1604  */
1606  {
1607  int cnt,i,j;
1608  boost::shared_ptr<MultiRegions::ExpList>
1609  expList=((m_linsys.lock())->GetLocMat()).lock();
1610  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
1611  StdRegions::VarCoeffMap vVarCoeffMap;
1612  StdRegions::StdExpansionSharedPtr locExpansion;
1613  locExpansion = expList->GetExp(0);
1614 
1615  DNekScalBlkMatSharedPtr RtetBlk, RprismBlk;
1616  DNekScalBlkMatSharedPtr RTtetBlk, RTprismBlk;
1617 
1618  DNekScalMatSharedPtr Rprismoriginal;
1619  DNekScalMatSharedPtr RTprismoriginal;
1620  DNekMatSharedPtr Rtettmp, RTtettmp, Rhextmp, RThextmp, Rprismtmp, RTprismtmp ;
1621 
1622  /*
1623  * Set up a Tetrahral & prismatic element which comprises
1624  * equilateral triangles as all faces for the tet and the end faces
1625  * for the prism. Using these elements a new expansion is created
1626  * (which is the same as the expansion specified in the input
1627  * file).
1628  */
1632 
1633  //Expansion as specified in the input file - here we need to alter
1634  //this so we can read in different exapansions for different element
1635  //types
1636  int nummodes=locExpansion->GetBasisNumModes(0);
1637 
1638  //Bases for Tetrahedral element
1639  const LibUtilities::BasisKey TetBa(
1640  LibUtilities::eModified_A, nummodes,
1642  const LibUtilities::BasisKey TetBb(
1643  LibUtilities::eModified_B, nummodes,
1645  const LibUtilities::BasisKey TetBc(
1646  LibUtilities::eModified_C, nummodes,
1648 
1649  //Create reference tetrahedral expansion
1651 
1653  ::AllocateSharedPtr(TetBa,TetBb,TetBc,
1654  tetgeom);
1655 
1656  //Bases for prismatic element
1657  const LibUtilities::BasisKey PrismBa(
1658  LibUtilities::eModified_A, nummodes,
1660  const LibUtilities::BasisKey PrismBb(
1661  LibUtilities::eModified_A, nummodes,
1663  const LibUtilities::BasisKey PrismBc(
1664  LibUtilities::eModified_B, nummodes,
1666 
1667  //Create reference prismatic expansion
1669 
1671  ::AllocateSharedPtr(PrismBa,PrismBb,PrismBc,
1672  prismgeom);
1673 
1674  //Bases for prismatic element
1675  const LibUtilities::BasisKey HexBa(
1676  LibUtilities::eModified_A, nummodes,
1678  const LibUtilities::BasisKey HexBb(
1679  LibUtilities::eModified_A, nummodes,
1681  const LibUtilities::BasisKey HexBc(
1682  LibUtilities::eModified_A, nummodes,
1684 
1685  //Create reference prismatic expansion
1687 
1689  ::AllocateSharedPtr(HexBa,HexBb,HexBc,
1690  hexgeom);
1691 
1692 
1693  // retrieve variable coefficient
1694  if(m_linSysKey.GetNVarCoeffs() > 0)
1695  {
1696  StdRegions::VarCoeffMap::const_iterator x;
1697  cnt = expList->GetPhys_Offset(0);
1698  for (x = m_linSysKey.GetVarCoeffs().begin();
1699  x != m_linSysKey.GetVarCoeffs().end(); ++x)
1700  {
1701  vVarCoeffMap[x->first] = x->second + cnt;
1702  }
1703  }
1704 
1705  StdRegions::MatrixType PreconR,PreconRT;
1706 
1707  if(m_linSysKey.GetMatrixType() == StdRegions::eMass)
1708  {
1709  PreconR = StdRegions::ePreconRMass;
1710  PreconRT = StdRegions::ePreconRTMass;
1711  }
1712  else
1713  {
1714  PreconR = StdRegions::ePreconR;
1715  PreconRT = StdRegions::ePreconRT;
1716  }
1717 
1718 
1719 
1720  /*
1721  * Matrix keys - for each element type there are two matrix keys
1722  * corresponding to the transformation matrix R and its transpose
1723  */
1724 
1725 
1726  //Matrix keys for tetrahedral element transformation matrix
1728  (PreconR, LibUtilities::eTetrahedron,
1729  *TetExp, m_linSysKey.GetConstFactors(),
1730  vVarCoeffMap);
1731 
1732  //Matrix keys for tetrahedral transposed transformation matrix
1734  (PreconRT, LibUtilities::eTetrahedron,
1735  *TetExp, m_linSysKey.GetConstFactors(),
1736  vVarCoeffMap);
1737 
1738  //Matrix keys for prismatic element transformation matrix
1740  (PreconR, LibUtilities::ePrism,
1741  *PrismExp, m_linSysKey.GetConstFactors(),
1742  vVarCoeffMap);
1743 
1744  //Matrix keys for prismatic element transposed transformation matrix
1745  LocalRegions::MatrixKey PrismRT
1746  (PreconRT, LibUtilities::ePrism,
1747  *PrismExp, m_linSysKey.GetConstFactors(),
1748  vVarCoeffMap);
1749 
1750  //Matrix keys for hexahedral element transformation matrix
1752  (PreconR, LibUtilities::eHexahedron,
1753  *HexExp, m_linSysKey.GetConstFactors(),
1754  vVarCoeffMap);
1755 
1756  //Matrix keys for hexahedral element transposed transformation
1757  //matrix
1759  (PreconRT, LibUtilities::eHexahedron,
1760  *HexExp, m_linSysKey.GetConstFactors(),
1761  vVarCoeffMap);
1762 
1763  /*
1764  * Create transformation matrices for the tetrahedral element
1765  */
1766 
1767  //Get tetrahedral transformation matrix
1768  m_Rtet = TetExp->GetLocMatrix(TetR);
1769 
1770  //Get tetrahedral transposed transformation matrix
1771  m_RTtet = TetExp->GetLocMatrix(TetRT);
1772 
1773  // Using the transformation matrix and the inverse transformation
1774  // matrix create the inverse matrices
1775  Rtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1776 
1777  //Inverse transposed transformation matrix
1778  RTtettmp=TetExp->BuildInverseTransformationMatrix(m_Rtet);
1779  RTtettmp->Transpose();
1780 
1782  ::AllocateSharedPtr(1.0,Rtettmp);
1784  ::AllocateSharedPtr(1.0,RTtettmp);
1785 
1786  /*
1787  * Create transformation matrices for the hexahedral element
1788  */
1789 
1790  //Get hexahedral transformation matrix
1791  m_Rhex = HexExp->GetLocMatrix(HexR);
1792  //Get hexahedral transposed transformation matrix
1793  m_RThex = HexExp->GetLocMatrix(HexRT);
1794 
1795  // Using the transformation matrix and the inverse transformation
1796  // matrix create the inverse matrices
1797  Rhextmp=HexExp->BuildInverseTransformationMatrix(m_Rhex);
1798  //Inverse transposed transformation matrix
1799  RThextmp=HexExp->BuildInverseTransformationMatrix(m_Rhex);
1800  RThextmp->Transpose();
1801 
1803  ::AllocateSharedPtr(1.0,Rhextmp);
1805  ::AllocateSharedPtr(1.0,RThextmp);
1806 
1807  /*
1808  * Create transformation matrices for the prismatic element
1809  */
1810 
1811  //Get prism transformation matrix
1812  Rprismoriginal = PrismExp->GetLocMatrix(PrismR);
1813  //Get prism transposed transformation matrix
1814  RTprismoriginal = PrismExp->GetLocMatrix(PrismRT);
1815 
1816  unsigned int nRows=Rprismoriginal->GetRows();
1817  NekDouble zero=0.0;
1819  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1821  AllocateSharedPtr(nRows,nRows,zero,eFULL);
1822  NekDouble Rvalue, RTvalue;
1823 
1824  //Copy values from the prism transformation matrix
1825  for(i=0; i<nRows; ++i)
1826  {
1827  for(j=0; j<nRows; ++j)
1828  {
1829  Rvalue=(*Rprismoriginal)(i,j);
1830  RTvalue=(*RTprismoriginal)(i,j);
1831  Rtmpprism->SetValue(i,j,Rvalue);
1832  RTtmpprism->SetValue(i,j,RTvalue);
1833  }
1834  }
1835 
1836  //Replace triangular faces and edges of the prims transformation
1837  //matrix with the corresponding values of the tetrahedral
1838  //transformation matrix.
1839  ModifyPrismTransformationMatrix(TetExp,PrismExp,Rtmpprism,RTtmpprism);
1840 
1842  ::AllocateSharedPtr(1.0,Rtmpprism);
1843 
1845  ::AllocateSharedPtr(1.0,RTtmpprism);
1846 
1847  //Inverse transformation matrix
1848  Rprismtmp=PrismExp->BuildInverseTransformationMatrix(m_Rprism);
1849 
1850  //Inverse transposed transformation matrix
1851  RTprismtmp=PrismExp->BuildInverseTransformationMatrix(m_Rprism);
1852  RTprismtmp->Transpose();
1853 
1855  ::AllocateSharedPtr(1.0,Rprismtmp);
1856 
1858  ::AllocateSharedPtr(1.0,RTprismtmp);
1859  }
1860 
1861  /**
1862  * \brief Modify the prism transformation matrix to align with the
1863  * tetrahedral modes.
1864  *
1865  * This routine replaces the edge and triangular face components of the
1866  * prismatic vertex transformation matrices \f$\mathbf{R}_{ve}\f$ and
1867  * \f$\mathbf{R}_{vf}\f$ with the corresponding components from the
1868  * tetrahedral transformation matrices. Additionally, triangular face
1869  * components in the prismatic edge transformation matrix
1870  * \f$\mathbf{R}_{ef}\f$ with the corresponding component from the
1871  * tetrahedral transformation matrix.
1872  */
1876  DNekMatSharedPtr Rmodprism,
1877  DNekMatSharedPtr RTmodprism)
1878  {
1879  NekDouble Rvalue, RTvalue;
1880  int i, j;
1881 
1882  //For a tet element the bottom face is made up of the following:
1883  //vertices: 0, 1 and 2 edges: 0, 1 and 2 face: 0. We first need to
1884  //determine the mode locations of these vertices, edges and face so
1885  //we can extract the correct values from the tetrahedral R matrix.
1886 
1887  //These are the vertex mode locations of R which need to be replaced
1888  //in the prism element
1889  int TetVertex0=TetExp->GetVertexMap(0);
1890  int TetVertex1=TetExp->GetVertexMap(1);
1891  int TetVertex2=TetExp->GetVertexMap(2);
1892  int TetVertex3=TetExp->GetVertexMap(3);
1893 
1894 
1895  //These are the edge mode locations of R which need to be replaced
1896  //in the prism element
1897  Array<OneD, unsigned int> TetEdge0=TetExp->GetEdgeInverseBoundaryMap(0);
1898  Array<OneD, unsigned int> TetEdge1=TetExp->GetEdgeInverseBoundaryMap(1);
1899  Array<OneD, unsigned int> TetEdge2=TetExp->GetEdgeInverseBoundaryMap(2);
1900  Array<OneD, unsigned int> TetEdge3=TetExp->GetEdgeInverseBoundaryMap(3);
1901  Array<OneD, unsigned int> TetEdge4=TetExp->GetEdgeInverseBoundaryMap(4);
1902  Array<OneD, unsigned int> TetEdge5=TetExp->GetEdgeInverseBoundaryMap(5);
1903 
1904  //These are the face mode locations of R which need to be replaced
1905  //in the prism element
1906  Array<OneD, unsigned int> TetFace=TetExp->GetFaceInverseBoundaryMap(1);
1907 
1908  //Prism vertex modes
1909  int PrismVertex0=PrismExp->GetVertexMap(0);
1910  int PrismVertex1=PrismExp->GetVertexMap(1);
1911  int PrismVertex2=PrismExp->GetVertexMap(2);
1912  int PrismVertex3=PrismExp->GetVertexMap(3);
1913  int PrismVertex4=PrismExp->GetVertexMap(4);
1914  int PrismVertex5=PrismExp->GetVertexMap(5);
1915 
1916  //Prism edge modes
1917  Array<OneD, unsigned int> PrismEdge0=
1918  PrismExp->GetEdgeInverseBoundaryMap(0);
1919  Array<OneD, unsigned int> PrismEdge1=
1920  PrismExp->GetEdgeInverseBoundaryMap(1);
1921  Array<OneD, unsigned int> PrismEdge2=
1922  PrismExp->GetEdgeInverseBoundaryMap(2);
1923  Array<OneD, unsigned int> PrismEdge3=
1924  PrismExp->GetEdgeInverseBoundaryMap(3);
1925  Array<OneD, unsigned int> PrismEdge4=
1926  PrismExp->GetEdgeInverseBoundaryMap(4);
1927  Array<OneD, unsigned int> PrismEdge5=
1928  PrismExp->GetEdgeInverseBoundaryMap(5);
1929  Array<OneD, unsigned int> PrismEdge6=
1930  PrismExp->GetEdgeInverseBoundaryMap(6);
1931  Array<OneD, unsigned int> PrismEdge7=
1932  PrismExp->GetEdgeInverseBoundaryMap(7);
1933  Array<OneD, unsigned int> PrismEdge8=
1934  PrismExp->GetEdgeInverseBoundaryMap(8);
1935 
1936  //Prism face 1 & 3 face modes
1937  Array<OneD, unsigned int> PrismFace1=
1938  PrismExp->GetFaceInverseBoundaryMap(1);
1939  Array<OneD, unsigned int> PrismFace3=
1940  PrismExp->GetFaceInverseBoundaryMap(3);
1941  Array<OneD, unsigned int> PrismFace0=
1942  PrismExp->GetFaceInverseBoundaryMap(0);
1943  Array<OneD, unsigned int> PrismFace2=
1944  PrismExp->GetFaceInverseBoundaryMap(2);
1945  Array<OneD, unsigned int> PrismFace4=
1946  PrismExp->GetFaceInverseBoundaryMap(4);
1947 
1948  //vertex 0 edge 0 3 & 4
1949  for(i=0; i< PrismEdge0.num_elements(); ++i)
1950  {
1951  Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
1952  Rmodprism->SetValue(PrismVertex0,PrismEdge0[i],Rvalue);
1953  Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
1954  Rmodprism->SetValue(PrismVertex0,PrismEdge3[i],Rvalue);
1955  Rvalue=(*m_Rtet)(TetVertex0,TetEdge3[i]);
1956  Rmodprism->SetValue(PrismVertex0,PrismEdge4[i],Rvalue);
1957 
1958  //transposed values
1959  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
1960  RTmodprism->SetValue(PrismEdge0[i],PrismVertex0,RTvalue);
1961  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
1962  RTmodprism->SetValue(PrismEdge3[i],PrismVertex0,RTvalue);
1963  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex0);
1964  RTmodprism->SetValue(PrismEdge4[i],PrismVertex0,RTvalue);
1965  }
1966 
1967  //vertex 1 edge 0 1 & 5
1968  for(i=0; i< PrismEdge1.num_elements(); ++i)
1969  {
1970  Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1971  Rmodprism->SetValue(PrismVertex1,PrismEdge0[i],Rvalue);
1972  Rvalue=(*m_Rtet)(TetVertex1,TetEdge1[i]);
1973  Rmodprism->SetValue(PrismVertex1,PrismEdge1[i],Rvalue);
1974  Rvalue=(*m_Rtet)(TetVertex1,TetEdge4[i]);
1975  Rmodprism->SetValue(PrismVertex1,PrismEdge5[i],Rvalue);
1976 
1977  //transposed values
1978  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
1979  RTmodprism->SetValue(PrismEdge0[i],PrismVertex1,RTvalue);
1980  RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex1);
1981  RTmodprism->SetValue(PrismEdge1[i],PrismVertex1,RTvalue);
1982  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex1);
1983  RTmodprism->SetValue(PrismEdge5[i],PrismVertex1,RTvalue);
1984  }
1985 
1986  //vertex 2 edge 1 2 & 6
1987  for(i=0; i< PrismEdge2.num_elements(); ++i)
1988  {
1989  Rvalue=(*m_Rtet)(TetVertex2,TetEdge1[i]);
1990  Rmodprism->SetValue(PrismVertex2,PrismEdge1[i],Rvalue);
1991  Rvalue=(*m_Rtet)(TetVertex1,TetEdge0[i]);
1992  Rmodprism->SetValue(PrismVertex2,PrismEdge2[i],Rvalue);
1993  Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
1994  Rmodprism->SetValue(PrismVertex2,PrismEdge6[i],Rvalue);
1995 
1996  //transposed values
1997  RTvalue=(*m_RTtet)(TetEdge1[i],TetVertex2);
1998  RTmodprism->SetValue(PrismEdge1[i],PrismVertex2,RTvalue);
1999  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex1);
2000  RTmodprism->SetValue(PrismEdge2[i],PrismVertex2,RTvalue);
2001  RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2002  RTmodprism->SetValue(PrismEdge6[i],PrismVertex2,RTvalue);
2003  }
2004 
2005  //vertex 3 edge 3 2 & 7
2006  for(i=0; i< PrismEdge3.num_elements(); ++i)
2007  {
2008  Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2009  Rmodprism->SetValue(PrismVertex3,PrismEdge3[i],Rvalue);
2010  Rvalue=(*m_Rtet)(TetVertex0,TetEdge0[i]);
2011  Rmodprism->SetValue(PrismVertex3,PrismEdge2[i],Rvalue);
2012  Rvalue=(*m_Rtet)(TetVertex2,TetEdge5[i]);
2013  Rmodprism->SetValue(PrismVertex3,PrismEdge7[i],Rvalue);
2014 
2015  //transposed values
2016  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2017  RTmodprism->SetValue(PrismEdge3[i],PrismVertex3,RTvalue);
2018  RTvalue=(*m_RTtet)(TetEdge0[i],TetVertex0);
2019  RTmodprism->SetValue(PrismEdge2[i],PrismVertex3,RTvalue);
2020  RTvalue=(*m_RTtet)(TetEdge5[i],TetVertex2);
2021  RTmodprism->SetValue(PrismEdge7[i],PrismVertex3,RTvalue);
2022  }
2023 
2024  //vertex 4 edge 4 5 & 8
2025  for(i=0; i< PrismEdge4.num_elements(); ++i)
2026  {
2027  Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2028  Rmodprism->SetValue(PrismVertex4,PrismEdge4[i],Rvalue);
2029  Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2030  Rmodprism->SetValue(PrismVertex4,PrismEdge5[i],Rvalue);
2031  Rvalue=(*m_Rtet)(TetVertex0,TetEdge2[i]);
2032  Rmodprism->SetValue(PrismVertex4,PrismEdge8[i],Rvalue);
2033 
2034  //transposed values
2035  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2036  RTmodprism->SetValue(PrismEdge4[i],PrismVertex4,RTvalue);
2037  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2038  RTmodprism->SetValue(PrismEdge5[i],PrismVertex4,RTvalue);
2039  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex0);
2040  RTmodprism->SetValue(PrismEdge8[i],PrismVertex4,RTvalue);
2041  }
2042 
2043  //vertex 5 edge 6 7 & 8
2044  for(i=0; i< PrismEdge5.num_elements(); ++i)
2045  {
2046  Rvalue=(*m_Rtet)(TetVertex3,TetEdge3[i]);
2047  Rmodprism->SetValue(PrismVertex5,PrismEdge6[i],Rvalue);
2048  Rvalue=(*m_Rtet)(TetVertex3,TetEdge4[i]);
2049  Rmodprism->SetValue(PrismVertex5,PrismEdge7[i],Rvalue);
2050  Rvalue=(*m_Rtet)(TetVertex2,TetEdge2[i]);
2051  Rmodprism->SetValue(PrismVertex5,PrismEdge8[i],Rvalue);
2052 
2053  //transposed values
2054  RTvalue=(*m_RTtet)(TetEdge3[i],TetVertex3);
2055  RTmodprism->SetValue(PrismEdge6[i],PrismVertex5,RTvalue);
2056  RTvalue=(*m_RTtet)(TetEdge4[i],TetVertex3);
2057  RTmodprism->SetValue(PrismEdge7[i],PrismVertex5,RTvalue);
2058  RTvalue=(*m_RTtet)(TetEdge2[i],TetVertex2);
2059  RTmodprism->SetValue(PrismEdge8[i],PrismVertex5,RTvalue);
2060  }
2061 
2062  // face 1 vertices 0 1 4
2063  for(i=0; i< PrismFace1.num_elements(); ++i)
2064  {
2065  Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2066  Rmodprism->SetValue(PrismVertex0,PrismFace1[i],Rvalue);
2067  Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2068  Rmodprism->SetValue(PrismVertex1,PrismFace1[i],Rvalue);
2069  Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2070  Rmodprism->SetValue(PrismVertex4,PrismFace1[i],Rvalue);
2071 
2072  //transposed values
2073  RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2074  RTmodprism->SetValue(PrismFace1[i],PrismVertex0,RTvalue);
2075  RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2076  RTmodprism->SetValue(PrismFace1[i],PrismVertex1,RTvalue);
2077  RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2078  RTmodprism->SetValue(PrismFace1[i],PrismVertex4,RTvalue);
2079  }
2080 
2081  // face 3 vertices 2, 3 & 5
2082  for(i=0; i< PrismFace3.num_elements(); ++i)
2083  {
2084  Rvalue=(*m_Rtet)(TetVertex1,TetFace[i]);
2085  Rmodprism->SetValue(PrismVertex2,PrismFace3[i],Rvalue);
2086  Rvalue=(*m_Rtet)(TetVertex0,TetFace[i]);
2087  Rmodprism->SetValue(PrismVertex3,PrismFace3[i],Rvalue);
2088  Rvalue=(*m_Rtet)(TetVertex3,TetFace[i]);
2089  Rmodprism->SetValue(PrismVertex5,PrismFace3[i],Rvalue);
2090 
2091  //transposed values
2092  RTvalue=(*m_RTtet)(TetFace[i],TetVertex1);
2093  RTmodprism->SetValue(PrismFace3[i],PrismVertex2,RTvalue);
2094  RTvalue=(*m_RTtet)(TetFace[i],TetVertex0);
2095  RTmodprism->SetValue(PrismFace3[i],PrismVertex3,RTvalue);
2096  RTvalue=(*m_RTtet)(TetFace[i],TetVertex3);
2097  RTmodprism->SetValue(PrismFace3[i],PrismVertex5,RTvalue);
2098  }
2099 
2100  // Face 1 edge 0 4 5
2101  for(i=0; i< PrismFace1.num_elements(); ++i)
2102  {
2103  for(j=0; j<PrismEdge0.num_elements(); ++j)
2104  {
2105  Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2106  Rmodprism->SetValue(PrismEdge0[j],PrismFace1[i],Rvalue);
2107  Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2108  Rmodprism->SetValue(PrismEdge4[j],PrismFace1[i],Rvalue);
2109  Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2110  Rmodprism->SetValue(PrismEdge5[j],PrismFace1[i],Rvalue);
2111 
2112  //transposed values
2113  RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2114  RTmodprism->SetValue(PrismFace1[i],PrismEdge0[j],RTvalue);
2115  RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2116  RTmodprism->SetValue(PrismFace1[i],PrismEdge4[j],RTvalue);
2117  RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2118  RTmodprism->SetValue(PrismFace1[i],PrismEdge5[j],RTvalue);
2119  }
2120  }
2121 
2122  // Face 3 edge 2 6 7
2123  for(i=0; i< PrismFace3.num_elements(); ++i)
2124  {
2125  for(j=0; j<PrismEdge2.num_elements(); ++j)
2126  {
2127  Rvalue=(*m_Rtet)(TetEdge0[j],TetFace[i]);
2128  Rmodprism->SetValue(PrismEdge2[j],PrismFace3[i],Rvalue);
2129  Rvalue=(*m_Rtet)(TetEdge4[j],TetFace[i]);
2130  Rmodprism->SetValue(PrismEdge6[j],PrismFace3[i],Rvalue);
2131  Rvalue=(*m_Rtet)(TetEdge3[j],TetFace[i]);
2132  Rmodprism->SetValue(PrismEdge7[j],PrismFace3[i],Rvalue);
2133 
2134  RTvalue=(*m_RTtet)(TetFace[i],TetEdge0[j]);
2135  RTmodprism->SetValue(PrismFace3[i],PrismEdge2[j],RTvalue);
2136  RTvalue=(*m_RTtet)(TetFace[i],TetEdge4[j]);
2137  RTmodprism->SetValue(PrismFace3[i],PrismEdge6[j],RTvalue);
2138  RTvalue=(*m_RTtet)(TetFace[i],TetEdge3[j]);
2139  RTmodprism->SetValue(PrismFace3[i],PrismEdge7[j],RTvalue);
2140  }
2141  }
2142  }
2143 
2144  }
2145 }
2146 
2147 
2148 
2149 
2150 
2151 
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:198
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:239
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:644
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
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:271
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:228
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
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:227
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:709
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< PrismExp > PrismExpSharedPtr
Definition: PrismExp.h:222
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.
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
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:60
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:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:52
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:183
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215