Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AssemblyMapCG.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AssemblyMapCG.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: C0-continuous Local to Global mapping routines, base class
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/AssemblyMap/AssemblyMapCG.h>
37 #include <MultiRegions/ExpList.h>
38 #include <LocalRegions/Expansion.h>
39 #include <LocalRegions/Expansion2D.h>
40 #include <LocalRegions/Expansion3D.h>
41 
42 
43 #include <boost/config.hpp>
44 #include <boost/graph/adjacency_list.hpp>
45 #include <boost/graph/cuthill_mckee_ordering.hpp>
46 #include <boost/graph/properties.hpp>
47 #include <boost/graph/bandwidth.hpp>
48 
49 namespace Nektar
50 {
51  namespace MultiRegions
52  {
53  /**
54  * @class AssemblyMapCG
55  * Mappings are created for three possible global solution types:
56  * - Direct full matrix
57  * - Direct static condensation
58  * - Direct multi-level static condensation
59  * In the latter case, mappings are created recursively for the
60  * different levels of static condensation.
61  *
62  * These mappings are used by GlobalLinSys to generate the global
63  * system.
64  */
65 
66  /**
67  *
68  */
69  AssemblyMapCG::AssemblyMapCG(
70  const LibUtilities::SessionReaderSharedPtr &pSession,
71  const std::string variable):
72  AssemblyMap(pSession,variable)
73  {
74  pSession->LoadParameter("MaxStaticCondLevel",m_maxStaticCondLevel,100);
75  }
76 
77 
78  /**
79  *
80  */
81  AssemblyMapCG::AssemblyMapCG(
82  const LibUtilities::SessionReaderSharedPtr &pSession,
83  const int numLocalCoeffs,
84  const ExpList &locExp):
85  AssemblyMap(pSession)
86  {
87  ASSERTL0(false,"AssemblyMapCG: you need to instantiate dimension-specific derived class.");
88  }
89 
90 
91 
92  /**
93  *
94  */
95  AssemblyMapCG::~AssemblyMapCG()
96  {
97  }
98 
99 
100 
101  /**
102  * Given faceOrient of a local element to its local face and
103  * perFaceOrient which states the alignment of one periodic
104  * face to the other global face determine a new faceOrient
105  * that takes this local element face to the global/unique
106  * face
107  */
108  StdRegions::Orientation DeterminePeriodicFaceOrient(
109  StdRegions::Orientation faceOrient,
110  StdRegions::Orientation perFaceOrient)
111  {
112 
113  StdRegions::Orientation returnval = faceOrient;
114 
115  if(perFaceOrient != StdRegions::eDir1FwdDir1_Dir2FwdDir2)
116  {
117  int tmp1 = (int)faceOrient - 5;
118  int tmp2 = (int)perFaceOrient - 5;
119 
120  int flipDir1Map [8] = {2,3,0,1,6,7,4,5};
121  int flipDir2Map [8] = {1,0,3,2,5,4,7,6};
122  int transposeMap[8] = {4,5,6,7,0,2,1,3};
123 
124  // Transpose orientation
125  if (tmp2 > 3)
126  {
127  tmp1 = transposeMap[tmp1];
128  }
129 
130  // Reverse orientation in direction 1.
131  if (tmp2 == 2 || tmp2 == 3 || tmp2 == 6 || tmp2 == 7)
132  {
133  tmp1 = flipDir1Map[tmp1];
134  }
135 
136  // Reverse orientation in direction 2
137  if (tmp2 % 2 == 1)
138  {
139  tmp1 = flipDir2Map[tmp1];
140  }
141 
142  returnval = (StdRegions::Orientation)(tmp1+5);
143  }
144  return returnval;
145  }
146 
147 
148  /**
149  * Sets up the global to universal mapping of degrees of freedom across
150  * processors.
151  */
152  void AssemblyMapCG::SetUpUniversalC0ContMap(
153  const ExpList &locExp,
154  const PeriodicMap &perVerts,
155  const PeriodicMap &perEdges,
156  const PeriodicMap &perFaces)
157  {
158  LocalRegions::ExpansionSharedPtr locExpansion;
159  int nDim = 0;
160  int nVert = 0;
161  int nEdge = 0;
162  int nFace = 0;
163  int maxEdgeDof = 0;
164  int maxFaceDof = 0;
165  int maxIntDof = 0;
166  int dof = 0;
167  int cnt;
168  int i,j,k;
169  int meshVertId;
170  int meshEdgeId;
171  int meshFaceId;
172  int elementId;
173  int vGlobalId;
174  int maxBndGlobalId = 0;
175  StdRegions::Orientation edgeOrient;
176  StdRegions::Orientation faceOrient;
177  Array<OneD, unsigned int> edgeInteriorMap;
178  Array<OneD, int> edgeInteriorSign;
179  Array<OneD, unsigned int> faceInteriorMap;
180  Array<OneD, int> faceInteriorSign;
181  Array<OneD, unsigned int> interiorMap;
182  PeriodicMap::const_iterator pIt;
183 
184  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
185  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
186 
187  m_globalToUniversalMap = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
188  m_globalToUniversalMapUnique = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
189  m_globalToUniversalBndMap = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
190  m_globalToUniversalBndMapUnique = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
191 
192  // Loop over all the elements in the domain to gather mesh data
193  for(i = 0; i < locExpVector.size(); ++i)
194  {
195  locExpansion = locExpVector[i];
196  nVert += locExpansion->GetNverts();
197  nEdge += locExpansion->GetNedges();
198  nFace += locExpansion->GetNfaces();
199  // Loop over all edges (and vertices) of element i
200  for(j = 0; j < locExpansion->GetNedges(); ++j)
201  {
202  dof = locExpansion->GetEdgeNcoeffs(j)-2;
203  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
204  }
205  for(j = 0; j < locExpansion->GetNfaces(); ++j)
206  {
207  dof = locExpansion->GetFaceIntNcoeffs(j);
208  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
209  }
210  locExpansion->GetInteriorMap(interiorMap);
211  dof = interiorMap.num_elements();
212  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
213  }
214 
215  // Tell other processes about how many dof we have
216  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
217  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
218  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
219  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
220  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
221  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
222 
223  // Assemble global to universal mapping for this process
224  for(i = 0; i < locExpVector.size(); ++i)
225  {
226  locExpansion = locExpVector[i];
227  nDim = locExpansion->GetShapeDimension();
228  cnt = locExp.GetCoeff_Offset(i);
229 
230  // Loop over all vertices of element i
231  for(j = 0; j < locExpansion->GetNverts(); ++j)
232  {
233  meshVertId = locExpansion->GetGeom()->GetVid(j);
234  vGlobalId = m_localToGlobalMap[cnt+locExpansion->GetVertexMap(j)];
235 
236  pIt = perVerts.find(meshVertId);
237  if (pIt != perVerts.end())
238  {
239  for (k = 0; k < pIt->second.size(); ++k)
240  {
241  meshVertId = min(meshVertId, pIt->second[k].id);
242  }
243  }
244 
245  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
246  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
247  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
248  }
249 
250  // Loop over all edges of element i
251  for(j = 0; j < locExpansion->GetNedges(); ++j)
252  {
253  meshEdgeId = locExpansion->GetGeom()->GetEid(j);
254  pIt = perEdges.find(meshEdgeId);
255  if (pIt != perEdges.end())
256  {
257  for (k = 0; k < pIt->second.size(); ++k)
258  {
259  meshEdgeId = min(meshEdgeId, pIt->second[k].id);
260  }
261  }
262 
263  edgeOrient = locExpansion->GetGeom()->GetEorient(j);
264  locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
265  dof = locExpansion->GetEdgeNcoeffs(j)-2;
266 
267  // Set the global DOF's for the interior modes of edge j
268  for(k = 0; k < dof; ++k)
269  {
270  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
271  m_globalToUniversalMap[vGlobalId]
272  = nVert + meshEdgeId * maxEdgeDof + k + 1;
273  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
274  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
275  }
276  }
277 
278  // Loop over all faces of element i
279  for(j = 0; j < locExpansion->GetNfaces(); ++j)
280  {
281  faceOrient = locExpansion->as<LocalRegions::Expansion3D>()
282  ->GetGeom3D()->GetFaceOrient(j);
283 
284  meshFaceId = locExpansion->GetGeom()->GetFid(j);
285 
286  pIt = perFaces.find(meshFaceId);
287  if (pIt != perFaces.end())
288  {
289  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
290  {
291  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
292  }
293  meshFaceId = min(meshFaceId, pIt->second[0].id);
294  }
295 
296 
297  locExpansion->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
298  dof = locExpansion->GetFaceIntNcoeffs(j);
299 
300 
301  for(k = 0; k < dof; ++k)
302  {
303  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
304  m_globalToUniversalMap[vGlobalId]
305  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
306  + k + 1;
307  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
308 
309  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
310  }
311  }
312 
313  // Add interior DOFs to complete universal numbering
314  locExpansion->GetInteriorMap(interiorMap);
315  dof = interiorMap.num_elements();
316  elementId = (locExpansion->GetGeom())->GetGlobalID();
317  for (k = 0; k < dof; ++k)
318  {
319  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
320  m_globalToUniversalMap[vGlobalId]
321  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
322  }
323  }
324 
325  // Set up the GSLib universal assemble mapping
326  // Internal DOF do not participate in any data
327  // exchange, so we keep these set to the special GSLib id=0 so
328  // they are ignored.
329  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
330  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
331  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
332  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
333  {
334  tmp[i] = m_globalToUniversalMap[i];
335  }
336 
337  m_gsh = Gs::Init(tmp, vCommRow);
338  m_bndGsh = Gs::Init(tmp2, vCommRow);
339  Gs::Unique(tmp, vCommRow);
340  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
341  {
342  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
343  }
344  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
345  {
346  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
347  }
348  }
349 
350  /**
351  * @brief Construct an AssemblyMapCG object which corresponds to the
352  * linear space of the current object.
353  *
354  * This function is used to create a linear-space assembly map, which is
355  * then used in the linear space preconditioner in the conjugate
356  * gradient solve.
357  */
358  AssemblyMapSharedPtr AssemblyMapCG::v_LinearSpaceMap(
359  const ExpList &locexp, GlobalSysSolnType solnType)
360  {
361  AssemblyMapCGSharedPtr returnval;
362 
363  int i, j;
364  int nverts = 0;
365  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
366  = locexp.GetExp();
367  int nelmts = exp->size();
368 
369  // Get Default Map and turn off any searched values.
370  returnval = MemoryManager<AssemblyMapCG>
371  ::AllocateSharedPtr(m_session);
372  returnval->m_solnType = solnType;
373  returnval->m_preconType = eNull;
374  returnval->m_maxStaticCondLevel = 0;
375  returnval->m_signChange = false;
376  returnval->m_comm = m_comm;
377 
378  // Count the number of vertices
379  for (i = 0; i < nelmts; ++i)
380  {
381  nverts += (*exp)[i]->GetNverts();
382  }
383 
384  returnval->m_numLocalCoeffs = nverts;
385  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
386 
387  // Store original global ids in this map
388  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
389 
390  int cnt = 0;
391  int cnt1 = 0;
392  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
393 
394  // Set up local to global map;
395  for (i = 0; i < nelmts; ++i)
396  {
397  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
398  {
399  returnval->m_localToGlobalMap[cnt] =
400  returnval->m_localToGlobalBndMap[cnt] =
401  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
402  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
403 
404 #if 1
405  // Set up numLocalDirBndCoeffs
406  if ((returnval->m_localToGlobalMap[cnt]) <
407  m_numGlobalDirBndCoeffs)
408  {
409  returnval->m_numLocalDirBndCoeffs++;
410  }
411 #endif
412  cnt++;
413  }
414  cnt1 += (*exp)[i]->GetNcoeffs();
415  }
416 
417  cnt = 0;
418  // Reset global numbering and count number of dofs
419  for (i = 0; i < m_numGlobalCoeffs; ++i)
420  {
421  if (GlobCoeffs[i] != -1)
422  {
423  GlobCoeffs[i] = cnt++;
424  }
425  }
426 
427  // Set up number of globalCoeffs;
428  returnval->m_numGlobalCoeffs = cnt;
429 
430  // Set up number of global Dirichlet boundary coefficients
431  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
432  {
433  if (GlobCoeffs[i] != -1)
434  {
435  returnval->m_numGlobalDirBndCoeffs++;
436  }
437  }
438 
439  // Set up global to universal map
440  if (m_globalToUniversalMap.num_elements())
441  {
442  LibUtilities::CommSharedPtr vCommRow
443  = m_session->GetComm()->GetRowComm();
444  int nglocoeffs = returnval->m_numGlobalCoeffs;
445  returnval->m_globalToUniversalMap
446  = Array<OneD, int> (nglocoeffs);
447  returnval->m_globalToUniversalMapUnique
448  = Array<OneD, int> (nglocoeffs);
449 
450  // Reset local to global map and setup universal map
451  for (i = 0; i < nverts; ++i)
452  {
453  cnt = returnval->m_localToGlobalMap[i];
454  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
455 
456  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
457  m_globalToUniversalMap[cnt];
458  }
459 
460  Nektar::Array<OneD, long> tmp(nglocoeffs);
461  Vmath::Zero(nglocoeffs, tmp, 1);
462  for (unsigned int i = 0; i < nglocoeffs; ++i)
463  {
464  tmp[i] = returnval->m_globalToUniversalMap[i];
465  }
466  returnval->m_gsh = Gs::Init(tmp, vCommRow);
467  Gs::Unique(tmp, vCommRow);
468  for (unsigned int i = 0; i < nglocoeffs; ++i)
469  {
470  returnval->m_globalToUniversalMapUnique[i]
471  = (tmp[i] >= 0 ? 1 : 0);
472  }
473  }
474  else // not sure this option is ever needed.
475  {
476  for (i = 0; i < nverts; ++i)
477  {
478  cnt = returnval->m_localToGlobalMap[i];
479  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
480  }
481  }
482 
483  return returnval;
484  }
485 
486  /**
487  * The bandwidth calculated here corresponds to what is referred to as
488  * half-bandwidth. If the elements of the matrix are designated as
489  * a_ij, it corresponds to the maximum value of |i-j| for non-zero
490  * a_ij. As a result, the value also corresponds to the number of
491  * sub- or super-diagonals.
492  *
493  * The bandwith can be calculated elementally as it corresponds to the
494  * maximal elemental bandwith (i.e. the maximal difference in global
495  * DOF index for every element).
496  *
497  * We caluclate here the bandwith of the full global system.
498  */
499  void AssemblyMapCG::CalculateFullSystemBandWidth()
500  {
501  int i,j;
502  int cnt = 0;
503  int locSize;
504  int maxId;
505  int minId;
506  int bwidth = -1;
507  for(i = 0; i < m_numPatches; ++i)
508  {
509  locSize = m_numLocalBndCoeffsPerPatch[i]+m_numLocalIntCoeffsPerPatch[i];
510  maxId = -1;
511  minId = m_numLocalCoeffs+1;
512  for(j = 0; j < locSize; j++)
513  {
514  if(m_localToGlobalMap[cnt+j] >= m_numGlobalDirBndCoeffs)
515  {
516  if(m_localToGlobalMap[cnt+j] > maxId)
517  {
518  maxId = m_localToGlobalMap[cnt+j];
519  }
520 
521  if(m_localToGlobalMap[cnt+j] < minId)
522  {
523  minId = m_localToGlobalMap[cnt+j];
524  }
525  }
526  }
527  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
528 
529  cnt+=locSize;
530  }
531 
532  m_fullSystemBandWidth = bwidth;
533  }
534 
535 
536  int AssemblyMapCG::v_GetLocalToGlobalMap(const int i) const
537  {
538  return m_localToGlobalMap[i];
539  }
540 
541  int AssemblyMapCG::v_GetGlobalToUniversalMap(const int i) const
542  {
543  return m_globalToUniversalMap[i];
544  }
545 
546  int AssemblyMapCG::v_GetGlobalToUniversalMapUnique(const int i) const
547  {
548  return m_globalToUniversalMapUnique[i];
549  }
550 
551  const Array<OneD,const int>&
552  AssemblyMapCG::v_GetLocalToGlobalMap(void)
553  {
554  return m_localToGlobalMap;
555  }
556 
557  const Array<OneD,const int>&
558  AssemblyMapCG::v_GetGlobalToUniversalMap(void)
559  {
560  return m_globalToUniversalMap;
561  }
562 
563  const Array<OneD,const int>&
564  AssemblyMapCG::v_GetGlobalToUniversalMapUnique(void)
565  {
566  return m_globalToUniversalMapUnique;
567  }
568 
569  NekDouble AssemblyMapCG::v_GetLocalToGlobalSign(
570  const int i) const
571  {
572  if(m_signChange)
573  {
574  return m_localToGlobalSign[i];
575  }
576  else
577  {
578  return 1.0;
579  }
580  }
581 
582  const Array<OneD, NekDouble>& AssemblyMapCG::v_GetLocalToGlobalSign() const
583  {
584  return m_localToGlobalSign;
585  }
586 
587  void AssemblyMapCG::v_LocalToGlobal(
588  const Array<OneD, const NekDouble>& loc,
589  Array<OneD, NekDouble>& global) const
590  {
591  Array<OneD, const NekDouble> local;
592  if(global.data() == loc.data())
593  {
594  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
595  }
596  else
597  {
598  local = loc; // create reference
599  }
600 
601 
602  if(m_signChange)
603  {
604  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
605  }
606  else
607  {
608  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
609  }
610 
611  // ensure all values are unique by calling a max
612  Gs::Gather(global, Gs::gs_max, m_gsh);
613  }
614 
615  void AssemblyMapCG::v_LocalToGlobal(
616  const NekVector<NekDouble>& loc,
617  NekVector< NekDouble>& global) const
618  {
619  LocalToGlobal(loc.GetPtr(),global.GetPtr());
620  }
621 
622  void AssemblyMapCG::v_GlobalToLocal(
623  const Array<OneD, const NekDouble>& global,
624  Array<OneD, NekDouble>& loc) const
625  {
626  Array<OneD, const NekDouble> glo;
627  if(global.data() == loc.data())
628  {
629  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
630  }
631  else
632  {
633  glo = global; // create reference
634  }
635 
636 
637  if(m_signChange)
638  {
639  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
640  }
641  else
642  {
643  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
644  }
645  }
646 
647  void AssemblyMapCG::v_GlobalToLocal(
648  const NekVector<NekDouble>& global,
649  NekVector< NekDouble>& loc) const
650  {
651  GlobalToLocal(global.GetPtr(),loc.GetPtr());
652  }
653 
654  void AssemblyMapCG::v_Assemble(
655  const Array<OneD, const NekDouble> &loc,
656  Array<OneD, NekDouble> &global) const
657  {
658  Array<OneD, const NekDouble> local;
659  if(global.data() == loc.data())
660  {
661  local = Array<OneD, NekDouble>(local.num_elements(),local.data());
662  }
663  else
664  {
665  local = loc; // create reference
666  }
667  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
668 
669  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
670 
671  if(m_signChange)
672  {
673  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
674  }
675  else
676  {
677  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
678  }
679  UniversalAssemble(global);
680  }
681 
682  void AssemblyMapCG::v_Assemble(
683  const NekVector<NekDouble>& loc,
684  NekVector< NekDouble>& global) const
685  {
686  Assemble(loc.GetPtr(),global.GetPtr());
687  }
688 
689  void AssemblyMapCG::v_UniversalAssemble(
690  Array<OneD, NekDouble>& pGlobal) const
691  {
692  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
693  }
694 
695  void AssemblyMapCG::v_UniversalAssemble(
696  NekVector< NekDouble>& pGlobal) const
697  {
698  UniversalAssemble(pGlobal.GetPtr());
699  }
700 
701  void AssemblyMapCG::v_UniversalAssemble(
702  Array<OneD, NekDouble>& pGlobal,
703  int offset) const
704  {
705  Array<OneD, NekDouble> tmp(offset);
706  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
707  UniversalAssemble(pGlobal);
708  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
709  }
710 
711  int AssemblyMapCG::v_GetFullSystemBandWidth() const
712  {
713  return m_fullSystemBandWidth;
714  }
715 
716  int AssemblyMapCG::v_GetNumNonDirVertexModes() const
717  {
718  return m_numNonDirVertexModes;
719  }
720 
721  int AssemblyMapCG::v_GetNumNonDirEdgeModes() const
722  {
723  return m_numNonDirEdgeModes;
724  }
725 
726  int AssemblyMapCG::v_GetNumNonDirFaceModes() const
727  {
728  return m_numNonDirFaceModes;
729  }
730 
731  int AssemblyMapCG::v_GetNumDirEdges() const
732  {
733  return m_numDirEdges;
734  }
735 
736  int AssemblyMapCG::v_GetNumDirFaces() const
737  {
738  return m_numDirFaces;
739  }
740 
741  int AssemblyMapCG::v_GetNumNonDirEdges() const
742  {
743  return m_numNonDirEdges;
744  }
745 
746  int AssemblyMapCG::v_GetNumNonDirFaces() const
747  {
748  return m_numNonDirFaces;
749  }
750 
751  const Array<OneD, const int>& AssemblyMapCG::v_GetExtraDirEdges()
752  {
753  return m_extraDirEdges;
754  }
755  } // namespace
756 } // namespace
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2