Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PreconditionerBlock.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PreconditionerBlock.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Block Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/BasicUtils/VDmathArray.hpp>
37 #include <MultiRegions/PreconditionerBlock.h>
38 #include <MultiRegions/GlobalMatrixKey.h>
39 #include <MultiRegions/GlobalLinSysIterativeStaticCond.h>
40 #include <MultiRegions/GlobalLinSys.h>
41 #include <LocalRegions/MatrixKey.h>
42 #include <LocalRegions/SegExp.h>
43 #include <math.h>
44 
45 namespace Nektar
46 {
47  namespace MultiRegions
48  {
49  /**
50  * Registers the class with the Factory.
51  */
52  string PreconditionerBlock::className
53  = GetPreconFactory().RegisterCreatorFunction(
54  "Block",
55  PreconditionerBlock::create,
56  "Block Diagonal Preconditioning");
57  /**
58  * @class Block Preconditioner
59  *
60  * This class implements Block preconditioning for the conjugate
61  * gradient matrix solver.
62  */
63 
64  PreconditionerBlock::PreconditionerBlock(
65  const boost::shared_ptr<GlobalLinSys> &plinsys,
66  const AssemblyMapSharedPtr &pLocToGloMap)
67  : Preconditioner(plinsys, pLocToGloMap),
68  m_linsys(plinsys),
69  m_preconType(pLocToGloMap->GetPreconType()),
70  m_locToGloMap(pLocToGloMap)
71  {
72  }
73 
74  void PreconditionerBlock::v_InitObject()
75  {
76  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
77  ASSERTL0(solvertype == MultiRegions::eIterativeStaticCond,"Solver type not valid");
78  }
79 
80 
81  void PreconditionerBlock::v_BuildPreconditioner()
82  {
83 
84  boost::shared_ptr<MultiRegions::ExpList>
85  expList=((m_linsys.lock())->GetLocMat()).lock();
86  StdRegions::StdExpansionSharedPtr locExpansion;
87  locExpansion = expList->GetExp(0);
88  int nDim = locExpansion->GetShapeDimension();
89 
90  if (nDim == 1)
91  {
92  ASSERTL0(0,"Unknown preconditioner");
93  }
94  else if (nDim == 2)
95  {
96  BlockPreconditioner2D();
97  }
98  else if (nDim ==3)
99  {
100  BlockPreconditioner3D();
101  }
102 
103  }
104 
105  /**
106  * \brief Construct a block preconditioner from
107  * \f$\mathbf{S}_{1}\f$
108  *
109  *\f[\mathbf{M}^{-1}=\left[\begin{array}{ccc}
110  * Diag[(\mathbf{S_{1}})_{vv}] & & \\ & (\mathbf{S}_{1})_{eb} & \\ & &
111  * (\mathbf{S}_{1})_{fb} \end{array}\right] \f]
112  *
113  * where \f$\mathbf{S}_{1}\f$ is the local schur complement matrix for
114  * each element.
115  */
116  void PreconditionerBlock::BlockPreconditioner2D()
117  {
118  boost::shared_ptr<MultiRegions::ExpList>
119  expList=((m_linsys.lock())->GetLocMat()).lock();
120  LocalRegions::ExpansionSharedPtr locExpansion;
121  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
122  StdRegions::VarCoeffMap vVarCoeffMap;
123  int i, j, k, nel;
124  int nVerts, nEdges;
125  int eid, n, cnt, nedgemodes;
126  NekDouble zero = 0.0;
127 
128  int vMap1, vMap2, sign1, sign2;
129  int m, v, eMap1, eMap2;
130  int offset, globalrow, globalcol;
131 
132  //matrix storage
133  MatrixStorage storage = eFULL;
134  MatrixStorage vertstorage = eDIAGONAL;
135  MatrixStorage blkmatStorage = eDIAGONAL;
136 
137  //local element static condensed matrices
138  DNekScalBlkMatSharedPtr loc_mat;
139  DNekScalMatSharedPtr bnd_mat;
140 
141  DNekMatSharedPtr VertBlk;
142 
143  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
144  int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
145 
146  //Vertex and edge preconditioner matrices
147  VertBlk = MemoryManager<DNekMat>::
148  AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
149 
150  Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
151  Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
152 
153  int n_exp = expList->GetNumElmts();
154  int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
155 
156  //set the number of blocks in the matrix
157  Array<OneD,unsigned int> n_blks(1+nNonDirEdgeIDs);
158  n_blks[0]=nNonDirVerts;
159 
160  map<int,int> edgeDirMap;
161  map<int,int> uniqueEdgeMap;
162 
163  //this should be of size total number of local edges
164  Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
165  Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
166 
167  const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
168  StdRegions::StdExpansion1DSharedPtr bndCondFaceExp;
169  LocalRegions::SegExpSharedPtr bndSegExp;
170  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>&
171  bndConditions = expList->GetBndConditions();
172 
173  int meshVertId;
174  int meshEdgeId;
175 
176  //Determine which boundary edges and faces have dirichlet values
177  for(i = 0; i < bndCondExp.num_elements(); i++)
178  {
179  cnt = 0;
180  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
181  {
182  bndSegExp = bndCondExp[i]->GetExp(j)
183  ->as<LocalRegions::SegExp>();
184  if (bndConditions[i]->GetBoundaryConditionType() ==
185  SpatialDomains::eDirichlet)
186  {
187  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
188  edgeDirMap[meshEdgeId] = 1;
189  }
190  }
191  }
192 
193  int dof=0;
194  int maxEdgeDof=0;
195  int nlocalNonDirEdges=0;
196 
197  int edgematrixlocation=0;
198  int ntotaledgeentries=0;
199 
200  // Loop over all the elements in the domain and compute max edge
201  // DOF. Reduce across all processes to get universal maximum.
202  for(cnt=n=0; n < n_exp; ++n)
203  {
204  nel = expList->GetOffset_Elmt_Id(n);
205 
206  locExpansion = expList->GetExp(nel);
207 
208  for (j = 0; j < locExpansion->GetNedges(); ++j)
209  {
210  dof = locExpansion->GetEdgeNcoeffs(j)-2;
211  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
212  meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()
213  ->GetGeom2D()->GetEid(j);
214 
215  if(edgeDirMap.count(meshEdgeId)==0)
216  {
217  if(uniqueEdgeMap.count(meshEdgeId)==0)
218  {
219  uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
220 
221  edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
222 
223  edgemodeoffset[edgematrixlocation]=dof*dof;
224 
225  ntotaledgeentries+=dof*dof;
226 
227  n_blks[1+edgematrixlocation++]=dof;
228 
229  }
230 
231  nlocalNonDirEdges+=dof*dof;
232  }
233  }
234  }
235 
236  m_comm = expList->GetComm()->GetRowComm();
237  m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
238 
239  //Allocate arrays for block to universal map (number of expansions * p^2)
240  Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
241 
242  Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
243 
244  //Allocate arrays to store matrices (number of expansions * p^2)
245  Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
246 
247  int edgematrixoffset=0;
248  int vGlobal;
249 
250  for(n=0; n < n_exp; ++n)
251  {
252  nel = expList->GetOffset_Elmt_Id(n);
253 
254  locExpansion = expList->GetExp(nel);
255 
256  //loop over the edges of the expansion
257  for(j = 0; j < locExpansion->GetNedges(); ++j)
258  {
259  //get mesh edge id
260  meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()
261  ->GetGeom2D()->GetEid(j);
262 
263  nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
264 
265  if(edgeDirMap.count(meshEdgeId)==0)
266  {
267  for(k=0; k<nedgemodes*nedgemodes; ++k)
268  {
269  vGlobal=edgeglobaloffset[uniqueEdgeMap[meshEdgeId]]+k;
270 
271 
272  localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
273 
274  EdgeBlockToUniversalMap[vGlobal]
275  = meshEdgeId * maxEdgeDof * maxEdgeDof + k + 1;
276  }
277  edgematrixoffset+=nedgemodes*nedgemodes;
278  }
279  }
280  }
281 
282  edgematrixoffset=0;
283 
284  m_blkMat = MemoryManager<DNekBlkMat>
285  ::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
286 
287  //Here we loop over the expansion and build the block low energy
288  //preconditioner as well as the block versions of the transformation
289  //matrices.
290  for(cnt=n=0; n < n_exp; ++n)
291  {
292  nel = expList->GetOffset_Elmt_Id(n);
293 
294  locExpansion = expList->GetExp(nel);
295 
296  nVerts=locExpansion->GetGeom()->GetNumVerts();
297  nEdges=locExpansion->GetGeom()->GetNumEdges();
298 
299  //Get statically condensed matrix
300  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
301 
302  //Extract boundary block (elemental S1)
303  bnd_mat=loc_mat->GetBlock(0,0);
304 
305  //offset by number of rows
306  offset = bnd_mat->GetRows();
307 
308  DNekScalMat &S=(*bnd_mat);
309 
310  //loop over vertices of the element and return the vertex map
311  //for each vertex
312  for (v=0; v<nVerts; ++v)
313  {
314  vMap1=locExpansion->GetVertexMap(v);
315 
316  //Get vertex map
317  globalrow = m_locToGloMap->
318  GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
319 
320  if(globalrow >= 0)
321  {
322  for (m=0; m<nVerts; ++m)
323  {
324  vMap2=locExpansion->GetVertexMap(m);
325 
326  //global matrix location (without offset due to
327  //dirichlet values)
328  globalcol = m_locToGloMap->
329  GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
330 
331  //offset for dirichlet conditions
332  if (globalcol == globalrow)
333  {
334  meshVertId = locExpansion->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetVid(v);
335 
336  //modal connectivity between elements
337  sign1 = m_locToGloMap->
338  GetLocalToGlobalBndSign(cnt + vMap1);
339  sign2 = m_locToGloMap->
340  GetLocalToGlobalBndSign(cnt + vMap2);
341 
342  vertArray[globalrow]
343  += sign1*sign2*S(vMap1,vMap2);
344 
345  VertBlockToUniversalMap[globalrow]
346  = meshVertId * maxEdgeDof * maxEdgeDof + 1;
347  }
348  }
349  }
350  }
351 
352  //loop over edges of the element and return the edge map
353  for (eid=0; eid<nEdges; ++eid)
354  {
355  nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
356 
357  DNekMatSharedPtr locMat =
358  MemoryManager<DNekMat>::AllocateSharedPtr
359  (nedgemodes,nedgemodes,zero,storage);
360 
361  meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(eid);
362  Array<OneD, unsigned int> edgemodearray =
363  locExpansion->GetEdgeInverseBoundaryMap(eid);
364 
365  if(edgeDirMap.count(meshEdgeId)==0)
366  {
367  for (v=0; v<nedgemodes; ++v)
368  {
369  eMap1=edgemodearray[v];
370 
371  for (m=0; m<nedgemodes; ++m)
372  {
373  eMap2=edgemodearray[m];
374 
375  //modal connectivity between elements
376  sign1 = m_locToGloMap->
377  GetLocalToGlobalBndSign(cnt + eMap1);
378  sign2 = m_locToGloMap->
379  GetLocalToGlobalBndSign(cnt + eMap2);
380 
381  NekDouble globalEdgeValue =
382  sign1*sign2*S(eMap1,eMap2);
383 
384  EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=
385  globalEdgeValue;
386  }
387  }
388  edgematrixoffset+=nedgemodes*nedgemodes;
389  }
390  }
391 
392  //offset for the expansion
393  cnt+=offset;
394  }
395 
396  //Assemble edge matrices of each process
397  Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries);
398  Vmath::Zero(ntotaledgeentries, GlobalEdgeBlock.get(), 1);
399  Vmath::Assmb(EdgeBlockArray.num_elements(),
400  EdgeBlockArray.get(),
401  localEdgeToGlobalMatrixMap.get(),
402  GlobalEdgeBlock.get());
403 
404  //Exchange vertex data over different processes
405  if(nNonDirVerts!=0)
406  {
407  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
408  Gs::Gather(vertArray, Gs::gs_add, tmp);
409  }
410 
411  //Exchange edge data over different processes
412  Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
413  Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
414 
415  // Populate vertex block
416  for (int i = 0; i < nNonDirVerts; ++i)
417  {
418  VertBlk->SetValue(i,i,1.0/vertArray[i]);
419  }
420 
421  //Set the first block to be the diagonal of the vertex space
422  m_blkMat->SetBlock(0,0, VertBlk);
423 
424  offset=0;
425  //Build the edge matrices from the vector
426  for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
427  {
428  DNekMatSharedPtr gmat =
429  MemoryManager<DNekMat>::AllocateSharedPtr
430  (nedgemodes,nedgemodes,zero,storage);
431 
432  for (v=0; v<nedgemodes; ++v)
433  {
434  for (m=0; m<nedgemodes; ++m)
435  {
436  NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
437  gmat->SetValue(v,m,EdgeValue);
438  }
439  }
440 
441  m_blkMat->SetBlock(1+loc,1+loc, gmat);
442 
443  offset+=edgemodeoffset[loc];
444  }
445 
446  int totblks=m_blkMat->GetNumberOfBlockRows();
447  for (i=1; i< totblks; ++i)
448  {
449  unsigned int nmodes=m_blkMat->GetNumberOfRowsInBlockRow(i);
450  DNekMatSharedPtr tmp_mat =
451  MemoryManager<DNekMat>::AllocateSharedPtr
452  (nmodes,nmodes,zero,storage);
453 
454  tmp_mat=m_blkMat->GetBlock(i,i);
455  tmp_mat->Invert();
456  m_blkMat->SetBlock(i,i,tmp_mat);
457  }
458  }
459 
460  /**
461  *
462  */
463  void PreconditionerBlock::BlockPreconditioner3D()
464  {
465  boost::shared_ptr<MultiRegions::ExpList>
466  expList=((m_linsys.lock())->GetLocMat()).lock();
467  LocalRegions::ExpansionSharedPtr locExpansion;
468  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
469  StdRegions::VarCoeffMap vVarCoeffMap;
470  int i, j, k, nel;
471  int nVerts, nEdges,nFaces;
472  int eid, fid, n, cnt, nedgemodes, nfacemodes;
473  NekDouble zero = 0.0;
474 
475  int vMap1, vMap2, sign1, sign2;
476  int m, v, eMap1, eMap2, fMap1, fMap2;
477  int offset, globalrow, globalcol;
478 
479  //matrix storage
480  MatrixStorage storage = eFULL;
481  MatrixStorage vertstorage = eDIAGONAL;
482  MatrixStorage blkmatStorage = eDIAGONAL;
483 
484  //local element static condensed matrices
485  DNekScalBlkMatSharedPtr loc_mat;
486  DNekScalMatSharedPtr bnd_mat;
487 
488  DNekMatSharedPtr VertBlk;
489 
490  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
491  int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
492 
493  //Vertex, edge and face preconditioner matrices
494  VertBlk = MemoryManager<DNekMat>::
495  AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
496 
497  Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
498  Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
499 
500  int n_exp = expList->GetNumElmts();
501  int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
502  int nNonDirFaceIDs=m_locToGloMap->GetNumNonDirFaces();
503 
504  //set the number of blocks in the matrix
505  Array<OneD,unsigned int> n_blks(1+nNonDirEdgeIDs+nNonDirFaceIDs);
506  n_blks[0]=nNonDirVerts;
507 
508  map<int,int> edgeDirMap;
509  map<int,int> faceDirMap;
510  map<int,int> uniqueEdgeMap;
511  map<int,int> uniqueFaceMap;
512 
513  //this should be of size total number of local edges
514  Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
515  Array<OneD, int> facemodeoffset(nNonDirFaceIDs,0);
516 
517  Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
518  Array<OneD, int> faceglobaloffset(nNonDirFaceIDs,0);
519 
520  const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
521  StdRegions::StdExpansion2DSharedPtr bndCondFaceExp;
522  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>& bndConditions = expList->GetBndConditions();
523 
524  int meshVertId;
525  int meshEdgeId;
526  int meshFaceId;
527 
528  const Array<OneD, const int> &extradiredges
529  = m_locToGloMap->GetExtraDirEdges();
530  for(i=0; i<extradiredges.num_elements(); ++i)
531  {
532  meshEdgeId=extradiredges[i];
533  edgeDirMap[meshEdgeId] = 1;
534  }
535 
536  //Determine which boundary edges and faces have dirichlet values
537  for(i = 0; i < bndCondExp.num_elements(); i++)
538  {
539  cnt = 0;
540  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
541  {
542  bndCondFaceExp = bndCondExp[i]->GetExp(j)->as<StdRegions::StdExpansion2D>();
543  if (bndConditions[i]->GetBoundaryConditionType() ==
544  SpatialDomains::eDirichlet)
545  {
546  for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
547  {
548  meshEdgeId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(k);
549  if(edgeDirMap.count(meshEdgeId) == 0)
550  {
551  edgeDirMap[meshEdgeId] = 1;
552  }
553  }
554  meshFaceId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetFid();
555  faceDirMap[meshFaceId] = 1;
556  }
557  }
558  }
559 
560  int dof=0;
561  int maxFaceDof=0;
562  int maxEdgeDof=0;
563  int nlocalNonDirEdges=0;
564  int nlocalNonDirFaces=0;
565 
566  int edgematrixlocation=0;
567  int ntotaledgeentries=0;
568 
569  // Loop over all the elements in the domain and compute max edge
570  // DOF. Reduce across all processes to get universal maximum.
571  for(cnt=n=0; n < n_exp; ++n)
572  {
573  nel = expList->GetOffset_Elmt_Id(n);
574 
575  locExpansion = expList->GetExp(nel);
576 
577  for (j = 0; j < locExpansion->GetNedges(); ++j)
578  {
579  dof = locExpansion->GetEdgeNcoeffs(j)-2;
580  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
581  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
582 
583  if(edgeDirMap.count(meshEdgeId)==0)
584  {
585  if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
586  {
587  uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
588 
589  edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
590 
591  edgemodeoffset[edgematrixlocation]=dof*dof;
592 
593  ntotaledgeentries+=dof*dof;
594 
595  n_blks[1+edgematrixlocation++]=dof;
596  }
597 
598  nlocalNonDirEdges+=dof*dof;
599  }
600  }
601  }
602 
603  int facematrixlocation=0;
604  int ntotalfaceentries=0;
605 
606  // Loop over all the elements in the domain and compute max face
607  // DOF. Reduce across all processes to get universal maximum.
608  for(cnt=n=0; n < n_exp; ++n)
609  {
610  nel = expList->GetOffset_Elmt_Id(n);
611 
612  locExpansion = expList->GetExp(nel);
613 
614  for (j = 0; j < locExpansion->GetNfaces(); ++j)
615  {
616  dof = locExpansion->GetFaceIntNcoeffs(j);
617  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
618 
619  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
620 
621  if(faceDirMap.count(meshFaceId)==0)
622  {
623  if(uniqueFaceMap.count(meshFaceId)==0 && dof > 0)
624  {
625  uniqueFaceMap[meshFaceId]=facematrixlocation;
626 
627  facemodeoffset[facematrixlocation]=dof*dof;
628 
629  faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
630 
631  ntotalfaceentries+=dof*dof;
632 
633  n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
634  }
635  nlocalNonDirFaces+=dof*dof;
636  }
637 
638  }
639  }
640 
641  m_comm = expList->GetComm();
642  m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
643  m_comm->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
644 
645  //Allocate arrays for block to universal map (number of expansions * p^2)
646  Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
647  Array<OneD, long> FaceBlockToUniversalMap(ntotalfaceentries,-1);
648 
649  Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
650  Array<OneD, int> localFaceToGlobalMatrixMap(nlocalNonDirFaces,-1);
651 
652  //Allocate arrays to store matrices (number of expansions * p^2)
653  Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
654  Array<OneD, NekDouble> FaceBlockArray(nlocalNonDirFaces,-1);
655 
656  int edgematrixoffset=0;
657  int facematrixoffset=0;
658  int vGlobal;
659 
660  for(n=0; n < n_exp; ++n)
661  {
662  nel = expList->GetOffset_Elmt_Id(n);
663 
664  locExpansion = expList->GetExp(nel);
665 
666  //loop over the edges of the expansion
667  for(j = 0; j < locExpansion->GetNedges(); ++j)
668  {
669  //get mesh edge id
670  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
671 
672  nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
673 
674  if(edgeDirMap.count(meshEdgeId)==0)
675  {
676  for(k=0; k<nedgemodes*nedgemodes; ++k)
677  {
678  vGlobal=edgeglobaloffset[uniqueEdgeMap[meshEdgeId]]+k;
679 
680 
681  localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
682 
683  EdgeBlockToUniversalMap[vGlobal]
684  = meshEdgeId * maxEdgeDof * maxEdgeDof + k + 1;
685  }
686  edgematrixoffset+=nedgemodes*nedgemodes;
687  }
688  }
689 
690  //loop over the faces of the expansion
691  for(j = 0; j < locExpansion->GetNfaces(); ++j)
692  {
693  //get mesh face id
694  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
695 
696  nfacemodes = locExpansion->GetFaceIntNcoeffs(j);
697 
698  //Check if face is has dirichlet values
699  if(faceDirMap.count(meshFaceId)==0)
700  {
701  for(k=0; k<nfacemodes*nfacemodes; ++k)
702  {
703  vGlobal=faceglobaloffset[uniqueFaceMap[meshFaceId]]+k;
704 
705  localFaceToGlobalMatrixMap[facematrixoffset+k]
706  = vGlobal;
707 
708  FaceBlockToUniversalMap[vGlobal]
709  = meshFaceId * maxFaceDof * maxFaceDof + k + 1;
710  }
711  facematrixoffset+=nfacemodes*nfacemodes;
712  }
713  }
714  }
715 
716  edgematrixoffset=0;
717  facematrixoffset=0;
718 
719  m_blkMat = MemoryManager<DNekBlkMat>
720  ::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
721 
722  //Here we loop over the expansion and build the block low energy
723  //preconditioner as well as the block versions of the transformation
724  //matrices.
725  for(cnt=n=0; n < n_exp; ++n)
726  {
727  nel = expList->GetOffset_Elmt_Id(n);
728 
729  locExpansion = expList->GetExp(nel);
730 
731  nVerts=locExpansion->GetGeom()->GetNumVerts();
732  nEdges=locExpansion->GetGeom()->GetNumEdges();
733  nFaces=locExpansion->GetGeom()->GetNumFaces();
734 
735  //Get statically condensed matrix
736  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
737 
738  //Extract boundary block (elemental S1)
739  bnd_mat=loc_mat->GetBlock(0,0);
740 
741  //offset by number of rows
742  offset = bnd_mat->GetRows();
743 
744  DNekScalMat &S=(*bnd_mat);
745 
746  //loop over vertices of the element and return the vertex map
747  //for each vertex
748  for (v=0; v<nVerts; ++v)
749  {
750  vMap1=locExpansion->GetVertexMap(v);
751 
752  //Get vertex map
753  globalrow = m_locToGloMap->
754  GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
755 
756  if(globalrow >= 0)
757  {
758  for (m=0; m<nVerts; ++m)
759  {
760  vMap2=locExpansion->GetVertexMap(m);
761 
762  //global matrix location (without offset due to
763  //dirichlet values)
764  globalcol = m_locToGloMap->
765  GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
766 
767  //offset for dirichlet conditions
768  if (globalcol == globalrow)
769  {
770  meshVertId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetVid(v);
771 
772  //modal connectivity between elements
773  sign1 = m_locToGloMap->
774  GetLocalToGlobalBndSign(cnt + vMap1);
775  sign2 = m_locToGloMap->
776  GetLocalToGlobalBndSign(cnt + vMap2);
777 
778  vertArray[globalrow]
779  += sign1*sign2*S(vMap1,vMap2);
780 
781  VertBlockToUniversalMap[globalrow]
782  = meshVertId * maxEdgeDof * maxEdgeDof + 1;
783  }
784  }
785  }
786  }
787 
788  //loop over edges of the element and return the edge map
789  for (eid=0; eid<nEdges; ++eid)
790  {
791  nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
792 
793  DNekMatSharedPtr locMat =
794  MemoryManager<DNekMat>::AllocateSharedPtr
795  (nedgemodes,nedgemodes,zero,storage);
796 
797  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(eid);
798  Array<OneD, unsigned int> edgemodearray = locExpansion->GetEdgeInverseBoundaryMap(eid);
799 
800  if(edgeDirMap.count(meshEdgeId)==0)
801  {
802  for (v=0; v<nedgemodes; ++v)
803  {
804  eMap1=edgemodearray[v];
805 
806  for (m=0; m<nedgemodes; ++m)
807  {
808  eMap2=edgemodearray[m];
809 
810  //modal connectivity between elements
811  sign1 = m_locToGloMap->
812  GetLocalToGlobalBndSign(cnt + eMap1);
813  sign2 = m_locToGloMap->
814  GetLocalToGlobalBndSign(cnt + eMap2);
815 
816  NekDouble globalEdgeValue = sign1*sign2*S(eMap1,eMap2);
817 
818  EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=globalEdgeValue;
819  }
820  }
821  edgematrixoffset+=nedgemodes*nedgemodes;
822  }
823  }
824 
825  //loop over faces of the element and return the face map
826  for (fid=0; fid<nFaces; ++fid)
827  {
828  nfacemodes = locExpansion->GetFaceIntNcoeffs(fid);
829 
830  DNekMatSharedPtr locMat =
831  MemoryManager<DNekMat>::AllocateSharedPtr
832  (nfacemodes,nfacemodes,zero,storage);
833 
834  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(fid);
835 
836  Array<OneD, unsigned int> facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid);
837 
838  if(faceDirMap.count(meshFaceId)==0)
839  {
840  for (v=0; v<nfacemodes; ++v)
841  {
842  fMap1=facemodearray[v];
843 
844  for (m=0; m<nfacemodes; ++m)
845  {
846  fMap2=facemodearray[m];
847 
848  //modal connectivity between elements
849  sign1 = m_locToGloMap->
850  GetLocalToGlobalBndSign(cnt + fMap1);
851  sign2 = m_locToGloMap->
852  GetLocalToGlobalBndSign(cnt + fMap2);
853 
854  // Get the face-face value from the low energy matrix (S2)
855  NekDouble globalFaceValue = sign1*sign2*S(fMap1,fMap2);
856 
857  //local face value to global face value
858  FaceBlockArray[facematrixoffset+v*nfacemodes+m]=globalFaceValue;
859  }
860  }
861  facematrixoffset+=nfacemodes*nfacemodes;
862  }
863  }
864 
865  //offset for the expansion
866  cnt+=offset;
867  }
868 
869  //Assemble edge matrices of each process
870  Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries);
871  Vmath::Zero(ntotaledgeentries, GlobalEdgeBlock.get(), 1);
872  Vmath::Assmb(EdgeBlockArray.num_elements(),
873  EdgeBlockArray.get(),
874  localEdgeToGlobalMatrixMap.get(),
875  GlobalEdgeBlock.get());
876 
877  //Assemble face matrices of each process
878  Array<OneD, NekDouble> GlobalFaceBlock(ntotalfaceentries);
879  Vmath::Zero(ntotalfaceentries, GlobalFaceBlock.get(), 1);
880  Vmath::Assmb(FaceBlockArray.num_elements(),
881  FaceBlockArray.get(),
882  localFaceToGlobalMatrixMap.get(),
883  GlobalFaceBlock.get());
884 
885  //Exchange vertex data over different processes
886  if(nNonDirVerts!=0)
887  {
888  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
889  Gs::Gather(vertArray, Gs::gs_add, tmp);
890  }
891 
892  //Exchange edge data over different processes
893  Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
894  Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
895 
896  //Exchange face data over different processes
897  Gs::gs_data *tmp2 = Gs::Init(FaceBlockToUniversalMap, m_comm);
898  Gs::Gather(GlobalFaceBlock, Gs::gs_add, tmp2);
899 
900  // Populate vertex block
901  for (int i = 0; i < nNonDirVerts; ++i)
902  {
903  VertBlk->SetValue(i,i,1.0/vertArray[i]);
904  }
905 
906  //Set the first block to be the diagonal of the vertex space
907  m_blkMat->SetBlock(0,0, VertBlk);
908 
909  offset=0;
910  //Build the edge matrices from the vector
911  for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
912  {
913  DNekMatSharedPtr gmat =
914  MemoryManager<DNekMat>::AllocateSharedPtr
915  (nedgemodes,nedgemodes,zero,storage);
916 
917  for (v=0; v<nedgemodes; ++v)
918  {
919  for (m=0; m<nedgemodes; ++m)
920  {
921  NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
922  gmat->SetValue(v,m,EdgeValue);
923  }
924  }
925 
926  m_blkMat->SetBlock(1+loc,1+loc, gmat);
927 
928  offset+=edgemodeoffset[loc];
929  }
930 
931  offset=0;
932  //Build the face matrices from the vector
933  for(int loc=0; loc<nNonDirFaceIDs; ++loc)
934  {
935  nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
936 
937  DNekMatSharedPtr gmat =
938  MemoryManager<DNekMat>::AllocateSharedPtr
939  (nfacemodes,nfacemodes,zero,storage);
940 
941  for (v=0; v<nfacemodes; ++v)
942  {
943  for (m=0; m<nfacemodes; ++m)
944  {
945  NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
946  gmat->SetValue(v,m,FaceValue);
947  }
948  }
949 
950  m_blkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
951 
952  offset+=facemodeoffset[loc];
953  }
954 
955 
956  int totblks=m_blkMat->GetNumberOfBlockRows();
957  for (i=1; i< totblks; ++i)
958  {
959  unsigned int nmodes=m_blkMat->GetNumberOfRowsInBlockRow(i);
960  DNekMatSharedPtr tmp_mat =
961  MemoryManager<DNekMat>::AllocateSharedPtr
962  (nmodes,nmodes,zero,storage);
963 
964  tmp_mat=m_blkMat->GetBlock(i,i);
965  tmp_mat->Invert();
966 
967  m_blkMat->SetBlock(i,i,tmp_mat);
968  }
969  }
970 
971  /**
972  *
973  */
974  void PreconditionerBlock::v_DoPreconditioner(
975  const Array<OneD, NekDouble>& pInput,
976  Array<OneD, NekDouble>& pOutput)
977  {
978  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
979  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
980  int nNonDir = nGlobal-nDir;
981  DNekBlkMat &M = (*m_blkMat);
982 
983  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
984  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
985  z = M * r;
986  }
987  }
988 }
989 
990 
991 
992 
993 
994 
Generated on Sun Mar 15 2015 00:11:33 for Nektar++ by   doxygen 1.8.1.2