Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AssemblyMapCG2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AssemblyMapCG2D.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 assembly mappings specific to 2D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/MultiRegions.hpp>
37 #include <MultiRegions/AssemblyMap/AssemblyMapCG2D.h>
38 #include <MultiRegions/ExpList.h>
39 #include <LocalRegions/SegExp.h>
40 #include <LocalRegions/Expansion2D.h>
41 
42 #include <boost/config.hpp>
43 #include <boost/graph/adjacency_list.hpp>
44 #include <boost/graph/cuthill_mckee_ordering.hpp>
45 #include <boost/graph/properties.hpp>
46 #include <boost/graph/bandwidth.hpp>
47 
48 #include <iomanip>
49 
50 namespace Nektar
51 {
52  namespace MultiRegions
53  {
54  /**
55  * @class AssemblyMapCG2D
56  * Mappings are created for three possible global solution types:
57  * - Direct full matrix
58  * - Direct static condensation
59  * - Direct multi-level static condensation
60  * In the latter case, mappings are created recursively for the
61  * different levels of static condensation.
62  *
63  * These mappings are used by GlobalLinSys to generate the global
64  * system.
65  */
66 
67  /**
68  *
69  */
70  AssemblyMapCG2D::AssemblyMapCG2D(
71  const LibUtilities::SessionReaderSharedPtr &pSession,
72  const std::string variable):
73  AssemblyMapCG(pSession,variable)
74  {
75  }
76 
77 
78 
79 
80  /**
81  *
82  */
83  AssemblyMapCG2D::AssemblyMapCG2D(
84  const LibUtilities::SessionReaderSharedPtr &pSession,
85  const int numLocalCoeffs,
86  const ExpList &locExp,
87  const Array<OneD, const ExpListSharedPtr> &bndCondExp,
88  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
89  &bndConditions,
90  const PeriodicMap& periodicVertsId,
91  const PeriodicMap& periodicEdgesId,
92  const bool checkIfSystemSingular,
93  const std::string variable) :
94  AssemblyMapCG(pSession,variable)
95  {
96  SetUp2DExpansionC0ContMap(numLocalCoeffs,
97  locExp,
98  bndCondExp,
99  bndConditions,
100  periodicVertsId,
101  periodicEdgesId,
102  checkIfSystemSingular);
103 
104  CalculateBndSystemBandWidth();
105  CalculateFullSystemBandWidth();
106  }
107 
108 
109  /**
110  *
111  */
112  AssemblyMapCG2D::AssemblyMapCG2D(
113  const LibUtilities::SessionReaderSharedPtr &pSession,
114  const int numLocalCoeffs,
115  const ExpList &locExp,
116  const std::string variable):
117  AssemblyMapCG(pSession,variable)
118  {
119  SetUp2DExpansionC0ContMap(numLocalCoeffs, locExp);
120  CalculateBndSystemBandWidth();
121  CalculateFullSystemBandWidth();
122  }
123 
124 
125  /**
126  *
127  */
128  AssemblyMapCG2D::~AssemblyMapCG2D()
129  {
130  }
131 
132 
133  /**
134  * Construction of the local->global map is achieved in
135  * several stages. A mesh vertex and mesh edge renumbering is
136  * constructed in #vertReorderedGraphVertId and
137  * #edgeReorderedGraphVertId through a call ot
138  * #SetUp2DGraphC0ContMap.
139  *
140  * The local numbering is then deduced in the following steps:
141  */
142  void AssemblyMapCG2D::SetUp2DExpansionC0ContMap(
143  const int numLocalCoeffs,
144  const ExpList &locExp,
145  const Array<OneD, const ExpListSharedPtr> &bndCondExp,
146  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
147  &bndConditions,
148  const PeriodicMap& periodicVertsId,
149  const PeriodicMap& periodicEdgesId,
150  const bool checkIfSystemSingular)
151  {
152  int i,j,k;
153  int cnt = 0,offset=0;
154  int meshVertId;
155  int meshEdgeId;
156  int bndEdgeCnt;
157  int globalId, nGraphVerts;
158  int nEdgeCoeffs;
159  int nEdgeInteriorCoeffs;
160  int firstNonDirGraphVertId;
161  int nLocBndCondDofs = 0;
162  int nLocDirBndCondDofs = 0;
163  int nExtraDirichlet = 0;
164  LocalRegions::Expansion2DSharedPtr locExpansion;
165  LocalRegions::SegExpSharedPtr bndSegExp;
166  LibUtilities::BasisType bType;
167  StdRegions::Orientation edgeOrient;
168  Array<OneD, unsigned int> edgeInteriorMap;
169  Array<OneD, int> edgeInteriorSign;
170  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
171  m_signChange = false;
172  m_systemSingular = false;
173  Array<OneD, map<int,int> > ReorderedGraphVertId(2);
174  Array<OneD, map<int,int> > Dofs(2);
175  BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
176  set<int> extraDirVerts;
177  LocalRegions::Expansion2DSharedPtr exp2d;
178 
179  for(i = 0; i < locExpVector.size(); ++i)
180  {
181  exp2d = locExpVector[i]->as<LocalRegions::Expansion2D>();
182  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
183  {
184  // Vert Dofs
185  Dofs[0][(exp2d->GetGeom2D())->GetVid(j)] = 1;
186  // Edge Dofs
187  Dofs[1][(exp2d->GetGeom2D())->GetEid(j)] =
188  exp2d->GetEdgeNcoeffs(j)-2;
189  }
190  }
191 
192  /**
193  * STEP 1: Wrap boundary conditions in an array (since
194  * routine is set up for multiple fields) and call the
195  * graph re-ordering subroutine to obtain the reordered
196  * values
197  */
198  Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> > bndConditionsVec(1,bndConditions);
199  nGraphVerts = SetUp2DGraphC0ContMap(locExp,
200  bndCondExp,bndConditionsVec,
201  periodicVertsId,
202  periodicEdgesId,
203  Dofs,
204  ReorderedGraphVertId,
205  firstNonDirGraphVertId,
206  nExtraDirichlet,
207  bottomUpGraph,
208  extraDirVerts,
209  checkIfSystemSingular);
210 
211  /**
212  * STEP 2: Count out the number of Dirichlet vertices and
213  * edges first
214  */
215  for(i = 0; i < bndCondExp.num_elements(); i++)
216  {
217  for(j = 0; j < bndCondExp[i]->GetExpSize(); j++)
218  {
219  bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
220  if(bndConditions[i]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
221  {
222  nLocDirBndCondDofs += bndSegExp->GetNcoeffs();
223  }
224  if( bndCondExp[i]->GetExp(j)==bndCondExp[bndCondExp.num_elements()-1]->GetExp(bndCondExp[bndCondExp.num_elements()-1]->GetExpSize()-1)
225  && i==(bndCondExp.num_elements()-1)
226  && nLocDirBndCondDofs ==0
227  && checkIfSystemSingular==true)
228  {
229  nLocDirBndCondDofs =1;
230  m_systemSingular=true;
231  }
232 
233  nLocBndCondDofs += bndSegExp->GetNcoeffs();
234  }
235  }
236  m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet;
237 
238  /**
239  * STEP 3: Set up an array which contains the offset information of
240  * the different graph vertices.
241  *
242  * This basically means to identify to how many global degrees of
243  * freedom the individual graph vertices correspond. Obviously,
244  * the graph vertices corresponding to the mesh-vertices account
245  * for a single global DOF. However, the graph vertices
246  * corresponding to the element edges correspond to N-2 global DOF
247  * where N is equal to the number of boundary modes on this edge.
248  */
249  Array<OneD, int> graphVertOffset(ReorderedGraphVertId[0].size()+
250  ReorderedGraphVertId[1].size()+1);
251  graphVertOffset[0] = 0;
252  m_numLocalBndCoeffs = 0;
253  for(i = 0; i < locExpVector.size(); ++i)
254  {
255  locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]->as<LocalRegions::Expansion2D>();
256  for(j = 0; j < locExpansion->GetNedges(); ++j)
257  {
258  nEdgeCoeffs = locExpansion->GetEdgeNcoeffs(j);
259  meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
260  meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
261  graphVertOffset[ReorderedGraphVertId[0][meshVertId]+1] = 1;
262  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]+1] = nEdgeCoeffs-2;
263 
264  bType = locExpansion->GetEdgeBasisType(j);
265  // need a sign vector for modal expansions if nEdgeCoeffs >=4
266  if( (nEdgeCoeffs >= 4)&&
267  ( (bType == LibUtilities::eModified_A)||
268  (bType == LibUtilities::eModified_B) ) )
269  {
270  m_signChange = true;
271  }
272  }
273  m_numLocalBndCoeffs += locExpVector[i]->NumBndryCoeffs();
274  }
275 
276  for(i = 1; i < graphVertOffset.num_elements(); i++)
277  {
278  graphVertOffset[i] += graphVertOffset[i-1];
279  }
280 
281  // Allocate the proper amount of space for the class-data and fill
282  // information that is already known
283  m_numLocalCoeffs = numLocalCoeffs;
284  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
285  m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
286  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
287  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int>(nLocBndCondDofs,-1);
288 
289  // If required, set up the sign-vector
290  if(m_signChange)
291  {
292  m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
293  m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
294  m_bndCondCoeffsToGlobalCoeffsSign = Array<OneD, NekDouble>(nLocBndCondDofs,1.0);
295  }
296  else
297  {
298  m_localToGlobalSign = NullNekDouble1DArray;
299  m_localToGlobalBndSign = NullNekDouble1DArray;
300  m_bndCondCoeffsToGlobalCoeffsSign = NullNekDouble1DArray;
301  }
302 
303  // Set up information for multi-level static condensation.
304  m_staticCondLevel = 0;
305  m_numPatches = locExpVector.size();
306  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
307  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
308  for(i = 0; i < m_numPatches; ++i)
309  {
310  int elmtid = locExp.GetOffset_Elmt_Id(i);
311  locExpansion = locExpVector[elmtid]->as<LocalRegions::Expansion2D>();
312  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
313  locExpVector[elmtid]->NumBndryCoeffs();
314  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
315  locExpVector[elmtid]->GetNcoeffs() -
316  locExpVector[elmtid]->NumBndryCoeffs();
317  }
318 
319  /**
320  * STEP 4: Now, all ingredients are ready to set up the actual
321  * local to global mapping.
322  *
323  * The remainder of the map consists of the element-interior
324  * degrees of freedom. This leads to the block-diagonal submatrix
325  * as each element-interior mode is globally orthogonal to modes
326  * in all other elements.
327  */
328  cnt = 0;
329  // Loop over all the elements in the domain
330  for(i = 0; i < locExpVector.size(); ++i)
331  {
332  locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
333 
334  cnt = locExp.GetCoeff_Offset(i);
335 
336  // Loop over all edges (and vertices) of element i
337  for(j = 0; j < locExpansion->GetNedges(); ++j)
338  {
339  PeriodicMap::const_iterator it;
340 
341  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
342  edgeOrient = (locExpansion->GetGeom2D())->GetEorient(j);
343  meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
344  meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
345 
346  /*
347  * Where the edge orientations of periodic edges matches,
348  * we reverse the vertices of each edge in
349  * DisContField2D::GetPeriodicEdges. We must therefore
350  * reverse the orientation of precisely one of the two
351  * edges so that the sign array is correctly populated.
352  */
353  it = periodicEdgesId.find(meshEdgeId);
354  if (it != periodicEdgesId.end() &&
355  it->second[0].orient == StdRegions::eBackwards &&
356  meshEdgeId == min(meshEdgeId, it->second[0].id))
357  {
358  edgeOrient = edgeOrient == StdRegions::eForwards ?
359  StdRegions::eBackwards :
360  StdRegions::eForwards;
361  }
362 
363  locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
364  // Set the global DOF for vertex j of element i
365  m_localToGlobalMap[cnt+locExpansion->GetVertexMap(j)] =
366  graphVertOffset[ReorderedGraphVertId[0][meshVertId]];
367 
368  // Set the global DOF's for the interior modes of edge j
369  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
370  {
371  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
372  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]]+k;
373  }
374 
375  // Fill the sign vector if required
376  if(m_signChange)
377  {
378  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
379  {
380  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
381  }
382  }
383  }
384 
385  cnt += locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs();
386  }
387 
388 
389  // Set up the mapping for the boundary conditions
390  offset = cnt = 0;
391  for(i = 0; i < bndCondExp.num_elements(); i++)
392  {
393  set<int> foundExtraVerts;
394  for(j = 0; j < bndCondExp[i]->GetExpSize(); j++)
395  {
396  bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
397 
398  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
399  for(k = 0; k < 2; k++)
400  {
401  meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
402  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]];
403 
404  set<int>::iterator iter = extraDirVerts.find(meshVertId);
405  if (iter != extraDirVerts.end() &&
406  foundExtraVerts.count(meshVertId) == 0)
407  {
408  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
409  bndSegExp->GetVertexMap(k);
410  int gid = graphVertOffset[
411  ReorderedGraphVertId[0][meshVertId]];
412  ExtraDirDof t(loc, gid, 1.0);
413  m_extraDirDofs[i].push_back(t);
414  foundExtraVerts.insert(meshVertId);
415  }
416  }
417 
418  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
419  bndEdgeCnt = 0;
420  for(k = 0; k < bndSegExp->GetNcoeffs(); k++)
421  {
422  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
423  {
424  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
425  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]]+bndEdgeCnt;
426  bndEdgeCnt++;
427  }
428  }
429  }
430  offset += bndCondExp[i]->GetNcoeffs();
431  }
432 
433  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
434  m_numGlobalBndCoeffs = globalId;
435 
436 
437  /**
438  * STEP 5: The boundary condition mapping is generated from the
439  * same vertex renumbering.
440  */
441  cnt=0;
442  for(i = 0; i < m_numLocalCoeffs; ++i)
443  {
444  if(m_localToGlobalMap[i] == -1)
445  {
446  m_localToGlobalMap[i] = globalId++;
447  }
448  else
449  {
450  if(m_signChange)
451  {
452  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
453  }
454  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
455  }
456  }
457  m_numGlobalCoeffs = globalId;
458 
459  SetUpUniversalC0ContMap(locExp, periodicVertsId, periodicEdgesId);
460 
461  // Since we can now have multiple entries to m_extraDirDofs due to
462  // periodic boundary conditions we make a call to work out the
463  // multiplicity of all entries and invert
464  Array<OneD, NekDouble> valence(m_numGlobalBndCoeffs,0.0);
465 
466  // Fill in Dirichlet coefficients that are to be sent to other
467  // processors with a value of 1
468  map<int, vector<ExtraDirDof> >::iterator Tit;
469 
470  // Generate valence for extraDirDofs
471  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
472  {
473  for (i = 0; i < Tit->second.size(); ++i)
474  {
475  valence[Tit->second[i].get<1>()] = 1.0;
476  }
477  }
478 
479  UniversalAssembleBnd(valence);
480 
481  // Set third argument of tuple to inverse of valence.
482  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
483  {
484  for (i = 0; i < Tit->second.size(); ++i)
485  {
486  boost::get<2>(Tit->second.at(i)) = 1.0/valence[Tit->second.at(i).get<1>()];
487  }
488  }
489 
490  // Set up the local to global map for the next level when using
491  // multi-level static condensation
492  if ((m_solnType == eDirectMultiLevelStaticCond ||
493  m_solnType == eIterativeMultiLevelStaticCond ||
494  m_solnType == eXxtMultiLevelStaticCond) && nGraphVerts)
495  {
496  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1) &&
497  m_staticCondLevel < m_maxStaticCondLevel)
498  {
499  Array<OneD, int> vwgts_perm(
500  Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
501  for(i = 0; i < locExpVector.size(); ++i)
502  {
503  int eid = locExp.GetOffset_Elmt_Id(i);
504  locExpansion = locExpVector[eid]->as<LocalRegions::Expansion2D>();
505  for(j = 0; j < locExpansion->GetNverts(); ++j)
506  {
507  meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
508  meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
509 
510  if(ReorderedGraphVertId[0][meshVertId] >=
511  firstNonDirGraphVertId)
512  {
513  vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
514  firstNonDirGraphVertId] =
515  Dofs[0][meshVertId];
516  }
517 
518  if(ReorderedGraphVertId[1][meshEdgeId] >=
519  firstNonDirGraphVertId)
520  {
521  vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
522  firstNonDirGraphVertId] =
523  Dofs[1][meshEdgeId];
524  }
525  }
526  }
527 
528  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
529 
530  m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
531  AllocateSharedPtr(this,bottomUpGraph);
532  }
533  }
534 
535  m_hash = boost::hash_range(
536  m_localToGlobalMap.begin(), m_localToGlobalMap.end());
537 
538  // Add up hash values if parallel
539  int hash = m_hash;
540  m_comm->GetRowComm()->AllReduce(hash,
541  LibUtilities::ReduceSum);
542  m_hash = hash;
543  }
544 
545  /**
546  * The only unique identifiers of the vertices and edges of the mesh are
547  * the vertex id and the mesh id (stored in their corresponding Geometry
548  * object). However, setting up a global numbering based on these id's
549  * will not lead to a suitable or optimal numbering. Mainly because:
550  * - we want the Dirichlet DOF's to be listed first
551  * - we want an optimal global numbering of the remaining DOF's
552  * (strategy still need to be defined but can for example be: minimum
553  * bandwith or minimum fill-in of the resulting global system matrix)
554  *
555  * The vertices and egdes therefore need to be rearranged which is
556  * perofrmed in in the following way: The vertices and edges of the mesh
557  * are considered as vertices of a graph (in a computer science
558  * terminology, equivalently, they can also be considered as boundary
559  * degrees of freedom, whereby all boundary modes of a single edge are
560  * considered as a single DOF). We then will use different algorithms to
561  * reorder the graph-vertices.
562  *
563  * In the following we use a boost graph object to store this graph the
564  * first template parameter (=OutEdgeList) is chosen to be of type
565  * std::set. Similarly we also use a std::set to hold the adjacency
566  * information. A similar edge might exist multiple times and so to
567  * prevent the definition of parallel edges, we use std::set
568  * (=boost::setS) rather than std::vector (=boost::vecS).
569  *
570  * Two different containers are used to store the graph vertex id's of
571  * the different mesh vertices and edges. They are implemented as a STL
572  * map such that the graph vertex id can later be retrieved by the
573  * unique mesh vertex or edge id's which serve as the key of the map.
574  *
575  * Therefore, the algorithm proceeds as follows:
576  */
577  int AssemblyMapCG2D::SetUp2DGraphC0ContMap(
578  const ExpList &locExp,
579  const Array<OneD, const ExpListSharedPtr> &bndCondExp,
580  const Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> > &bndConditions,
581  const PeriodicMap& periodicVerts,
582  const PeriodicMap& periodicEdges,
583  Array<OneD, map<int,int> > &Dofs,
584  Array<OneD, map<int,int> > &ReorderedGraphVertId,
585  int &firstNonDirGraphVertId,
586  int &nExtraDirichlet,
587  BottomUpSubStructuredGraphSharedPtr &bottomUpGraph,
588  set<int> &extraDirVerts,
589  const bool checkIfSystemSingular,
590  int mdswitch,
591  bool doInteriorMap)
592  {
593  int i,j,k,l;
594  int cnt = 0;
595  int meshVertId, meshVertId2;
596  int meshEdgeId;
597  int graphVertId = 0;
598  StdRegions::StdExpansion2DSharedPtr locExpansion;
599  LocalRegions::SegExpSharedPtr bndSegExp;
600  MultiRegions::ExpList0DSharedPtr bndVertExp;
601  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
602  map<int,int>::iterator mapIt;
603  map<int,int>::const_iterator mapConstIt;
604  bool systemSingular = true;
605  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
606  PeriodicMap::const_iterator pIt;
607 
608  /**
609  * STEP 1: Order the Dirichlet vertices and edges first
610  */
611  for(i = 0; i < bndCondExp.num_elements(); i++)
612  {
613  // Check to see if any value on edge has Dirichlet value.
614  cnt = 0;
615  for(k = 0; k < bndConditions.num_elements(); ++k)
616  {
617  if (bndConditions[k][i]->GetBoundaryConditionType() ==
618  SpatialDomains::eDirichlet)
619  {
620  cnt++;
621  }
622  if (bndConditions[k][i]->GetBoundaryConditionType() !=
623  SpatialDomains::eNeumann)
624  {
625  systemSingular = false;
626  }
627  }
628 
629  // If all boundaries are Dirichlet take out of mask
630  if(cnt == bndConditions.num_elements())
631  {
632  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
633  {
634  bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
635  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
636  ReorderedGraphVertId[1][meshEdgeId] = graphVertId++;
637  for(k = 0; k < 2; k++)
638  {
639  meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
640  if(ReorderedGraphVertId[0].count(meshVertId) == 0)
641  {
642  ReorderedGraphVertId[0][meshVertId] =
643  graphVertId++;
644  }
645  }
646  }
647  }
648  }
649 
650  /**
651  * STEP 1.4: Check for singular system and add pinning Dirichlet
652  * vertex
653  */
654  // Check between processes if the whole system is singular
655  int n = m_comm->GetSize();
656  int p = m_comm->GetRank();
657  int s = systemSingular ? 1 : 0;
658  vCommRow->AllReduce(s, LibUtilities::ReduceMin);
659  systemSingular = (s == 1 ? true : false);
660 
661  // Count the number of boundary regions on each process
662  Array<OneD, int> bccounts(n, 0);
663  bccounts[p] = bndCondExp.num_elements();
664  vCommRow->AllReduce(bccounts, LibUtilities::ReduceSum);
665 
666  // Find the process rank with the maximum number of boundary regions
667  int maxBCIdx = Vmath::Imax(n, bccounts, 1);
668 
669  // If the system is singular, the process with the maximum number of
670  // BCs will set a Dirichlet vertex to make system non-singular.
671  // Note: we find the process with maximum boundary regions to ensure
672  // we do not try to set a Dirichlet vertex on a partition with no
673  // intersection with the boundary.
674  meshVertId = 0;
675  if(systemSingular == true && checkIfSystemSingular && maxBCIdx == p)
676  {
677  if (m_session->DefinesParameter("SingularElement"))
678  {
679  int s_eid;
680  m_session->LoadParameter("SingularElement", s_eid);
681 
682  ASSERTL1(s_eid < locExpVector.size(),
683  "SingularElement Parameter is too large");
684 
685  meshVertId = locExpVector[s_eid]
686  ->as<LocalRegions::Expansion2D>()
687  ->GetGeom2D()->GetVid(0);
688  }
689  else if (m_session->DefinesParameter("SingularVertex"))
690  {
691  m_session->LoadParameter("SingularVertex", meshVertId);
692  }
693  else if (bndCondExp.num_elements() == 0)
694  {
695  // All boundaries are periodic.
696  meshVertId = locExpVector[0]
697  ->as<LocalRegions::Expansion2D>()
698  ->GetGeom2D()->GetVid(0);
699  }
700  else
701  {
702  bndSegExp = bndCondExp[bndCondExp.num_elements()-1]
703  ->GetExp(0)->as<LocalRegions::SegExp>();
704  meshVertId = bndSegExp->GetGeom1D()->GetVid(0);
705  }
706 
707  if (ReorderedGraphVertId[0].count(meshVertId) == 0)
708  {
709  ReorderedGraphVertId[0][meshVertId] = graphVertId++;
710  }
711  }
712 
713  vCommRow->AllReduce(meshVertId, LibUtilities::ReduceSum);
714 
715  // When running in parallel, we need to ensure that the singular
716  // mesh vertex is communicated to any periodic vertices, otherwise
717  // the system may diverge.
718  if(systemSingular == true && checkIfSystemSingular)
719  {
720  // Firstly, we check that no other processors have this
721  // vertex. If they do, then we mark the vertex as also being
722  // Dirichlet.
723  if (maxBCIdx != p)
724  {
725  if (Dofs[0].count(meshVertId) > 0)
726  {
727  if (ReorderedGraphVertId[0].count(meshVertId) == 0)
728  {
729  ReorderedGraphVertId[0][meshVertId] = graphVertId++;
730  }
731  }
732  }
733 
734  // In the case that meshVertId is periodic with other vertices,
735  // this process and all other processes need to make sure that
736  // the periodic vertices are also marked as Dirichlet.
737  int gId;
738 
739  // At least one process (maxBCidx) will have already associated
740  // a graphVertId with meshVertId. Others won't even have any of
741  // the vertices. The logic below is designed to handle both
742  // cases.
743  if (ReorderedGraphVertId[0].count(meshVertId) == 0)
744  {
745  gId = -1;
746  }
747  else
748  {
749  gId = ReorderedGraphVertId[0][meshVertId];
750  }
751 
752  for (pIt = periodicVerts.begin();
753  pIt != periodicVerts.end(); ++pIt)
754  {
755  // Either the vertex is local to this processor (in which
756  // case it will be in the pIt->first position) or else
757  // meshVertId might be contained within another processor's
758  // vertex list. The if statement below covers both cases. If
759  // we find it, set as Dirichlet with the vertex id gId.
760  if (pIt->first == meshVertId)
761  {
762  ReorderedGraphVertId[0][meshVertId] =
763  gId < 0 ? graphVertId++ : gId;
764 
765  for (i = 0; i < pIt->second.size(); ++i)
766  {
767  if (pIt->second[i].isLocal)
768  {
769  ReorderedGraphVertId[0][pIt->second[i].id] =
770  gId;
771  }
772  }
773  }
774  else
775  {
776  bool found = false;
777  for (i = 0; i < pIt->second.size(); ++i)
778  {
779  if (pIt->second[i].id == meshVertId)
780  {
781  found = true;
782  break;
783  }
784  }
785 
786  if (found)
787  {
788  ReorderedGraphVertId[0][pIt->first] =
789  gId < 0 ? graphVertId++ : gId;
790 
791  for (i = 0; i < pIt->second.size(); ++i)
792  {
793  if (pIt->second[i].isLocal)
794  {
795  ReorderedGraphVertId[0][
796  pIt->second[i].id] = gId;
797  }
798  }
799  }
800  }
801  }
802  }
803 
804  /**
805  * STEP 1.5: Exchange Dirichlet mesh vertices between processes and
806  * check for singular problems.
807  */
808  // Collate information on Dirichlet vertices from all processes
809  Array<OneD, int> counts (n, 0);
810  Array<OneD, int> offsets(n, 0);
811  counts[p] = ReorderedGraphVertId[0].size();
812  vCommRow->AllReduce(counts, LibUtilities::ReduceSum);
813 
814  for (i = 1; i < n; ++i)
815  {
816  offsets[i] = offsets[i-1] + counts[i-1];
817  }
818 
819  int nTot = Vmath::Vsum(n,counts,1);
820  Array<OneD, int> vertexlist(nTot, 0);
821  std::map<int, int>::iterator it;
822  for (it = ReorderedGraphVertId[0].begin(), i = 0;
823  it != ReorderedGraphVertId[0].end();
824  ++it, ++i)
825  {
826  vertexlist[offsets[p] + i] = it->first;
827  }
828  vCommRow->AllReduce(vertexlist, LibUtilities::ReduceSum);
829 
830  map<int, int> extraDirVertIds;
831 
832  // Ensure Dirchlet vertices are consistently recorded between
833  // processes (e.g. Dirichlet region meets Neumann region across a
834  // partition boundary requires vertex on partition to be Dirichlet).
835  for (i = 0; i < n; ++i)
836  {
837  if (i == p)
838  {
839  continue;
840  }
841 
842  for(j = 0; j < locExpVector.size(); j++)
843  {
844  locExpansion = boost::dynamic_pointer_cast<
845  StdRegions::StdExpansion2D>(
846  locExpVector[locExp.GetOffset_Elmt_Id(j)]);
847 
848  for(k = 0; k < locExpansion->GetNverts(); k++)
849  {
850  meshVertId = locExpansion
851  ->as<LocalRegions::Expansion2D>()
852  ->GetGeom2D()->GetVid(k);
853  if (ReorderedGraphVertId[0].count(meshVertId) != 0)
854  {
855  continue;
856  }
857 
858  for (l = 0; l < counts[i]; ++l)
859  {
860  if (vertexlist[offsets[i]+l] == meshVertId)
861  {
862  extraDirVertIds[meshVertId] = i;
863  ReorderedGraphVertId[0][meshVertId] =
864  graphVertId++;
865  nExtraDirichlet++;
866  }
867  }
868  }
869  }
870  }
871 
872  for (i = 0; i < n; ++i)
873  {
874  counts [i] = 0;
875  offsets[i] = 0;
876  }
877 
878  counts[p] = extraDirVertIds.size();
879  vCommRow->AllReduce(counts, LibUtilities::ReduceSum);
880  nTot = Vmath::Vsum(n, counts, 1);
881 
882  offsets[0] = 0;
883 
884  for (i = 1; i < n; ++i)
885  {
886  offsets[i] = offsets[i-1] + counts[i-1];
887  }
888 
889  Array<OneD, int> vertids (nTot, 0);
890  Array<OneD, int> vertprocs(nTot, 0);
891 
892  for (it = extraDirVertIds.begin(), i = 0;
893  it != extraDirVertIds.end(); ++it, ++i)
894  {
895  vertids [offsets[p]+i] = it->first;
896  vertprocs[offsets[p]+i] = it->second;
897  }
898 
899  vCommRow->AllReduce(vertids, LibUtilities::ReduceSum);
900  vCommRow->AllReduce(vertprocs, LibUtilities::ReduceSum);
901 
902  for (i = 0; i < nTot; ++i)
903  {
904  if (m_comm->GetRank() != vertprocs[i])
905  {
906  continue;
907  }
908 
909  extraDirVerts.insert(vertids[i]);
910  }
911 
912  firstNonDirGraphVertId = graphVertId;
913 
914  /**
915  * STEP 2: Now order all other vertices and edges in the graph and
916  * create a temporary numbering of domain-interior vertices and
917  * edges.
918  */
919  typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
920  BoostGraph boostGraphObj;
921 
922  int tempGraphVertId = 0;
923  int nVerts;
924  int vertCnt;
925  int edgeCnt;
926  int localOffset=0;
927  int nTotalVerts=0;
928  map<int, int> vertTempGraphVertId;
929  map<int, int> edgeTempGraphVertId;
930  map<int, int> intTempGraphVertId;
931  Array<OneD, int> localVerts;
932  Array<OneD, int> localEdges;
933  Array<OneD, int> localinterior;
934 
935  m_numLocalBndCoeffs = 0;
936 
937  /// - Local periodic vertices.
938  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
939  {
940  meshVertId = pIt->first;
941 
942  // This periodic vertex is already Dirichlet.
943  if (ReorderedGraphVertId[0].count(pIt->first) != 0)
944  {
945  for (i = 0; i < pIt->second.size(); ++i)
946  {
947  meshVertId2 = pIt->second[i].id;
948  if (ReorderedGraphVertId[0].count(meshVertId2) == 0 &&
949  pIt->second[i].isLocal)
950  {
951  ReorderedGraphVertId[0][meshVertId2] =
952  ReorderedGraphVertId[0][meshVertId];
953  }
954  }
955  continue;
956  }
957 
958  // One of the attached vertices is Dirichlet.
959  bool isDirichlet = false;
960  for (i = 0; i < pIt->second.size(); ++i)
961  {
962  if (!pIt->second[i].isLocal)
963  {
964  continue;
965  }
966 
967  meshVertId2 = pIt->second[i].id;
968  if (ReorderedGraphVertId[0].count(meshVertId2) > 0)
969  {
970  isDirichlet = true;
971  break;
972  }
973  }
974 
975  if (isDirichlet)
976  {
977  ReorderedGraphVertId[0][meshVertId] =
978  ReorderedGraphVertId[0][pIt->second[i].id];
979 
980  for (j = 0; j < pIt->second.size(); ++j)
981  {
982  meshVertId2 = pIt->second[i].id;
983  if (j == i || !pIt->second[j].isLocal ||
984  ReorderedGraphVertId[0].count(meshVertId2) > 0)
985  {
986  continue;
987  }
988 
989  ReorderedGraphVertId[0][meshVertId2] =
990  ReorderedGraphVertId[0][pIt->second[i].id];
991  }
992 
993  continue;
994  }
995 
996  // Otherwise, see if a vertex ID has already been set.
997  for (i = 0; i < pIt->second.size(); ++i)
998  {
999  if (!pIt->second[i].isLocal)
1000  {
1001  continue;
1002  }
1003 
1004  if (vertTempGraphVertId.count(pIt->second[i].id) > 0)
1005  {
1006  break;
1007  }
1008  }
1009 
1010  if (i == pIt->second.size())
1011  {
1012  vertTempGraphVertId[meshVertId] = tempGraphVertId++;
1013  }
1014  else
1015  {
1016  vertTempGraphVertId[meshVertId] = vertTempGraphVertId[pIt->second[i].id];
1017  }
1018  }
1019 
1020  /// - Local periodic edges. !! SJS: This will need pushing
1021  /// - further down for block precodnitioner
1022  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1023  {
1024  if (!pIt->second[0].isLocal)
1025  {
1026  // The edge mapped to is on another process.
1027  meshEdgeId = pIt->first;
1028  ASSERTL0(ReorderedGraphVertId[1].count(meshEdgeId) == 0,
1029  "This periodic boundary edge has been specified before");
1030  edgeTempGraphVertId[meshEdgeId] = tempGraphVertId++;
1031  }
1032  else if (pIt->first < pIt->second[0].id)
1033  {
1034  ASSERTL0(ReorderedGraphVertId[1].count(pIt->first) == 0,
1035  "This periodic boundary edge has been specified before");
1036  ASSERTL0(ReorderedGraphVertId[1].count(pIt->second[0].id) == 0,
1037  "This periodic boundary edge has been specified before");
1038 
1039  edgeTempGraphVertId[pIt->first] = tempGraphVertId;
1040  edgeTempGraphVertId[pIt->second[0].id] = tempGraphVertId++;
1041  }
1042  }
1043 
1044  /// - All other vertices and edges
1045  int elmtid;
1046  for(i = 0; i < locExpVector.size(); ++i)
1047  {
1048  elmtid = locExp.GetOffset_Elmt_Id(i);
1049  if((locExpansion = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D>(
1050  locExpVector[elmtid])))
1051  {
1052  m_numLocalBndCoeffs += locExpansion->NumBndryCoeffs();
1053 
1054  nTotalVerts += locExpansion->GetNverts();
1055  }
1056  }
1057 
1058  // Store the temporary graph vertex
1059  // id's of all element edges and
1060  // vertices in these 2 arrays below
1061  localVerts = Array<OneD, int>(nTotalVerts,-1);
1062  localEdges = Array<OneD, int>(nTotalVerts,-1);
1063 
1064  for(i = 0; i < locExpVector.size(); ++i)
1065  {
1066  elmtid = locExp.GetOffset_Elmt_Id(i);
1067  if((locExpansion = locExpVector[elmtid]
1068  ->as<StdRegions::StdExpansion2D>()))
1069  {
1070  vertCnt = 0;
1071  nVerts = locExpansion->GetNverts();
1072  for(j = 0; j < nVerts; ++j)
1073  {
1074  meshVertId = locExpansion
1075  ->as<LocalRegions::Expansion2D>()
1076  ->GetGeom2D()->GetVid(j);
1077  if(ReorderedGraphVertId[0].count(meshVertId) == 0)
1078  {
1079  // non-periodic & non-Dirichlet vertex
1080  if(vertTempGraphVertId.count(meshVertId) == 0)
1081  {
1082  boost::add_vertex(boostGraphObj);
1083  vertTempGraphVertId[meshVertId] = tempGraphVertId++;
1084  }
1085  localVerts[localOffset + vertCnt++] = vertTempGraphVertId[meshVertId];
1086  }
1087  }
1088  }
1089  localOffset+=nVerts;
1090  }
1091 
1092  m_numNonDirVertexModes=tempGraphVertId;
1093 
1094  localOffset=0;
1095  for(i = 0; i < locExpVector.size(); ++i)
1096  {
1097  elmtid = locExp.GetOffset_Elmt_Id(i);
1098  if((locExpansion = locExpVector[elmtid]
1099  ->as<StdRegions::StdExpansion2D>()))
1100  {
1101  edgeCnt = 0;
1102  nVerts = locExpansion->GetNverts();
1103  for(j = 0; j < nVerts; ++j)
1104  {
1105  meshEdgeId = locExpansion
1106  ->as<LocalRegions::Expansion2D>()
1107  ->GetGeom2D()->GetEid(j);
1108  if(ReorderedGraphVertId[1].count(meshEdgeId) == 0)
1109  {
1110  // non-periodic & non-Dirichlet edge
1111  if(edgeTempGraphVertId.count(meshEdgeId) == 0)
1112  {
1113  boost::add_vertex(boostGraphObj);
1114  edgeTempGraphVertId[meshEdgeId] = tempGraphVertId++;
1115  }
1116  localEdges[localOffset + edgeCnt++] = edgeTempGraphVertId[meshEdgeId];
1117  }
1118  }
1119  }
1120  localOffset+=nVerts;
1121  }
1122 
1123  // Number of dirichlet edges
1124  m_numNonDirEdges = edgeTempGraphVertId.size();
1125 
1126  if(doInteriorMap)
1127  {
1128  for(i = 0; i < locExpVector.size(); ++i)
1129  {
1130  elmtid = locExp.GetOffset_Elmt_Id(i);
1131  if((locExpansion = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D>(
1132  locExpVector[elmtid])))
1133  {
1134 
1135  boost::add_vertex(boostGraphObj);
1136  intTempGraphVertId[elmtid] = tempGraphVertId++;
1137  }
1138  }
1139  }
1140 
1141  localOffset=0;
1142  for(i = 0; i < locExpVector.size(); ++i)
1143  {
1144  elmtid = locExp.GetOffset_Elmt_Id(i);
1145  if((locExpansion = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D>(
1146  locExpVector[elmtid])))
1147  {
1148  nVerts = locExpansion->GetNverts();
1149  // Now loop over all local edges and vertices of this
1150  // element and define that all other edges and vertices of
1151  // this element are adjacent to them.
1152  for(j = 0; j < nVerts; j++)
1153  {
1154  if(localVerts[j+localOffset]==-1)
1155  {
1156  break;
1157  }
1158 
1159  for(k = 0; k < nVerts; k++)
1160  {
1161  if(localVerts[k+localOffset]==-1)
1162  {
1163  break;
1164  }
1165  if(k!=j)
1166  {
1167  boost::add_edge( (size_t) localVerts[j+localOffset],
1168  (size_t) localVerts[k+localOffset],boostGraphObj);
1169  }
1170  }
1171 
1172  for(k = 0; k < nVerts; k++)
1173  {
1174  if(localEdges[k+localOffset]==-1)
1175  {
1176  break;
1177  }
1178  boost::add_edge( (size_t) localVerts[j+localOffset],
1179  (size_t) localEdges[k+localOffset],boostGraphObj);
1180  }
1181 
1182  if(doInteriorMap)
1183  {
1184  boost::add_edge( (size_t) localVerts[j+localOffset],
1185  (size_t) intTempGraphVertId[elmtid],boostGraphObj);
1186  }
1187  }
1188 
1189  for(j = 0; j < nVerts; j++)
1190  {
1191  if(localEdges[j+localOffset]==-1)
1192  {
1193  break;
1194  }
1195  for(k = 0; k < nVerts; k++)
1196  {
1197  if(localEdges[k+localOffset]==-1)
1198  {
1199  break;
1200  }
1201  if(k!=j)
1202  {
1203  boost::add_edge( (size_t) localEdges[j+localOffset],
1204  (size_t) localEdges[k+localOffset],boostGraphObj);
1205  }
1206  }
1207  for(k = 0; k < nVerts; k++)
1208  {
1209  if(localVerts[k+localOffset]==-1)
1210  {
1211  break;
1212  }
1213  boost::add_edge( (size_t) localEdges[j+localOffset],
1214  (size_t) localVerts[k+localOffset],boostGraphObj);
1215  }
1216 
1217  if(doInteriorMap)
1218  {
1219  boost::add_edge( (size_t) localEdges[j+localOffset],
1220  (size_t) intTempGraphVertId[elmtid], boostGraphObj);
1221  }
1222  }
1223 
1224  if(doInteriorMap)
1225  {
1226  for(j = 0; j < nVerts; j++)
1227  {
1228  if(localVerts[j+localOffset]==-1)
1229  {
1230  break;
1231  }
1232  boost::add_edge( (size_t) intTempGraphVertId[elmtid],
1233  (size_t) localVerts[j+localOffset], boostGraphObj);
1234  }
1235 
1236  for(j = 0; j < nVerts; j++)
1237  {
1238  if(localEdges[j+localOffset]==-1)
1239  {
1240  break;
1241  }
1242 
1243  boost::add_edge( (size_t) intTempGraphVertId[elmtid],
1244  (size_t) localEdges[j+localOffset], boostGraphObj);
1245  }
1246  }
1247 
1248 
1249  }
1250  else
1251  {
1252  ASSERTL0(false,"dynamic cast to a local 2D expansion failed");
1253  }
1254  localOffset+=nVerts;
1255  }
1256 
1257  // Container to store vertices of the graph which correspond to
1258  // degrees of freedom along the boundary.
1259  set<int> partVerts;
1260 
1261  if (m_solnType == eIterativeMultiLevelStaticCond)
1262  {
1263  vector<long> procVerts, procEdges;
1264  set <int> foundVerts, foundEdges;
1265 
1266  // Loop over element and construct the procVerts and procEdges
1267  // vectors, which store the geometry IDs of mesh vertices and
1268  // edges respectively which are local to this process.
1269  for(i = cnt = 0; i < locExpVector.size(); ++i)
1270  {
1271  elmtid = locExp.GetOffset_Elmt_Id(i);
1272  if((locExpansion = boost::dynamic_pointer_cast<
1273  StdRegions::StdExpansion2D>(locExpVector[elmtid])))
1274  {
1275  for (j = 0; j < locExpansion->GetNverts(); ++j, ++cnt)
1276  {
1277  int vid = locExpansion
1278  ->as<LocalRegions::Expansion2D>()
1279  ->GetGeom2D()->GetVid(j)+1;
1280  int eid = locExpansion
1281  ->as<LocalRegions::Expansion2D>()
1282  ->GetGeom2D()->GetEid(j)+1;
1283 
1284  if (foundVerts.count(vid) == 0)
1285  {
1286  procVerts.push_back(vid);
1287  foundVerts.insert(vid);
1288  }
1289 
1290  if (foundEdges.count(eid) == 0)
1291  {
1292  procEdges.push_back(eid);
1293  foundEdges.insert(eid);
1294  }
1295  }
1296  }
1297  else
1298  {
1299  ASSERTL0(false,
1300  "dynamic cast to a local 2D expansion failed");
1301  }
1302  }
1303 
1304  int unique_verts = foundVerts.size();
1305  int unique_edges = foundEdges.size();
1306 
1307  // Now construct temporary GS objects. These will be used to
1308  // populate the arrays tmp3 and tmp4 with the multiplicity of
1309  // the vertices and edges respectively to identify those
1310  // vertices and edges which are located on partition boundary.
1311  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1312  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1313  Gs::gs_data *tmp1 = Gs::Init(vertArray, m_comm);
1314  Gs::gs_data *tmp2 = Gs::Init(edgeArray, m_comm);
1315  Array<OneD, NekDouble> tmp3(unique_verts, 1.0);
1316  Array<OneD, NekDouble> tmp4(unique_edges, 1.0);
1317  Gs::Gather(tmp3, Gs::gs_add, tmp1);
1318  Gs::Gather(tmp4, Gs::gs_add, tmp2);
1319 
1320  // Finally, fill the partVerts set with all non-Dirichlet
1321  // vertices which lie on a partition boundary.
1322  for (i = 0; i < unique_verts; ++i)
1323  {
1324  if (tmp3[i] > 1.0)
1325  {
1326  if (ReorderedGraphVertId[0].count(procVerts[i]-1) == 0)
1327  {
1328  partVerts.insert(vertTempGraphVertId[procVerts[i]-1]);
1329  }
1330  }
1331  }
1332 
1333  for (i = 0; i < unique_edges; ++i)
1334  {
1335  if (tmp4[i] > 1.0)
1336  {
1337  if (ReorderedGraphVertId[1].count(procEdges[i]-1) == 0)
1338  {
1339  partVerts.insert(edgeTempGraphVertId[procEdges[i]-1]);
1340  }
1341  }
1342  }
1343  }
1344 
1345  /**
1346  * STEP 3: Reorder graph for optimisation.
1347  */
1348  int nGraphVerts = tempGraphVertId;
1349  Array<OneD, int> perm(nGraphVerts);
1350  Array<OneD, int> iperm(nGraphVerts);
1351 
1352  if(nGraphVerts)
1353  {
1354  switch(m_solnType)
1355  {
1356  case eDirectFullMatrix:
1357  case eIterativeFull:
1358  case eIterativeStaticCond:
1359  case eXxtFullMatrix:
1360  case eXxtStaticCond:
1361  case ePETScFullMatrix:
1362  case ePETScStaticCond:
1363  {
1364  NoReordering(boostGraphObj,perm,iperm);
1365  }
1366  break;
1367  case eDirectStaticCond:
1368  {
1369  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1370  }
1371  break;
1372  case eDirectMultiLevelStaticCond:
1373  case eIterativeMultiLevelStaticCond:
1374  case eXxtMultiLevelStaticCond:
1375  case ePETScMultiLevelStaticCond:
1376  {
1377  MultiLevelBisectionReordering(boostGraphObj,perm,iperm,bottomUpGraph,partVerts,mdswitch);
1378  }
1379  break;
1380  default:
1381  {
1382  ASSERTL0(false,"Unrecognised solution type: " + std::string(MultiRegions::GlobalSysSolnTypeMap[m_solnType]));
1383  }
1384  }
1385  }
1386 
1387  // For parallel multi-level static condensation determine the lowest
1388  // static condensation level amongst processors.
1389  if (m_solnType == eIterativeMultiLevelStaticCond)
1390  {
1391  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1392  vCommRow->AllReduce(m_lowestStaticCondLevel,
1393  LibUtilities::ReduceMax);
1394  }
1395  else
1396  {
1397  m_lowestStaticCondLevel = 0;
1398  }
1399 
1400  /**
1401  * STEP 4: Fill the #vertReorderedGraphVertId and
1402  * #edgeReorderedGraphVertId with the optimal ordering from boost.
1403  */
1404  for(mapIt = vertTempGraphVertId.begin(); mapIt != vertTempGraphVertId.end(); mapIt++)
1405  {
1406  ReorderedGraphVertId[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
1407  }
1408  for(mapIt = edgeTempGraphVertId.begin(); mapIt != edgeTempGraphVertId.end(); mapIt++)
1409  {
1410  ReorderedGraphVertId[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
1411  }
1412 
1413  if(doInteriorMap)
1414  {
1415  for(mapIt = intTempGraphVertId.begin(); mapIt != intTempGraphVertId.end(); mapIt++)
1416  {
1417  ReorderedGraphVertId[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
1418  }
1419  }
1420 
1421  return nGraphVerts;
1422  }
1423  }
1424 }
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2