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