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