Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AssemblyMap.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AssemblyMap.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: Assembly (e.g. local to global) base mapping routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/AssemblyMap/AssemblyMap.h>
37 
38 namespace Nektar
39 {
40  namespace MultiRegions
41  {
42  /**
43  * @class AssemblyMap
44  * This class acts as a base class for constructing mappings between
45  * local, global and boundary degrees of freedom. It holds the storage
46  * for the maps and provides the accessors needed to retrieve them.
47  *
48  * There are two derived classes: AssemblyMapCG and
49  * AssemblyMapDG. These perform the actual construction of the
50  * maps within their specific contexts.
51  *
52  */
53 
54  /// Rounds a double precision number to an integer.
55  int RoundNekDoubleToInt(NekDouble x)
56  {
57  return int(x > 0.0 ? x + 0.5 : x - 0.5);
58  }
59 
60  /// Rounds an array of double precision numbers to integers.
61  void RoundNekDoubleToInt(const Array<OneD,const NekDouble> inarray, Array<OneD,int> outarray)
62  {
63  int size = inarray.num_elements();
64  ASSERTL1(outarray.num_elements()>=size,"Array sizes not compatible");
65 
66  NekDouble x;
67  for(int i = 0; i < size; i++)
68  {
69  x = inarray[i];
70  outarray[i] = int(x > 0.0 ? x + 0.5 : x - 0.5);
71  }
72  }
73 
74  /**
75  * Initialises an empty mapping.
76  */
77  AssemblyMap::AssemblyMap():
78  m_session(),
79  m_comm(),
80  m_hash(0),
81  m_numLocalBndCoeffs(0),
82  m_numGlobalBndCoeffs(0),
83  m_numLocalDirBndCoeffs(0),
84  m_numGlobalDirBndCoeffs(0),
85  m_solnType(eNoSolnType),
86  m_bndSystemBandWidth(0),
87  m_successiveRHS(0),
88  m_gsh(0),
89  m_bndGsh(0)
90  {
91  }
92 
93  AssemblyMap::AssemblyMap(
94  const LibUtilities::SessionReaderSharedPtr &pSession,
95  const std::string variable):
96  m_session(pSession),
97  m_comm(pSession->GetComm()),
98  m_hash(0),
99  m_numLocalBndCoeffs(0),
100  m_numGlobalBndCoeffs(0),
101  m_numLocalDirBndCoeffs(0),
102  m_numGlobalDirBndCoeffs(0),
103  m_bndSystemBandWidth(0),
104  m_successiveRHS(0),
105  m_gsh(0),
106  m_bndGsh(0)
107  {
108  // Default value from Solver Info
109  m_solnType = pSession->GetSolverInfoAsEnum<GlobalSysSolnType>(
110  "GlobalSysSoln");
111  m_preconType = pSession->GetSolverInfoAsEnum<PreconditionerType>(
112  "Preconditioner");
113 
114  // Override values with data from GlobalSysSolnInfo section
115  if(pSession->DefinesGlobalSysSolnInfo(variable, "GlobalSysSoln"))
116  {
117  std::string sysSoln = pSession->GetGlobalSysSolnInfo(variable,
118  "GlobalSysSoln");
119  m_solnType = pSession->GetValueAsEnum<GlobalSysSolnType>(
120  "GlobalSysSoln", sysSoln);
121  }
122 
123  if(pSession->DefinesGlobalSysSolnInfo(variable, "Preconditioner"))
124  {
125  std::string precon = pSession->GetGlobalSysSolnInfo(variable,
126  "Preconditioner");
127  m_preconType = pSession->GetValueAsEnum<PreconditionerType>(
128  "Preconditioner", precon);
129  }
130 
131  if(pSession->DefinesGlobalSysSolnInfo(variable,
132  "IterativeSolverTolerance"))
133  {
134  m_iterativeTolerance = boost::lexical_cast<NekDouble>(
135  pSession->GetGlobalSysSolnInfo(variable,
136  "IterativeSolverTolerance").c_str());
137  }
138  else
139  {
140  pSession->LoadParameter("IterativeSolverTolerance",
141  m_iterativeTolerance,
142  NekConstants::kNekIterativeTol);
143  }
144 
145 
146  if(pSession->DefinesGlobalSysSolnInfo(variable,"SuccessiveRHS"))
147  {
148  m_successiveRHS = boost::lexical_cast<int>(
149  pSession->GetGlobalSysSolnInfo(variable,
150  "SuccessiveRHS").c_str());
151  }
152  else
153  {
154  pSession->LoadParameter("SuccessiveRHS",
155  m_successiveRHS,0);
156  }
157 
158  }
159 
160  /**
161  * Create a new level of mapping using the information in
162  * multiLevelGraph and performing the following steps:
163  */
164  AssemblyMap::AssemblyMap(
165  AssemblyMap* oldLevelMap,
166  const BottomUpSubStructuredGraphSharedPtr& multiLevelGraph):
167  m_session(oldLevelMap->m_session),
168  m_comm(oldLevelMap->GetComm()),
169  m_hash(0),
170  m_solnType(oldLevelMap->m_solnType),
171  m_preconType(oldLevelMap->m_preconType),
172  m_iterativeTolerance(oldLevelMap->m_iterativeTolerance),
173  m_successiveRHS(oldLevelMap->m_successiveRHS),
174  m_gsh(oldLevelMap->m_gsh),
175  m_bndGsh(oldLevelMap->m_bndGsh),
176  m_lowestStaticCondLevel(oldLevelMap->m_lowestStaticCondLevel)
177  {
178  int i;
179  int j;
180  int cnt;
181 
182  //--------------------------------------------------------------
183  // -- Extract information from the input argument
184  int numGlobalBndCoeffsOld = oldLevelMap->GetNumGlobalBndCoeffs();
185  int numGlobalDirBndCoeffsOld = oldLevelMap->GetNumGlobalDirBndCoeffs();
186  int numLocalBndCoeffsOld = oldLevelMap->GetNumLocalBndCoeffs();
187  int numLocalDirBndCoeffsOld = oldLevelMap->GetNumLocalDirBndCoeffs();
188  bool signChangeOld = oldLevelMap->GetSignChange();
189 
190  int staticCondLevelOld = oldLevelMap->GetStaticCondLevel();
191  int numPatchesOld = oldLevelMap->GetNumPatches();
192  GlobalSysSolnType solnTypeOld = oldLevelMap->GetGlobalSysSolnType();
193  const Array<OneD, const unsigned int>& numLocalBndCoeffsPerPatchOld = oldLevelMap->GetNumLocalBndCoeffsPerPatch();
194  //--------------------------------------------------------------
195 
196  //--------------------------------------------------------------
197  int newLevel = staticCondLevelOld+1;
198  /** - STEP 1: setup a mask array to determine to which patch
199  * of the new level every patch of the current
200  * level belongs. To do so we make four arrays,
201  * #gloPatchMask, #globHomPatchMask,
202  * #locPatchMask_NekDouble and #locPatchMask.
203  * These arrays are then used to check which local
204  * dofs of the old level belong to which patch of
205  * the new level
206  */
207  Array<OneD, NekDouble> globPatchMask (numGlobalBndCoeffsOld,-1.0);
208  Array<OneD, NekDouble> globHomPatchMask (globPatchMask+numGlobalDirBndCoeffsOld);
209  Array<OneD, NekDouble> locPatchMask_NekDouble(numLocalBndCoeffsOld,-3.0);
210  Array<OneD, int> locPatchMask (numLocalBndCoeffsOld);
211 
212  // Fill the array globPatchMask as follows:
213  // - The first part (i.e. the glob bnd dofs) is filled with the
214  // value -1
215  // - The second part (i.e. the glob interior dofs) is numbered
216  // according to the patch it belongs to (i.e. dofs in first block
217  // all are numbered 0, the second block numbered are 1, etc...)
218  multiLevelGraph->MaskPatches(newLevel,globHomPatchMask);
219 
220  // Map from Global Dofs to Local Dofs
221  // As a result, we know for each local dof whether
222  // it is mapped to the boundary of the next level, or to which
223  // patch. Based upon this, we can than later associate every patch
224  // of the current level with a patch in the next level.
225  oldLevelMap->GlobalToLocalBndWithoutSign(globPatchMask,locPatchMask_NekDouble);
226 
227  // Convert the result to an array of integers rather than doubles
228  RoundNekDoubleToInt(locPatchMask_NekDouble,locPatchMask);
229 
230  /** - STEP 2: We calculate how many local bnd dofs of the
231  * old level belong to the boundaries of each patch at
232  * the new level. We need this information to set up the
233  * mapping between different levels.
234  */
235 
236  // Retrieve the number of patches at the next level
237  int numPatchesWithIntNew = multiLevelGraph->GetNpatchesWithInterior(newLevel);
238  int numPatchesNew = numPatchesWithIntNew;
239 
240  // Allocate memory to store the number of local dofs associated to
241  // each of elemental boundaries of these patches
242  map<int, int> numLocalBndCoeffsPerPatchNew;
243  for(int i = 0; i < numPatchesNew; i++)
244  {
245  numLocalBndCoeffsPerPatchNew[i] = 0;
246  }
247 
248  int minval;
249  int maxval;
250  int curPatch;
251  for(i = cnt = 0; i < numPatchesOld; i++)
252  {
253  // For every patch at the current level, the mask array
254  // locPatchMask should be filled with
255  // - the same (positive) number for each entry (which will
256  // correspond to the patch at the next level it belongs to)
257  // - the same (positive) number for each entry, except some
258  // entries that are -1 (the enties correspond to -1, will be
259  // mapped to the local boundary of the next level patch given
260  // by the positive number)
261  // - -1 for all entries. In this case, we will make an
262  // additional patch only consisting of boundaries at the next
263  // level
264  minval = *min_element(&locPatchMask[cnt],
265  &locPatchMask[cnt]+numLocalBndCoeffsPerPatchOld[i]);
266  maxval = *max_element(&locPatchMask[cnt],
267  &locPatchMask[cnt]+numLocalBndCoeffsPerPatchOld[i]);
268  ASSERTL0((minval==maxval)||(minval==-1),"These values should never be the same");
269 
270  if(maxval == -1)
271  {
272  curPatch = numPatchesNew;
273  numLocalBndCoeffsPerPatchNew[curPatch] = 0;
274  numPatchesNew++;
275  }
276  else
277  {
278  curPatch = maxval;
279  }
280 
281  for(j = 0; j < numLocalBndCoeffsPerPatchOld[i]; j++ )
282  {
283  ASSERTL0((locPatchMask[cnt]==maxval)||(locPatchMask[cnt]==minval),
284  "These values should never be the same");
285  if(locPatchMask[cnt] == -1)
286  {
287  ++numLocalBndCoeffsPerPatchNew[curPatch];
288  }
289  cnt++;
290  }
291  }
292 
293  // Count how many local dofs of the old level are mapped
294  // to the local boundary dofs of the new level
295  m_numLocalBndCoeffs = 0;
296  m_numPatches = numLocalBndCoeffsPerPatchNew.size();
297  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
298  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches,0u);
299  for(int i = 0; i < m_numPatches; i++)
300  {
301  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) numLocalBndCoeffsPerPatchNew[i];
302  m_numLocalBndCoeffs += numLocalBndCoeffsPerPatchNew[i];
303  }
304  multiLevelGraph->GetNintDofsPerPatch(newLevel,m_numLocalIntCoeffsPerPatch);
305 
306  // Also initialise some more data members
307  m_solnType = solnTypeOld;
308  ASSERTL1(m_solnType==eDirectMultiLevelStaticCond
309  ||m_solnType==eIterativeMultiLevelStaticCond
310  ||m_solnType==eXxtMultiLevelStaticCond
311  ||m_solnType==ePETScMultiLevelStaticCond,
312  "This method should only be called for in "
313  "case of multi-level static condensation.");
314  m_staticCondLevel = newLevel;
315  m_signChange = signChangeOld;
316  m_numLocalDirBndCoeffs = numLocalDirBndCoeffsOld;
317  m_numGlobalDirBndCoeffs = numGlobalDirBndCoeffsOld;
318  m_numGlobalBndCoeffs = multiLevelGraph->GetInteriorOffset(newLevel) +
319  m_numGlobalDirBndCoeffs;
320  m_numGlobalCoeffs = multiLevelGraph->GetNumGlobalDofs(newLevel) +
321  m_numGlobalDirBndCoeffs;
322  m_localToGlobalBndMap = Array<OneD,int>(m_numLocalBndCoeffs);
323  if(m_signChange)
324  {
325  m_localToGlobalBndSign = Array<OneD,NekDouble>(m_numLocalBndCoeffs);
326  }
327 
328  m_patchMapFromPrevLevel = MemoryManager<PatchMap>::AllocateSharedPtr(numLocalBndCoeffsOld);
329 
330  m_globalToUniversalBndMap = Array<OneD, int>(
331  m_numGlobalBndCoeffs, oldLevelMap->GetGlobalToUniversalBndMap());
332  m_globalToUniversalBndMapUnique = Array<OneD, int>(
333  m_numGlobalBndCoeffs, oldLevelMap->GetGlobalToUniversalBndMapUnique());
334 
335  // Set up an offset array that denotes the offset of the local
336  // boundary degrees of freedom of the next level
337  Array<OneD, int> numLocalBndCoeffsPerPatchOffset(m_numPatches+1,0);
338  for(int i = 1; i < m_numPatches+1; i++)
339  {
340  numLocalBndCoeffsPerPatchOffset[i] += numLocalBndCoeffsPerPatchOffset[i-1] + numLocalBndCoeffsPerPatchNew[i-1];
341  }
342 
343  int additionalPatchCnt = numPatchesWithIntNew;
344  int newid;
345  int blockid;
346  bool isBndDof;
347  NekDouble sign;
348  Array<OneD, int> bndDofPerPatchCnt(m_numPatches,0);
349  for(i = cnt = 0; i < numPatchesOld; i++)
350  {
351  minval = *min_element(&locPatchMask[cnt],&locPatchMask[cnt]+numLocalBndCoeffsPerPatchOld[i]);
352  maxval = *max_element(&locPatchMask[cnt],&locPatchMask[cnt]+numLocalBndCoeffsPerPatchOld[i]);
353  ASSERTL0((minval==maxval)||(minval==-1),"These values should never be the same");
354 
355  if(maxval == -1)
356  {
357  curPatch = additionalPatchCnt;
358  additionalPatchCnt++;
359  }
360  else
361  {
362  curPatch = maxval;
363  }
364 
365  for(j = 0; j < numLocalBndCoeffsPerPatchOld[i]; j++ )
366  {
367  ASSERTL0((locPatchMask[cnt]==maxval)||(locPatchMask[cnt]==minval),
368  "These values should never be the same");
369 
370  sign = oldLevelMap->GetLocalToGlobalBndSign(cnt);
371 
372  if(locPatchMask[cnt] == -1)
373  {
374  newid = numLocalBndCoeffsPerPatchOffset[curPatch];
375 
376  m_localToGlobalBndMap[newid] = oldLevelMap->GetLocalToGlobalBndMap(cnt);
377  if(m_signChange)
378  {
379  m_localToGlobalBndSign[ newid ] = sign;
380  }
381 
382  blockid = bndDofPerPatchCnt[curPatch];
383  isBndDof = true;
384 
385 
386  numLocalBndCoeffsPerPatchOffset[curPatch]++;
387  bndDofPerPatchCnt[curPatch]++;
388  }
389  else
390  {
391  newid = oldLevelMap->GetLocalToGlobalBndMap(cnt) -
392  m_numGlobalBndCoeffs+m_numLocalBndCoeffs;
393 
394  blockid = oldLevelMap->GetLocalToGlobalBndMap(cnt)-
395  m_numGlobalDirBndCoeffs - multiLevelGraph->GetInteriorOffset(newLevel,curPatch);
396  isBndDof = false;
397  }
398 
399  sign = isBndDof?1.0:sign;
400 
401  m_patchMapFromPrevLevel->SetPatchMap(cnt,curPatch,blockid,isBndDof,sign);
402 
403  cnt++;
404  }
405  }
406  CalculateBndSystemBandWidth();
407 
408  // Postprocess the computed information - Update the old
409  // level with the mapping to new evel
410  // oldLevelMap->SetLocalBndToNextLevelMap(oldLocalBndToNewLevelMap,oldLocalBndToNewLevelSign);
411  // - Construct the next level mapping object
412  if(m_staticCondLevel < (multiLevelGraph->GetNlevels()-1))
413  {
414  m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::AllocateSharedPtr(this,multiLevelGraph);
415  }
416  }
417 
418 
419  AssemblyMap::~AssemblyMap(void)
420  {
421  }
422 
423 
424  /**
425  * The bandwidth calculated corresponds to what is referred to as
426  * half-bandwidth. If the elements of the matrix are designated as
427  * a_ij, it corresponds to the maximum value of |i-j| for non-zero
428  * a_ij. As a result, the value also corresponds to the number of
429  * sub- or super-diagonals.
430  *
431  * The bandwith can be calculated elementally as it corresponds to the
432  * maximal elemental bandwith (i.e. the maximal difference in global
433  * DOF index for every element).
434  *
435  * We here calculate the bandwith of the global boundary system (as
436  * used for static condensation).
437  */
438  void AssemblyMap::CalculateBndSystemBandWidth()
439  {
440  int i,j;
441  int cnt = 0;
442  int locSize;
443  int maxId;
444  int minId;
445  int bwidth = -1;
446  for(i = 0; i < m_numPatches; ++i)
447  {
448  locSize = m_numLocalBndCoeffsPerPatch[i];
449  maxId = -1;
450  minId = m_numLocalBndCoeffs+1;
451  for(j = 0; j < locSize; j++)
452  {
453  if(m_localToGlobalBndMap[cnt+j] >= m_numGlobalDirBndCoeffs)
454  {
455  if(m_localToGlobalBndMap[cnt+j] > maxId)
456  {
457  maxId = m_localToGlobalBndMap[cnt+j];
458  }
459 
460  if(m_localToGlobalBndMap[cnt+j] < minId)
461  {
462  minId = m_localToGlobalBndMap[cnt+j];
463  }
464  }
465  }
466  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
467 
468  cnt+=locSize;
469  }
470 
471  m_bndSystemBandWidth = bwidth;
472  }
473 
474 
475  int AssemblyMap::v_GetLocalToGlobalMap(const int i) const
476  {
477  ASSERTL0(false, "Not defined for this type of mapping.");
478  return 0;
479  }
480 
481  int AssemblyMap::v_GetGlobalToUniversalMap(const int i) const
482  {
483  ASSERTL0(false, "Not defined for this type of mapping.");
484  return 0;
485  }
486 
487  int AssemblyMap::v_GetGlobalToUniversalMapUnique(const int i) const
488  {
489  ASSERTL0(false, "Not defined for this type of mapping.");
490  return 0;
491  }
492 
493  const Array<OneD,const int>& AssemblyMap::v_GetLocalToGlobalMap()
494  {
495  ASSERTL0(false, "Not defined for this type of mapping.");
496  static Array<OneD,const int> result;
497  return result;
498  }
499 
500  const Array<OneD, const int>& AssemblyMap::v_GetGlobalToUniversalMap()
501  {
502  ASSERTL0(false, "Not defined for this type of mapping.");
503  static Array<OneD, const int> result;
504  return result;
505  }
506 
507  const Array<OneD, const int>& AssemblyMap::v_GetGlobalToUniversalMapUnique()
508  {
509  ASSERTL0(false, "Not defined for this type of mapping.");
510  static Array<OneD, const int> result;
511  return result;
512  }
513 
514  NekDouble AssemblyMap::v_GetLocalToGlobalSign(const int i) const
515  {
516  ASSERTL0(false, "Not defined for this type of mapping.");
517  return 0.0;
518  }
519 
520  const Array<OneD, NekDouble>& AssemblyMap::v_GetLocalToGlobalSign() const
521  {
522  ASSERTL0(false, "Not defined for this type of mapping.");
523  static Array<OneD, NekDouble> result;
524  return result;
525  }
526 
527  void AssemblyMap::v_LocalToGlobal(
528  const Array<OneD, const NekDouble>& loc,
529  Array<OneD, NekDouble>& global) const
530  {
531  ASSERTL0(false, "Not defined for this type of mapping.");
532  }
533 
534  void AssemblyMap::v_LocalToGlobal(
535  const NekVector<NekDouble>& loc,
536  NekVector< NekDouble>& global) const
537  {
538  ASSERTL0(false, "Not defined for this type of mapping.");
539  }
540 
541  void AssemblyMap::v_GlobalToLocal(
542  const Array<OneD, const NekDouble>& global,
543  Array<OneD, NekDouble>& loc) const
544  {
545  ASSERTL0(false, "Not defined for this type of mapping.");
546  }
547 
548  void AssemblyMap::v_GlobalToLocal(
549  const NekVector<NekDouble>& global,
550  NekVector< NekDouble>& loc) const
551  {
552  ASSERTL0(false, "Not defined for this type of mapping.");
553  }
554 
555  void AssemblyMap::v_Assemble(
556  const Array<OneD, const NekDouble> &loc,
557  Array<OneD, NekDouble> &global) const
558  {
559  ASSERTL0(false, "Not defined for this type of mapping.");
560  }
561 
562  void AssemblyMap::v_Assemble(
563  const NekVector<NekDouble>& loc,
564  NekVector< NekDouble>& global) const
565  {
566  ASSERTL0(false, "Not defined for this type of mapping.");
567  }
568 
569  void AssemblyMap::v_UniversalAssemble(
570  Array<OneD, NekDouble>& pGlobal) const
571  {
572  // Do nothing here since multi-level static condensation uses a
573  // AssemblyMap and thus will call this routine in serial.
574  }
575 
576  void AssemblyMap::v_UniversalAssemble(
577  NekVector< NekDouble>& pGlobal) const
578  {
579  // Do nothing here since multi-level static condensation uses a
580  // AssemblyMap and thus will call this routine in serial.
581  }
582 
583  void AssemblyMap::v_UniversalAssemble(
584  Array<OneD, NekDouble>& pGlobal,
585  int offset) const
586  {
587  // Do nothing here since multi-level static condensation uses a
588  // AssemblyMap and thus will call this routine in serial.
589  }
590 
591  int AssemblyMap::v_GetFullSystemBandWidth() const
592  {
593  ASSERTL0(false, "Not defined for this type of mapping.");
594  return 0;
595  }
596 
597  int AssemblyMap::v_GetNumNonDirVertexModes() const
598  {
599  ASSERTL0(false, "Not defined for this type of mapping.");
600  return 0;
601  }
602 
603  int AssemblyMap::v_GetNumNonDirEdgeModes() const
604  {
605  ASSERTL0(false, "Not defined for this type of mapping.");
606  return 0;
607  }
608 
609  int AssemblyMap::v_GetNumNonDirFaceModes() const
610  {
611  ASSERTL0(false, "Not defined for this type of mapping.");
612  return 0;
613  }
614 
615  int AssemblyMap::v_GetNumDirEdges() const
616  {
617  ASSERTL0(false, "Not defined for this type of mapping.");
618  return 0;
619  }
620 
621  int AssemblyMap::v_GetNumDirFaces() const
622  {
623  ASSERTL0(false, "Not defined for this type of mapping.");
624  return 0;
625  }
626 
627  int AssemblyMap::v_GetNumNonDirEdges() const
628  {
629  ASSERTL0(false, "Not defined for this type of mapping.");
630  return 0;
631  }
632 
633  int AssemblyMap::v_GetNumNonDirFaces() const
634  {
635  ASSERTL0(false, "Not defined for this type of mapping.");
636  return 0;
637  }
638 
639  const Array<OneD, const int>& AssemblyMap::v_GetExtraDirEdges()
640  {
641  ASSERTL0(false, "Not defined for this type of mapping.");
642  static Array<OneD, const int> result;
643  return result;
644  }
645 
646  boost::shared_ptr<AssemblyMap> AssemblyMap::v_LinearSpaceMap(
647  const ExpList &locexp, GlobalSysSolnType solnType)
648  {
649  ASSERTL0(false, "Not defined for this sub class");
650  static boost::shared_ptr<AssemblyMap> result;
651  return result;
652  }
653 
654  LibUtilities::CommSharedPtr AssemblyMap::GetComm()
655  {
656  return m_comm;
657  }
658 
659  size_t AssemblyMap::GetHash() const
660  {
661  return m_hash;
662  }
663 
664  int AssemblyMap::GetLocalToGlobalMap(const int i) const
665  {
666  return v_GetLocalToGlobalMap(i);
667  }
668 
669  int AssemblyMap::GetGlobalToUniversalMap(const int i) const
670  {
671  return v_GetGlobalToUniversalMap(i);
672  }
673 
674  int AssemblyMap::GetGlobalToUniversalMapUnique(const int i) const
675  {
676  return v_GetGlobalToUniversalMapUnique(i);
677  }
678 
679  const Array<OneD,const int>& AssemblyMap::GetLocalToGlobalMap()
680  {
681  return v_GetLocalToGlobalMap();
682  }
683 
684  const Array<OneD, const int>& AssemblyMap::GetGlobalToUniversalMap()
685  {
686  return v_GetGlobalToUniversalMap();
687  }
688 
689  const Array<OneD, const int>& AssemblyMap::GetGlobalToUniversalMapUnique()
690  {
691  return v_GetGlobalToUniversalMapUnique();
692  }
693 
694  NekDouble AssemblyMap::GetLocalToGlobalSign(const int i) const
695  {
696  return v_GetLocalToGlobalSign(i);
697  }
698 
699  const Array<OneD, NekDouble>& AssemblyMap::GetLocalToGlobalSign() const
700  {
701  return v_GetLocalToGlobalSign();
702  }
703 
704  void AssemblyMap::LocalToGlobal(
705  const Array<OneD, const NekDouble>& loc,
706  Array<OneD, NekDouble>& global) const
707  {
708  v_LocalToGlobal(loc,global);
709  }
710 
711  void AssemblyMap::LocalToGlobal(
712  const NekVector<NekDouble>& loc,
713  NekVector< NekDouble>& global) const
714  {
715  v_LocalToGlobal(loc,global);
716  }
717 
718  void AssemblyMap::GlobalToLocal(
719  const Array<OneD, const NekDouble>& global,
720  Array<OneD, NekDouble>& loc) const
721  {
722  v_GlobalToLocal(global,loc);
723  }
724 
725  void AssemblyMap::GlobalToLocal(
726  const NekVector<NekDouble>& global,
727  NekVector< NekDouble>& loc) const
728  {
729  v_GlobalToLocal(global,loc);
730  }
731 
732  void AssemblyMap::Assemble(
733  const Array<OneD, const NekDouble> &loc,
734  Array<OneD, NekDouble> &global) const
735  {
736  v_Assemble(loc,global);
737  }
738 
739  void AssemblyMap::Assemble(
740  const NekVector<NekDouble>& loc,
741  NekVector< NekDouble>& global) const
742  {
743  v_Assemble(loc,global);
744  }
745 
746  void AssemblyMap::UniversalAssemble(
747  Array<OneD, NekDouble>& pGlobal) const
748  {
749  v_UniversalAssemble(pGlobal);
750  }
751 
752  void AssemblyMap::UniversalAssemble(
753  NekVector< NekDouble>& pGlobal) const
754  {
755  v_UniversalAssemble(pGlobal);
756  }
757 
758  void AssemblyMap::UniversalAssemble(
759  Array<OneD, NekDouble>& pGlobal,
760  int offset) const
761  {
762  v_UniversalAssemble(pGlobal, offset);
763  }
764 
765  int AssemblyMap::GetFullSystemBandWidth() const
766  {
767  return v_GetFullSystemBandWidth();
768  }
769 
770  int AssemblyMap::GetNumNonDirVertexModes() const
771  {
772  return v_GetNumNonDirVertexModes();
773  }
774 
775  int AssemblyMap::GetNumNonDirEdgeModes() const
776  {
777  return v_GetNumNonDirEdgeModes();
778  }
779 
780  int AssemblyMap::GetNumNonDirFaceModes() const
781  {
782  return v_GetNumNonDirFaceModes();
783  }
784 
785  int AssemblyMap::GetNumDirEdges() const
786  {
787  return v_GetNumDirEdges();
788  }
789 
790  int AssemblyMap::GetNumDirFaces() const
791  {
792  return v_GetNumDirFaces();
793  }
794 
795  int AssemblyMap::GetNumNonDirEdges() const
796  {
797  return v_GetNumNonDirEdges();
798  }
799 
800  int AssemblyMap::GetNumNonDirFaces() const
801  {
802  return v_GetNumNonDirFaces();
803  }
804 
805  const Array<OneD, const int>& AssemblyMap::GetExtraDirEdges()
806  {
807  return v_GetExtraDirEdges();
808  }
809 
810  boost::shared_ptr<AssemblyMap> AssemblyMap::LinearSpaceMap(const ExpList &locexp, GlobalSysSolnType solnType)
811  {
812  return v_LinearSpaceMap(locexp, solnType);
813  }
814 
815  int AssemblyMap::GetLocalToGlobalBndMap(const int i) const
816  {
817  return m_localToGlobalBndMap[i];
818  }
819 
820  const Array<OneD,const int>&
821  AssemblyMap::GetLocalToGlobalBndMap(void)
822  {
823  return m_localToGlobalBndMap;
824  }
825 
826  bool AssemblyMap::GetSignChange()
827  {
828  return m_signChange;
829  }
830 
831 
832  Array<OneD, const NekDouble>
833  AssemblyMap::GetLocalToGlobalBndSign(void) const
834  {
835  return m_localToGlobalBndSign;
836  }
837 
838  const Array<OneD, const int>& AssemblyMap::GetGlobalToUniversalBndMap()
839  {
840  return m_globalToUniversalBndMap;
841  }
842 
843  const Array<OneD, const int>& AssemblyMap::GetGlobalToUniversalBndMapUnique()
844  {
845  return m_globalToUniversalBndMapUnique;
846  }
847 
848  NekDouble AssemblyMap::GetLocalToGlobalBndSign(const int i) const
849  {
850  if(m_signChange)
851  {
852  return m_localToGlobalBndSign[i];
853  }
854  else
855  {
856  return 1.0;
857  }
858  }
859 
860 
861  int AssemblyMap::GetBndCondCoeffsToGlobalCoeffsMap(
862  const int i)
863  {
864  return m_bndCondCoeffsToGlobalCoeffsMap[i];
865  }
866 
867 
868  int AssemblyMap::GetBndCondTraceToGlobalTraceMap(
869  const int i)
870  {
871  ASSERTL1(i < m_bndCondTraceToGlobalTraceMap.num_elements(),
872  "Index out of range.");
873  return m_bndCondTraceToGlobalTraceMap[i];
874  }
875 
876  const Array<OneD, const int> &AssemblyMap
877  ::GetBndCondTraceToGlobalTraceMap()
878  {
879  return m_bndCondTraceToGlobalTraceMap;
880  }
881 
882  NekDouble AssemblyMap::GetBndCondCoeffsToGlobalCoeffsSign(const int i)
883  {
884  if(m_signChange)
885  {
886  return m_bndCondCoeffsToGlobalCoeffsSign[i];
887  }
888  else
889  {
890  return 1.0;
891  }
892  }
893 
894  const Array<OneD,const int>&
895  AssemblyMap::GetBndCondCoeffsToGlobalCoeffsMap()
896  {
897  return m_bndCondCoeffsToGlobalCoeffsMap;
898  }
899 
900 
901  int AssemblyMap::GetNumGlobalDirBndCoeffs() const
902  {
903  return m_numGlobalDirBndCoeffs;
904  }
905 
906 
907  int AssemblyMap::GetNumLocalDirBndCoeffs() const
908  {
909  return m_numLocalDirBndCoeffs;
910  }
911 
912  int AssemblyMap::GetNumLocalBndCoeffs() const
913  {
914  return m_numLocalBndCoeffs;
915  }
916 
917  int AssemblyMap::GetNumGlobalBndCoeffs() const
918  {
919  return m_numGlobalBndCoeffs;
920  }
921 
922  int AssemblyMap::GetNumLocalCoeffs() const
923  {
924  return m_numLocalCoeffs;
925  }
926 
927  int AssemblyMap::GetNumGlobalCoeffs() const
928  {
929  return m_numGlobalCoeffs;
930  }
931 
932  bool AssemblyMap::GetSingularSystem() const
933  {
934  return m_systemSingular;
935  }
936 
937  void AssemblyMap::GlobalToLocalBnd(
938  const NekVector<NekDouble>& global,
939  NekVector<NekDouble>& loc,
940  int offset) const
941  {
942  GlobalToLocalBnd(global.GetPtr(), loc.GetPtr(), offset);
943  }
944 
945 
946  void AssemblyMap::GlobalToLocalBnd(
947  const NekVector<NekDouble>& global,
948  NekVector<NekDouble>& loc) const
949  {
950  GlobalToLocalBnd(global.GetPtr(), loc.GetPtr());
951  }
952 
953 
954  void AssemblyMap::GlobalToLocalBnd(
955  const Array<OneD, const NekDouble>& global,
956  Array<OneD,NekDouble>& loc, int offset) const
957  {
958  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
959  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs-offset,"Global vector is not of correct dimension");
960 
961  // offset input data by length "offset" for Dirichlet boundary conditions.
962  Array<OneD,NekDouble> tmp(m_numGlobalBndCoeffs,0.0);
963  Vmath::Vcopy(m_numGlobalBndCoeffs-offset, global.get(), 1, tmp.get() + offset, 1);
964 
965  if(m_signChange)
966  {
967  Vmath::Gathr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(), tmp.get(), m_localToGlobalBndMap.get(), loc.get());
968  }
969  else
970  {
971  Vmath::Gathr(m_numLocalBndCoeffs, tmp.get(), m_localToGlobalBndMap.get(), loc.get());
972  }
973  }
974 
975 
976  void AssemblyMap::GlobalToLocalBnd(
977  const Array<OneD, const NekDouble>& global,
978  Array<OneD,NekDouble>& loc) const
979  {
980  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
981  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs,"Global vector is not of correct dimension");
982 
983  if(m_signChange)
984  {
985  Vmath::Gathr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(), global.get(), m_localToGlobalBndMap.get(), loc.get());
986  }
987  else
988  {
989  Vmath::Gathr(m_numLocalBndCoeffs, global.get(), m_localToGlobalBndMap.get(), loc.get());
990  }
991  }
992 
993  void AssemblyMap::LocalBndToGlobal(
994  const NekVector<NekDouble>& loc,
995  NekVector<NekDouble>& global,
996  int offset) const
997  {
998  LocalBndToGlobal(loc.GetPtr(), global.GetPtr(), offset);
999  }
1000 
1001 
1002  void AssemblyMap::LocalBndToGlobal(
1003  const Array<OneD, const NekDouble>& loc,
1004  Array<OneD,NekDouble>& global,
1005  int offset) const
1006  {
1007  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
1008  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs-offset,"Global vector is not of correct dimension");
1009 
1010  // offset input data by length "offset" for Dirichlet boundary conditions.
1011  Array<OneD,NekDouble> tmp(m_numGlobalBndCoeffs,0.0);
1012 
1013  if(m_signChange)
1014  {
1015  Vmath::Scatr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(), loc.get(), m_localToGlobalBndMap.get(), tmp.get());
1016  }
1017  else
1018  {
1019  Vmath::Scatr(m_numLocalBndCoeffs, loc.get(), m_localToGlobalBndMap.get(), tmp.get());
1020  }
1021 
1022  UniversalAssembleBnd(tmp);
1023  Vmath::Vcopy(m_numGlobalBndCoeffs-offset, tmp.get()+offset, 1, global.get(), 1);
1024  }
1025 
1026  void AssemblyMap::LocalBndToGlobal(
1027  const NekVector<NekDouble>& loc,
1028  NekVector<NekDouble>& global) const
1029  {
1030  LocalBndToGlobal(loc.GetPtr(), global.GetPtr());
1031  }
1032 
1033 
1034  void AssemblyMap::LocalBndToGlobal(
1035  const Array<OneD, const NekDouble>& loc,
1036  Array<OneD,NekDouble>& global) const
1037  {
1038  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
1039  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs,"Global vector is not of correct dimension");
1040 
1041  if(m_signChange)
1042  {
1043  Vmath::Scatr(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(), loc.get(), m_localToGlobalBndMap.get(), global.get());
1044  }
1045  else
1046  {
1047  Vmath::Scatr(m_numLocalBndCoeffs, loc.get(), m_localToGlobalBndMap.get(), global.get());
1048  }
1049  }
1050 
1051 
1052  void AssemblyMap::AssembleBnd(
1053  const NekVector<NekDouble>& loc,
1054  NekVector<NekDouble>& global, int offset) const
1055  {
1056  AssembleBnd(loc.GetPtr(), global.GetPtr(), offset);
1057  }
1058 
1059 
1060  void AssemblyMap::AssembleBnd(
1061  const NekVector<NekDouble>& loc,
1062  NekVector<NekDouble>& global) const
1063  {
1064  AssembleBnd(loc.GetPtr(), global.GetPtr());
1065  }
1066 
1067 
1068  void AssemblyMap::AssembleBnd(
1069  const Array<OneD,const NekDouble>& loc,
1070  Array<OneD, NekDouble>& global, int offset) const
1071  {
1072  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local array is not of correct dimension");
1073  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs-offset,"Global array is not of correct dimension");
1074  Array<OneD,NekDouble> tmp(m_numGlobalBndCoeffs,0.0);
1075 
1076  if(m_signChange)
1077  {
1078  Vmath::Assmb(m_numLocalBndCoeffs, m_localToGlobalBndSign.get(),loc.get(), m_localToGlobalBndMap.get(), tmp.get());
1079  }
1080  else
1081  {
1082  Vmath::Assmb(m_numLocalBndCoeffs,loc.get(), m_localToGlobalBndMap.get(), tmp.get());
1083  }
1084  UniversalAssembleBnd(tmp);
1085  Vmath::Vcopy(m_numGlobalBndCoeffs-offset, tmp.get() + offset, 1, global.get(), 1);
1086  }
1087 
1088 
1089  void AssemblyMap::AssembleBnd(
1090  const Array<OneD, const NekDouble>& loc,
1091  Array<OneD, NekDouble>& global) const
1092  {
1093  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
1094  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs,"Global vector is not of correct dimension");
1095 
1096  Vmath::Zero(m_numGlobalBndCoeffs, global.get(), 1);
1097 
1098  if(m_signChange)
1099  {
1100  Vmath::Assmb(m_numLocalBndCoeffs,m_localToGlobalBndSign.get(),
1101  loc.get(), m_localToGlobalBndMap.get(), global.get());
1102  }
1103  else
1104  {
1105  Vmath::Assmb(m_numLocalBndCoeffs,loc.get(), m_localToGlobalBndMap.get(), global.get());
1106  }
1107  UniversalAssembleBnd(global);
1108  }
1109 
1110  void AssemblyMap::UniversalAssembleBnd(
1111  Array<OneD, NekDouble>& pGlobal) const
1112  {
1113  ASSERTL1(pGlobal.num_elements() == m_numGlobalBndCoeffs,
1114  "Wrong size.");
1115  Gs::Gather(pGlobal, Gs::gs_add, m_bndGsh);
1116  }
1117 
1118  void AssemblyMap::UniversalAssembleBnd(
1119  NekVector< NekDouble>& pGlobal) const
1120  {
1121  UniversalAssembleBnd(pGlobal.GetPtr());
1122  }
1123 
1124  void AssemblyMap::UniversalAssembleBnd(
1125  Array<OneD, NekDouble>& pGlobal,
1126  int offset) const
1127  {
1128  Array<OneD, NekDouble> tmp(offset);
1129  if (offset > 0) Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
1130  UniversalAssembleBnd(pGlobal);
1131  if (offset > 0) Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
1132  }
1133 
1134  int AssemblyMap::GetBndSystemBandWidth() const
1135  {
1136  return m_bndSystemBandWidth;
1137  }
1138 
1139  int AssemblyMap::GetStaticCondLevel() const
1140  {
1141  return m_staticCondLevel;
1142  }
1143 
1144  int AssemblyMap::GetNumPatches() const
1145  {
1146  return m_numPatches;
1147  }
1148 
1149  const Array<OneD,const unsigned int>&
1150  AssemblyMap::GetNumLocalBndCoeffsPerPatch()
1151  {
1152  return m_numLocalBndCoeffsPerPatch;
1153  }
1154 
1155 
1156  const Array<OneD,const unsigned int>&
1157  AssemblyMap::GetNumLocalIntCoeffsPerPatch()
1158  {
1159  return m_numLocalIntCoeffsPerPatch;
1160  }
1161 
1162  const AssemblyMapSharedPtr
1163  AssemblyMap::GetNextLevelLocalToGlobalMap() const
1164  {
1165  return m_nextLevelLocalToGlobalMap;
1166  }
1167 
1168  const PatchMapSharedPtr&
1169  AssemblyMap::GetPatchMapFromPrevLevel(void)
1170  const
1171  {
1172  return m_patchMapFromPrevLevel;
1173  }
1174 
1175  bool AssemblyMap::AtLastLevel() const
1176  {
1177  return !( m_nextLevelLocalToGlobalMap.get() );
1178  }
1179 
1180 
1181  GlobalSysSolnType AssemblyMap::GetGlobalSysSolnType() const
1182  {
1183  return m_solnType;
1184  }
1185 
1186  PreconditionerType AssemblyMap::GetPreconType() const
1187  {
1188  return m_preconType;
1189  }
1190 
1191  NekDouble AssemblyMap::GetIterativeTolerance() const
1192  {
1193  return m_iterativeTolerance;
1194  }
1195 
1196  int AssemblyMap::GetSuccessiveRHS() const
1197  {
1198  return m_successiveRHS;
1199  }
1200 
1201  void AssemblyMap::GlobalToLocalBndWithoutSign(
1202  const Array<OneD, const NekDouble>& global,
1203  Array<OneD,NekDouble>& loc)
1204  {
1205  ASSERTL1(loc.num_elements() >= m_numLocalBndCoeffs,"Local vector is not of correct dimension");
1206  ASSERTL1(global.num_elements() >= m_numGlobalBndCoeffs,"Global vector is not of correct dimension");
1207 
1208  Vmath::Gathr(m_numLocalBndCoeffs, global.get(), m_localToGlobalBndMap.get(), loc.get());
1209  }
1210 
1211  void AssemblyMap::PrintStats(
1212  std::ostream &out, std::string variable) const
1213  {
1214  LibUtilities::CommSharedPtr vRowComm
1215  = m_session->GetComm()->GetRowComm();
1216  bool isRoot = vRowComm->GetRank() == 0;
1217  int n = vRowComm->GetSize();
1218  int i;
1219 
1220  // Determine number of global degrees of freedom.
1221  int globBndCnt = 0, globDirCnt = 0;
1222 
1223  for (i = 0; i < m_numGlobalBndCoeffs; ++i)
1224  {
1225  if (m_globalToUniversalBndMapUnique[i] > 0)
1226  {
1227  globBndCnt++;
1228 
1229  if (i < m_numGlobalDirBndCoeffs)
1230  {
1231  globDirCnt++;
1232  }
1233  }
1234  }
1235 
1236  int globCnt = m_numGlobalCoeffs - m_numGlobalBndCoeffs + globBndCnt;
1237 
1238  // Calculate maximum valency
1239  Array<OneD, NekDouble> tmpLoc (m_numLocalBndCoeffs, 1.0);
1240  Array<OneD, NekDouble> tmpGlob(m_numGlobalBndCoeffs, 0.0);
1241 
1242  Vmath::Assmb(m_numLocalBndCoeffs, tmpLoc.get(), m_localToGlobalBndMap.get(), tmpGlob.get());
1243  UniversalAssembleBnd(tmpGlob);
1244 
1245  int totGlobDof = globCnt;
1246  int totGlobBndDof = globBndCnt;
1247  int totGlobDirDof = globDirCnt;
1248  int totLocalDof = m_numLocalCoeffs;
1249  int totLocalBndDof = m_numLocalBndCoeffs;
1250  int totLocalDirDof = m_numLocalDirBndCoeffs;
1251 
1252  int meanValence = 0;
1253  int maxValence = 0;
1254  int minValence = 10000000;
1255  for (int i = 0; i < m_numGlobalBndCoeffs; ++i)
1256  {
1257  if (!m_globalToUniversalBndMapUnique[i])
1258  {
1259  continue;
1260  }
1261 
1262  if (tmpGlob[i] > maxValence)
1263  {
1264  maxValence = tmpGlob[i];
1265  }
1266  if (tmpGlob[i] < minValence)
1267  {
1268  minValence = tmpGlob[i];
1269  }
1270  meanValence += tmpGlob[i];
1271  }
1272 
1273  vRowComm->AllReduce(maxValence, LibUtilities::ReduceMax);
1274  vRowComm->AllReduce(minValence, LibUtilities::ReduceMin);
1275  vRowComm->AllReduce(meanValence, LibUtilities::ReduceSum);
1276  vRowComm->AllReduce(totGlobDof, LibUtilities::ReduceSum);
1277  vRowComm->AllReduce(totGlobBndDof, LibUtilities::ReduceSum);
1278  vRowComm->AllReduce(totGlobDirDof, LibUtilities::ReduceSum);
1279  vRowComm->AllReduce(totLocalDof, LibUtilities::ReduceSum);
1280  vRowComm->AllReduce(totLocalBndDof, LibUtilities::ReduceSum);
1281  vRowComm->AllReduce(totLocalDirDof, LibUtilities::ReduceSum);
1282 
1283  meanValence /= totGlobBndDof;
1284 
1285  if (isRoot)
1286  {
1287  out << "Assembly map statistics for field " << variable << ":"
1288  << endl;
1289  out << " - Number of local/global dof : "
1290  << totLocalDof << " " << totGlobDof << endl;
1291  out << " - Number of local/global boundary dof : "
1292  << totLocalBndDof << " " << totGlobBndDof << endl;
1293  out << " - Number of local/global Dirichlet dof : "
1294  << totLocalDirDof << " " << totGlobDirDof << endl;
1295  out << " - dof valency (min/max/mean) : "
1296  << minValence << " " << maxValence << " " << meanValence
1297  << endl;
1298 
1299  if (n > 1)
1300  {
1301  NekDouble mean = m_numLocalCoeffs, mean2 = mean * mean;
1302  NekDouble minval = mean, maxval = mean;
1303  Array<OneD, NekDouble> tmp(1);
1304 
1305  for (i = 1; i < n; ++i)
1306  {
1307  vRowComm->Recv(i, tmp);
1308  mean += tmp[0];
1309  mean2 += tmp[0]*tmp[0];
1310 
1311  if (tmp[0] > maxval)
1312  {
1313  maxval = tmp[0];
1314  }
1315  if (tmp[0] < minval)
1316  {
1317  minval = tmp[0];
1318  }
1319  }
1320 
1321  out << " - Local dof dist. (min/max/mean/dev) : "
1322  << minval << " " << maxval << " " << (mean / n) << " "
1323  << sqrt(mean2/n - mean*mean/n/n) << endl;
1324 
1325  vRowComm->Block();
1326 
1327  mean = minval = maxval = m_numLocalBndCoeffs;
1328  mean2 = mean * mean;
1329 
1330  for (i = 1; i < n; ++i)
1331  {
1332  vRowComm->Recv(i, tmp);
1333  mean += tmp[0];
1334  mean2 += tmp[0]*tmp[0];
1335 
1336  if (tmp[0] > maxval)
1337  {
1338  maxval = tmp[0];
1339  }
1340  if (tmp[0] < minval)
1341  {
1342  minval = tmp[0];
1343  }
1344  }
1345 
1346  out << " - Local bnd dof dist. (min/max/mean/dev) : "
1347  << minval << " " << maxval << " " << (mean / n) << " "
1348  << sqrt(mean2/n - mean*mean/n/n) << endl;
1349  }
1350  }
1351  else
1352  {
1353  Array<OneD, NekDouble> tmp(1);
1354  tmp[0] = m_numLocalCoeffs;
1355  vRowComm->Send(0, tmp);
1356  vRowComm->Block();
1357  tmp[0] = m_numLocalBndCoeffs;
1358  vRowComm->Send(0, tmp);
1359  }
1360  }
1361  } // namespace
1362 } // namespace
Generated on Wed Mar 4 2015 16:33:12 for Nektar++ by   doxygen 1.8.1.2