Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CoupledLocalToGlobalC0ContMap.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CoupledLcoalToGlobalC0ContMap.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: Wrapper class around the library
33 // LocalToGlobalC0ContMap class for use in the Couplied Linearised NS
34 // solver.
35 ///////////////////////////////////////////////////////////////////////////////
36 
39 #include <LocalRegions/SegExp.h>
43 
44 using namespace std;
45 
46 namespace Nektar
47 {
48  /**
49  * This is an vector extension of
50  * MultiRegions::AssemblyMapCG::SetUp2DExpansionC0ContMap related to the
51  * Linearised Navier Stokes problem
52  */
53  CoupledLocalToGlobalC0ContMap::CoupledLocalToGlobalC0ContMap(
56  const SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions,
58  const MultiRegions::ExpListSharedPtr &pressure,
59  const int nz_loc,
60  const bool CheckforSingularSys):
61  AssemblyMapCG(pSession)
62  {
63  int i,j,k,n;
64  int cnt = 0,offset=0;
65  int meshVertId;
66  int meshEdgeId;
67  int bndEdgeCnt;
68  int globalId;
69  int nEdgeCoeffs;
70  int nEdgeInteriorCoeffs;
71  int firstNonDirGraphVertId;
72  int nLocBndCondDofs = 0;
73  int nLocDirBndCondDofs = 0;
74  int nExtraDirichlet = 0;
78  StdRegions::Orientation edgeOrient;
79  Array<OneD, unsigned int> edgeInteriorMap;
80  Array<OneD, int> edgeInteriorSign;
81  int nvel = fields.num_elements();
82  MultiRegions::PeriodicMap::const_iterator pIt;
83 
84  const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
85  int id, diff;
86  int nel = fields[0]->GetNumElmts();
87 
88  MultiRegions::PeriodicMap periodicVerts;
89  MultiRegions::PeriodicMap periodicEdges;
90  MultiRegions::PeriodicMap periodicFaces;
91  vector<map<int,int> > ReorderedGraphVertId(3);
93  int staticCondLevel = 0;
94 
95  if(CheckforSingularSys) //all singularity checking by setting flag to true
96  {
97  m_systemSingular = true;
98  }
99  else // Turn off singular checking by setting flag to false
100  {
101  m_systemSingular = false;
102  }
103 
104  /**
105  * STEP 1: Wrap boundary conditions vector in an array
106  * (since routine is set up for multiple fields) and call
107  * the graph re-odering subroutine to obtain the reordered
108  * values
109  */
110 
111  // Obtain any periodic information and allocate default mapping array
112  fields[0]->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
113 
114 
115  const Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp = fields[0]->GetBndCondExpansions();
116 
118  for(i = 0; i < nvel; ++i)
119  {
120  bndConditionsVec[i] = fields[i]->GetBndConditions();
121  }
122 
123  map<int,int> IsDirVertDof;
124  map<int,int> IsDirEdgeDof;
126 
128  for(j = 0; j < bndCondExp.num_elements(); ++j)
129  {
130  map<int,int> BndExpVids;
131  // collect unique list of vertex ids for this expansion
132  for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
133  {
134  g = bndCondExp[j]->GetExp(k)->as<LocalRegions::Expansion1D>()
135  ->GetGeom1D();
136  BndExpVids[g->GetVid(0)] = g->GetVid(0);
137  BndExpVids[g->GetVid(1)] = g->GetVid(1);
138  }
139 
140  for(i = 0; i < nvel; ++i)
141  {
142  if(bndConditionsVec[i][j]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
143  {
144  // set number of Dirichlet conditions along edge
145  for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
146  {
147  IsDirEdgeDof[bndCondExp[j]->GetExp(k)
149  ->GetGeom1D()->GetEid()] += 1;
150  }
151 
152 
153  // Set number of Dirichlet conditions at vertices
154  // with a clamp on its maximum value being nvel to
155  // handle corners between expansions
156  for(mapIt = BndExpVids.begin(); mapIt != BndExpVids.end(); mapIt++)
157  {
158  id = IsDirVertDof[mapIt->second]+1;
159  IsDirVertDof[mapIt->second] = (id > nvel)?nvel:id;
160  }
161  }
162  else
163  {
164  // Check to see that edge normals have non-zero
165  // component in this direction since otherwise
166  // also can be singular.
167  /// @TODO: Fix this so that we can extract normals from edges
168  for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
169  {
172  = bndCondExp[j]->GetExp(k)
174  locnorm = loc_exp->GetLeftAdjacentElementExp()->GetEdgeNormal(loc_exp->GetLeftAdjacentElementEdge());
175  //locnorm = bndCondExp[j]->GetExp(k)->Get GetMetricInfo()->GetNormal();
176 
177  int ndir = locnorm.num_elements();
178  if(i < ndir) // account for Fourier version where n can be larger then ndir
179  {
180  for(int l = 0; l < locnorm[0].num_elements(); ++l)
181  {
182  if(fabs(locnorm[i][l]) > NekConstants::kNekZeroTol)
183  {
184  m_systemSingular = false;
185  break;
186  }
187  }
188  }
189  if(m_systemSingular == false)
190  {
191  break;
192  }
193  }
194  }
195  }
196  }
197 
198  Array<OneD, map<int,int> >Dofs(2);
199 
200  Array<OneD, int> AddMeanPressureToEdgeId(nel,-1);
201  int edgeId,vertId;
202 
203 
204  // special case of singular problem - need to fix one pressure
205  // dof to a dirichlet edge. Since we attached pressure dof to
206  // last velocity component of edge need to make sure this
207  // component is Dirichlet
208  if(m_systemSingular)
209  {
210  id = -1;
211  for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
212  {
213  if(bndConditionsVec[nvel-1][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
214  {
215  id = bndCondExp[i]->GetExp(0)
216  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
217  ->GetEid();
218  break;
219  }
220  }
221 
222  ASSERTL0(id != -1," Did not find an edge to attach singular pressure degree of freedom");
223 
224  // determine element with this edge id. There may be a
225  // more direct way of getting element from spatialDomains
226  for(i = 0; i < nel; ++i)
227  {
228  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
229  {
230  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
231  ->GetGeom2D())->GetEid(j);
232 
233  if(edgeId == id)
234  {
235  AddMeanPressureToEdgeId[i] = id;
236  break;
237  }
238  }
239 
240  if(AddMeanPressureToEdgeId[i] != -1)
241  {
242  break;
243  }
244  }
245  }
246 
247 
248  for(i = 0; i < nel; ++i)
249  {
250  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
251  {
252  vertId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
253  ->GetGeom2D())->GetVid(j);
254  if(Dofs[0].count(vertId) == 0)
255  {
256  Dofs[0][vertId] = nvel*nz_loc;
257 
258  // Adjust for a Dirichlet boundary condition to give number to be solved
259  if(IsDirVertDof.count(vertId) != 0)
260  {
261  Dofs[0][vertId] -= IsDirVertDof[vertId]*nz_loc;
262  }
263  }
264 
265  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
266  ->GetGeom2D())->GetEid(j);
267  if(Dofs[1].count(edgeId) == 0)
268  {
269  Dofs[1][edgeId] = nvel*(locExpVector[i]->GetEdgeNcoeffs(j)-2)*nz_loc;
270  }
271 
272  // Adjust for Dirichlet boundary conditions to give number to be solved
273  if(IsDirEdgeDof.count(edgeId) != 0)
274  {
275  Dofs[1][edgeId] -= IsDirEdgeDof[edgeId]*nz_loc*(locExpVector[i]->GetEdgeNcoeffs(j)-2);
276  }
277  }
278  }
279 
280  set<int> extraDirVerts, extraDirEdges;
281 
282  CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false,
283  periodicVerts, periodicEdges, periodicFaces,
284  ReorderedGraphVertId, bottomUpGraph, extraDirVerts,
285  extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, 4);
286  /*
287  SetUp2DGraphC0ContMap(*fields[0],
288  bndCondExp,
289  bndConditionsVec,
290  periodicVerts, periodicEdges,
291  Dofs, ReorderedGraphVertId,
292  firstNonDirGraphVertId, nExtraDirichlet,
293  bottomUpGraph, extraDir, false, 4);
294  */
295 
296  /**
297  * STEP 2a: Set the mean pressure modes to edges depending on
298  * type of direct solver technique;
299  */
300 
301  // determine which edge to add mean pressure dof based on
302  // ensuring that at least one pressure dof from an internal
303  // patch is associated with its boundary system
304  if(m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
305  {
306 
307 
308  FindEdgeIdToAddMeanPressure(ReorderedGraphVertId,
309  nel, locExpVector,
310  edgeId, vertId, firstNonDirGraphVertId, IsDirEdgeDof,
311  bottomUpGraph,
312  AddMeanPressureToEdgeId);
313  }
314 
315  // Set unset elmts to non-Dirichlet edges.
316  // special case of singular problem - need to fix one
317  // pressure dof to a dirichlet edge
318  for(i = 0; i < nel; ++i)
319  {
320  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
321  {
322  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
323  ->GetGeom2D())->GetEid(j);
324 
325  if(IsDirEdgeDof.count(edgeId) == 0) // interior edge
326  {
327  // setup AddMeanPressureToEdgeId to decide where to
328  // put pressure
329  if(AddMeanPressureToEdgeId[i] == -1)
330  {
331  AddMeanPressureToEdgeId[i] = edgeId;
332  }
333  }
334  }
335  ASSERTL0((AddMeanPressureToEdgeId[i] != -1),"Did not determine "
336  "an edge to attach mean pressure dof");
337  // Add the mean pressure degree of freedom to this edge
338  Dofs[1][AddMeanPressureToEdgeId[i]] += nz_loc;
339  }
340 
341  map<int,int> pressureEdgeOffset;
342 
343  /**
344  * STEP 2: Count out the number of Dirichlet vertices and edges first
345  */
346  for(i = 0; i < bndCondExp.num_elements(); i++)
347  {
348  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
349  {
350  bndSegExp = bndCondExp[i]->GetExp(j)
351  ->as<LocalRegions::SegExp>();
352  for(k = 0; k < nvel; ++k)
353  {
354  if(bndConditionsVec[k][i]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
355  {
356  nLocDirBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
357  }
358  nLocBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
359  }
360  }
361  }
362 
363  if(m_systemSingular)
364  {
365  m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet+nz_loc;
366  }
367  else
368  {
369  m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet;
370  }
371 
372  /**
373  * STEP 3: Set up an array which contains the offset information of
374  * the different graph vertices.
375  *
376  * This basically means to identify how many global degrees of
377  * freedom the individual graph vertices correspond. Obviously,
378  * the graph vertices corresponding to the mesh-vertices account
379  * for a single global DOF. However, the graph vertices
380  * corresponding to the element edges correspond to 2*(N-2) global DOF
381  * where N is equal to the number of boundary modes on this edge.
382  */
383  Array<OneD, int> graphVertOffset(nvel*nz_loc*(ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),0);
384  graphVertOffset[0] = 0;
385 
386  m_signChange = false;
387 
388  for(i = 0; i < nel; ++i)
389  {
390  locExpansion = locExpVector[i]->as<StdRegions::StdExpansion2D>();
391 
392  for(j = 0; j < locExpansion->GetNedges(); ++j)
393  {
394  nEdgeCoeffs = locExpansion->GetEdgeNcoeffs(j);
395  meshEdgeId = (locExpansion->as<LocalRegions::Expansion2D>()
396  ->GetGeom2D())->GetEid(j);
397  meshVertId = (locExpansion->as<LocalRegions::Expansion2D>()
398  ->GetGeom2D())->GetVid(j);
399 
400  for(k = 0; k < nvel*nz_loc; ++k)
401  {
402  graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+k] = 1;
403  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+k] = (nEdgeCoeffs-2);
404  }
405 
406  bType = locExpansion->GetEdgeBasisType(j);
407  // need a sign vector for modal expansions if nEdgeCoeffs >=4
408  if( (nEdgeCoeffs >= 4)&&
409  ( (bType == LibUtilities::eModified_A)||
410  (bType == LibUtilities::eModified_B) ) )
411  {
412  m_signChange = true;
413  }
414  }
415  }
416 
417  // Add mean pressure modes;
418  for(i = 0; i < nel; ++i)
419  {
420  graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1] += nz_loc;
421  //graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc] += nz_loc;
422  }
423 
424  // Negate the vertices and edges with only a partial
425  // Dirichlet conditon. Essentially we check to see if an edge
426  // has a mixed Dirichlet with Neumann/Robin Condition and if
427  // so negate the offset associated with this vertex.
428 
429  map<int,int> DirVertChk;
430 
431  for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
432  {
433  cnt = 0;
434  for(j = 0; j < nvel; ++j)
435  {
436  if(bndConditionsVec[j][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
437  {
438  cnt ++;
439  }
440  }
441 
442  // Case where partial Dirichlet boundary condition
443  if((cnt > 0)&&(cnt < nvel))
444  {
445  for(j = 0; j < nvel; ++j)
446  {
447  if(bndConditionsVec[j][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
448  {
449  //negate graph offsets which should be
450  //Dirichlet conditions
451  for(k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
452  {
453  // vertices with mix condition;
454  id = bndCondExp[i]->GetExp(k)
456  ->GetGeom1D()->GetVid(0);
457  if(DirVertChk.count(id*nvel+j) == 0)
458  {
459  DirVertChk[id*nvel+j] = 1;
460  for(n = 0; n < nz_loc; ++n)
461  {
462  graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc + n] *= -1;
463  }
464  }
465 
466  id = bndCondExp[i]->GetExp(k)
468  ->GetGeom1D()->GetVid(1);
469  if(DirVertChk.count(id*nvel+j) == 0)
470  {
471  DirVertChk[id*nvel+j] = 1;
472  for(n = 0; n < nz_loc; ++n)
473  {
474  graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc+n] *= -1;
475  }
476  }
477 
478  // edges with mixed id;
479  id = bndCondExp[i]->GetExp(k)
481  ->GetGeom1D()->GetEid();
482  for(n = 0; n < nz_loc; ++n)
483  {
484  graphVertOffset[ReorderedGraphVertId[1][id]*nvel*nz_loc+j*nz_loc +n] *= -1;
485  }
486  }
487  }
488  }
489  }
490  }
491 
492 
493  cnt = 0;
494  // assemble accumulative list of full Dirichlet values.
495  for(i = 0; i < firstNonDirGraphVertId*nvel*nz_loc; ++i)
496  {
497  diff = abs(graphVertOffset[i]);
498  graphVertOffset[i] = cnt;
499  cnt += diff;
500  }
501 
502  // set Dirichlet values with negative values to Dirichlet value
503  for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
504  {
505  if(graphVertOffset[i] < 0)
506  {
507  diff = -graphVertOffset[i];
508  graphVertOffset[i] = -cnt;
509  cnt += diff;
510  }
511  }
512 
513  // Accumulate all interior degrees of freedom with positive values
515 
516  // offset values
517  for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
518  {
519  if(graphVertOffset[i] >= 0)
520  {
521  diff = graphVertOffset[i];
522  graphVertOffset[i] = cnt;
523  cnt += diff;
524  }
525  }
526 
527  // Finally set negative entries (corresponding to Dirichlet
528  // values ) to be positive
529  for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
530  {
531  if(graphVertOffset[i] < 0)
532  {
533  graphVertOffset[i] = -graphVertOffset[i];
534  }
535  }
536 
537 
538  // Allocate the proper amount of space for the class-data and fill
539  // information that is already known
540  cnt = 0;
542  m_numLocalCoeffs = 0;
543 
544  for(i = 0; i < nel; ++i)
545  {
546  m_numLocalBndCoeffs += nz_loc*(nvel*locExpVector[i]->NumBndryCoeffs() + 1);
547  // add these coeffs up separately since
548  // pressure->GetNcoeffs can include the coefficient in
549  // multiple planes.
550  m_numLocalCoeffs += (pressure->GetExp(i)->GetNcoeffs()-1)*nz_loc;
551  }
552 
554 
555 
559 
560 
561  // Set default sign array.
564 
565  m_staticCondLevel = staticCondLevel;
566  m_numPatches = nel;
567 
570 
571  for(i = 0; i < nel; ++i)
572  {
573  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) nz_loc*(nvel*locExpVector[i]->NumBndryCoeffs() + 1);
574  m_numLocalIntCoeffsPerPatch[i] = (unsigned int) nz_loc*(pressure->GetExp(i)->GetNcoeffs()-1);
575  }
576 
577  /**
578  * STEP 4: Now, all ingredients are ready to set up the actual
579  * local to global mapping.
580  *
581  * The remainder of the map consists of the element-interior
582  * degrees of freedom. This leads to the block-diagonal submatrix
583  * as each element-interior mode is globally orthogonal to modes
584  * in all other elements.
585  */
586  cnt = 0;
587  int nv,velnbndry;
589 
590 
591  // Loop over all the elements in the domain in shuffled
592  // ordering (element type consistency)
593  for(i = 0; i < nel; ++i)
594  {
595  locExpansion = locExpVector[i]->as<StdRegions::StdExpansion2D>();
596 
597  velnbndry = locExpansion->NumBndryCoeffs();
598 
599  // require an inverse ordering of the bmap system to store
600  // local numbering system which takes matrix these
601  // matrices. Therefore get hold of elemental bmap and set
602  // up an inverse map
603  map<int,int> inv_bmap;
604  locExpansion->GetBoundaryMap(bmap);
605  for(j = 0; j < bmap.num_elements(); ++j)
606  {
607  inv_bmap[bmap[j]] = j;
608  }
609 
610  // Loop over all edges (and vertices) of element i
611  for(j = 0; j < locExpansion->GetNedges(); ++j)
612  {
613  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
614  edgeOrient = (locExpansion->as<LocalRegions::Expansion2D>()
615  ->GetGeom2D())->GetEorient(j);
616  meshEdgeId = (locExpansion->as<LocalRegions::Expansion2D>()
617  ->GetGeom2D())->GetEid(j);
618  meshVertId = (locExpansion->as<LocalRegions::Expansion2D>()
619  ->GetGeom2D())->GetVid(j);
620 
621  pIt = periodicEdges.find(meshEdgeId);
622 
623  // See if this edge is periodic. If it is, then we map all
624  // edges to the one with lowest ID, and align all
625  // coefficients to this edge orientation.
626  if (pIt != periodicEdges.end())
627  {
628  pair<int, StdRegions::Orientation> idOrient =
630  meshEdgeId, edgeOrient, pIt->second);
631  edgeOrient = idOrient.second;
632  }
633 
634  locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
635  // Set the global DOF for vertex j of element i
636 
637  for(nv = 0; nv < nvel*nz_loc; ++nv)
638  {
639  m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[locExpansion->GetVertexMap(j)]] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+ nv];
640  // Set the global DOF's for the interior modes of edge j
641  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
642  {
643  m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[edgeInteriorMap[k]]] = graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv]+k;
644  }
645  }
646 
647  // Fill the sign vector if required
648  if(m_signChange)
649  {
650  for(nv = 0; nv < nvel*nz_loc; ++nv)
651  {
652  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
653  {
654  m_localToGlobalSign[cnt+nv*velnbndry + inv_bmap[edgeInteriorMap[k]]] = (NekDouble) edgeInteriorSign[k];
655  }
656  }
657  }
658  }
659 
660  // use difference between two edges of the AddMeanPressureEdgeId to det nEdgeInteriorCoeffs.
661  nEdgeInteriorCoeffs = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc+1] - graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc];
662 
663  int psize = pressure->GetExp(i)->GetNcoeffs();
664  for(n = 0; n < nz_loc; ++n)
665  {
666  m_localToGlobalMap[cnt + nz_loc*nvel*velnbndry + n*psize] = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1]+nEdgeInteriorCoeffs + pressureEdgeOffset[AddMeanPressureToEdgeId[i]];
667 
668  pressureEdgeOffset[AddMeanPressureToEdgeId[i]] += 1;
669  }
670 
671  cnt += (velnbndry*nvel+ psize)*nz_loc;
672  }
673 
674  // Set up the mapping for the boundary conditions
675  offset = cnt = 0;
676  for(nv = 0; nv < nvel; ++nv)
677  {
678  for(i = 0; i < bndCondExp.num_elements(); i++)
679  {
680  for(n = 0; n < nz_loc; ++n)
681  {
682  int ncoeffcnt = 0;
683  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
684  {
685  bndSegExp = bndCondExp[i]->GetExp(j)
686  ->as<LocalRegions::SegExp>();
687 
688  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
689  for(k = 0; k < 2; k++)
690  {
691  meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
692  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+nv*nz_loc+n];
693  }
694 
695  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
696  bndEdgeCnt = 0;
697  nEdgeCoeffs = bndSegExp->GetNcoeffs();
698  for(k = 0; k < nEdgeCoeffs; k++)
699  {
700  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
701  {
702  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
703  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv*nz_loc+n]+bndEdgeCnt;
704  bndEdgeCnt++;
705  }
706  }
707  ncoeffcnt += nEdgeCoeffs;
708  }
709  // Note: Can not use bndCondExp[i]->GetNcoeffs()
710  // due to homogeneous extension not returning just
711  // the value per plane
712  offset += ncoeffcnt;
713  }
714  }
715  }
716 
717  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
718  m_numGlobalBndCoeffs = globalId;
719 
720  /**
721  * STEP 5: The boundary condition mapping is generated from the
722  * same vertex renumbering and fill in a unique interior map.
723  */
724  cnt=0;
725  for(i = 0; i < m_numLocalCoeffs; ++i)
726  {
727  if(m_localToGlobalMap[i] == -1)
728  {
729  m_localToGlobalMap[i] = globalId++;
730  }
731  else
732  {
733  if(m_signChange)
734  {
735  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
736  }
737  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
738  }
739  }
740  m_numGlobalCoeffs = globalId;
741 
742  // Set up the local to global map for the next level when using
743  // multi-level static condensation
744  if( m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond) )
745  {
746  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
747  {
748  Array<OneD, int> vwgts_perm(
749  Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
750  for(i = 0; i < locExpVector.size(); ++i)
751  {
752  locExpansion = locExpVector[i]
754  for(j = 0; j < locExpansion->GetNverts(); ++j)
755  {
756  meshEdgeId = (locExpansion
758  ->GetGeom2D())->GetEid(j);
759  meshVertId = (locExpansion
761  ->GetGeom2D())->GetVid(j);
762 
763  if(ReorderedGraphVertId[0][meshVertId] >=
764  firstNonDirGraphVertId)
765  {
766  vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
767  firstNonDirGraphVertId] =
768  Dofs[0][meshVertId];
769  }
770 
771  if(ReorderedGraphVertId[1][meshEdgeId] >=
772  firstNonDirGraphVertId)
773  {
774  vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
775  firstNonDirGraphVertId] =
776  Dofs[1][meshEdgeId];
777  }
778  }
779  }
780 
781  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
782 
784  AllocateSharedPtr(this,bottomUpGraph);
785  }
786  }
787  }
788 
789 
790 
791 void CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure(vector<map<int,int> > &ReorderedGraphVertId,
792  int &nel, const LocalRegions::ExpansionVector &locExpVector,
793  int &edgeId, int &vertId, int &firstNonDirGraphVertId, map<int,int> &IsDirEdgeDof,
795  Array<OneD, int> &AddMeanPressureToEdgeId)
796 {
797 
798  int i,j,k;
799 
800  // Make list of homogeneous graph edges to elmt mappings
801  Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(),2,-1);
802  map<int,int> HomGraphEdgeIdToEdgeId;
803 
804  for(i = 0; i < nel; ++i)
805  {
806  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
807  {
808  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
809  ->GetGeom2D())->GetEid(j);
810 
811  // note second condition stops us using mixed boundary condition
812  if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
813  && (IsDirEdgeDof.count(edgeId) == 0))
814  {
815  HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
816 
817  if(EdgeIdToElmts[edgeId][0] == -1)
818  {
819  EdgeIdToElmts[edgeId][0] = i;
820  }
821  else
822  {
823  EdgeIdToElmts[edgeId][1] = i;
824  }
825  }
826  }
827  }
828 
829 
831 
832  // Start at second to last level and find edge on boundary
833  // to attach element
834  int nlevels = bottomUpGraph->GetNlevels();
835 
836  // determine a default edge to attach pressure modes to
837  // which is part of the inner solve;
838  int defedge = -1;
839 
840  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
841  for(i = 0; i < bndgraphs.size(); ++i)
842  {
843  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
844 
845  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
846  {
847  // find edge in graph vert list
848  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
849  {
850  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
851 
852  if(defedge == -1)
853  {
854  defedge = edgeId;
855  break;
856  }
857  }
858  }
859  if(defedge != -1)
860  {
861  break;
862  }
863  }
864 
865  for(int n = 1; n < nlevels; ++n)
866  {
867  // produce a map with a key that is the element id
868  // that contains which next level patch it belongs to
869  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
870 
871  // Fill next level graph of adjacent elements and their level
872  map<int,int> ElmtInBndry;
873 
874  for(i = 0; i < bndgraphs.size(); ++i)
875  {
876  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
877 
878  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
879  {
880  // find edge in graph vert list
881  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
882  {
883  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
884 
885  if(EdgeIdToElmts[edgeId][0] != -1)
886  {
887  ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
888  }
889  if(EdgeIdToElmts[edgeId][1] != -1)
890  {
891  ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
892  }
893  }
894  }
895  }
896 
897  // Now search interior patches in this level for edges
898  // that share the same element as a boundary edge and
899  // assign this elmt that boundary edge
900  vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
901  for(i = 0; i < intgraphs.size(); ++i)
902  {
903  int GlobIdOffset = intgraphs[i]->GetIdOffset();
904  bool SetEdge = false;
905  int elmtid = 0;
906  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
907  {
908  // Check to see if graph vert is an edge
909  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
910  {
911  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
912 
913  for(k = 0; k < 2; ++k)
914  {
915  // relevant edge id
916  elmtid = EdgeIdToElmts[edgeId][k];
917 
918  if(elmtid != -1)
919  {
920  mapIt = ElmtInBndry.find(elmtid);
921 
922  if(mapIt != ElmtInBndry.end())
923  {
924  // now find a edge in the next level boundary graph
925  int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
926  for(int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
927  {
928  // find edge in graph vert list
929  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
930  {
931  //June 2012: commenting this condition apparently
932  //solved the bug caused by the edge reordering procedure
933 
934  //if(AddMeanPressureToEdgeId[elmtid] == -1)
935  //{
936 
937  //AddMeanPressureToEdgeId[elmtid] = HomGraphEdgeIdToEdgeId[GlobIdOffset1+l];
938  AddMeanPressureToEdgeId[elmtid] = defedge;
939 
940  //}
941  SetEdge = true;
942  break;
943  }
944  }
945  }
946  }
947  }
948  }
949  }
950 
951 
952  // if we have failed to find matching edge in next
953  // level patch boundary then set last found elmt
954  // associated to this interior patch to the
955  // default edget value
956  if(SetEdge == false)
957  {
958  if(elmtid == -1) // find an elmtid in patch
959  {
960  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
961  {
962  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
963  {
964  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
965  for(k = 0; k < 2; ++k)
966  {
967  // relevant edge id
968  elmtid = EdgeIdToElmts[edgeId][k];
969  if(elmtid != -1)
970  {
971  break;
972  }
973  }
974  }
975  if(elmtid != -1)
976  {
977  break;
978  }
979  }
980  }
981  if(AddMeanPressureToEdgeId[elmtid] == -1)
982  {
983  AddMeanPressureToEdgeId[elmtid] = defedge;
984  }
985  }
986  }
987  }
988 
989 }
990 
991 
992 
993 
994 }
995 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:322
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:316
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
int CreateGraph(const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
Principle Modified Functions .
Definition: BasisType.h:49
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
STL namespace.
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion.h:287
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:397
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:270
static const NekDouble kNekZeroTol
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:824
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:390
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:320
Principle Modified Functions .
Definition: BasisType.h:50
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
double NekDouble
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:392
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
Definition: AssemblyMap.h:350
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:354
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:314
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:386
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:352
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:284
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:305
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
void FindEdgeIdToAddMeanPressure(std::vector< std::map< int, int > > &ReorderedGraphVertId, int &nel, const LocalRegions::ExpansionVector &locExpVector, int &edgeId, int &vertId, int &firstNonDirGraphVertId, std::map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)
Constructs mappings for the C0 scalar continuous Galerkin formulation.
Definition: AssemblyMapCG.h:71
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:388