Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 namespace Nektar
45 {
46  /**
47  * This is an vector extension of
48  * MultiRegions::AssemblyMapCG::SetUp2DExpansionC0ContMap related to the
49  * Linearised Navier Stokes problem
50  */
54  const SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions,
55  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
56  const MultiRegions::ExpListSharedPtr &pressure,
57  const int nz_loc,
58  const bool CheckforSingularSys):
59  AssemblyMapCG(pSession)
60  {
61  int i,j,k,n;
62  int cnt = 0,offset=0;
63  int meshVertId;
64  int meshEdgeId;
65  int bndEdgeCnt;
66  int globalId;
67  int nEdgeCoeffs;
68  int nEdgeInteriorCoeffs;
69  int firstNonDirGraphVertId;
70  int nLocBndCondDofs = 0;
71  int nLocDirBndCondDofs = 0;
72  int nExtraDirichlet = 0;
76  StdRegions::Orientation edgeOrient;
77  Array<OneD, unsigned int> edgeInteriorMap;
78  Array<OneD, int> edgeInteriorSign;
79  int nvel = fields.num_elements();
80 
81  const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
82  int eid, id, diff;
83  int nel = fields[0]->GetNumElmts();
84 
85  MultiRegions::PeriodicMap periodicVerts;
86  MultiRegions::PeriodicMap periodicEdges;
87  MultiRegions::PeriodicMap periodicFaces;
88  vector<map<int,int> > ReorderedGraphVertId(3);
90  int staticCondLevel = 0;
91 
92  if(CheckforSingularSys) //all singularity checking by setting flag to true
93  {
94  m_systemSingular = true;
95  }
96  else // Turn off singular checking by setting flag to false
97  {
98  m_systemSingular = false;
99  }
100 
101  /**
102  * STEP 1: Wrap boundary conditions vector in an array
103  * (since routine is set up for multiple fields) and call
104  * the graph re-odering subroutine to obtain the reordered
105  * values
106  */
107 
108  // Obtain any periodic information and allocate default mapping array
109  fields[0]->GetPeriodicEntities(periodicVerts,periodicEdges,periodicFaces);
110 
111 
112  const Array<OneD, const MultiRegions::ExpListSharedPtr> bndCondExp = fields[0]->GetBndCondExpansions();
113 
114  Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> > bndConditionsVec(nvel);
115  for(i = 0; i < nvel; ++i)
116  {
117  bndConditionsVec[i] = fields[i]->GetBndConditions();
118  }
119 
120  map<int,int> IsDirVertDof;
121  map<int,int> IsDirEdgeDof;
123 
125  for(j = 0; j < bndCondExp.num_elements(); ++j)
126  {
127  map<int,int> BndExpVids;
128  // collect unique list of vertex ids for this expansion
129  for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
130  {
131  g = bndCondExp[j]->GetExp(k)->as<LocalRegions::Expansion1D>()
132  ->GetGeom1D();
133  BndExpVids[g->GetVid(0)] = g->GetVid(0);
134  BndExpVids[g->GetVid(1)] = g->GetVid(1);
135  }
136 
137  for(i = 0; i < nvel; ++i)
138  {
139  if(bndConditionsVec[i][j]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
140  {
141  // set number of Dirichlet conditions along edge
142  for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
143  {
144  IsDirEdgeDof[bndCondExp[j]->GetExp(k)
146  ->GetGeom1D()->GetEid()] += 1;
147  }
148 
149 
150  // Set number of Dirichlet conditions at vertices
151  // with a clamp on its maximum value being nvel to
152  // handle corners between expansions
153  for(mapIt = BndExpVids.begin(); mapIt != BndExpVids.end(); mapIt++)
154  {
155  id = IsDirVertDof[mapIt->second]+1;
156  IsDirVertDof[mapIt->second] = (id > nvel)?nvel:id;
157  }
158  }
159  else
160  {
161  // Check to see that edge normals have non-zero
162  // component in this direction since otherwise
163  // also can be singular.
164  /// @TODO: Fix this so that we can extract normals from edges
165  for(k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
166  {
167  Array<OneD,Array<OneD,NekDouble> > locnorm;
169  = bndCondExp[j]->GetExp(k)
171  locnorm = loc_exp->GetLeftAdjacentElementExp()->GetEdgeNormal(loc_exp->GetLeftAdjacentElementEdge());
172  //locnorm = bndCondExp[j]->GetExp(k)->Get GetMetricInfo()->GetNormal();
173 
174  int ndir = locnorm.num_elements();
175  if(i < ndir) // account for Fourier version where n can be larger then ndir
176  {
177  for(int l = 0; l < locnorm[0].num_elements(); ++l)
178  {
179  if(fabs(locnorm[i][l]) > NekConstants::kNekZeroTol)
180  {
181  m_systemSingular = false;
182  break;
183  }
184  }
185  }
186  if(m_systemSingular == false)
187  {
188  break;
189  }
190  }
191  }
192  }
193  }
194 
195  Array<OneD, map<int,int> >Dofs(2);
196 
197  Array<OneD, int> AddMeanPressureToEdgeId(nel,-1);
198  int edgeId,vertId;
199 
200 
201  // special case of singular problem - need to fix one pressure
202  // dof to a dirichlet edge. Since we attached pressure dof to
203  // last velocity component of edge need to make sure this
204  // component is Dirichlet
205  if(m_systemSingular)
206  {
207  id = -1;
208  for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
209  {
210  if(bndConditionsVec[nvel-1][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
211  {
212  id = bndCondExp[i]->GetExp(0)
213  ->as<LocalRegions::Expansion1D>()->GetGeom1D()
214  ->GetEid();
215  break;
216  }
217  }
218 
219  ASSERTL0(id != -1," Did not find an edge to attach singular pressure degree of freedom");
220 
221  // determine element with this edge id. There may be a
222  // more direct way of getting element from spatialDomains
223  for(i = 0; i < nel; ++i)
224  {
225  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
226  {
227  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
228  ->GetGeom2D())->GetEid(j);
229 
230  if(edgeId == id)
231  {
232  AddMeanPressureToEdgeId[i] = id;
233  break;
234  }
235  }
236 
237  if(AddMeanPressureToEdgeId[i] != -1)
238  {
239  break;
240  }
241  }
242  }
243 
244 
245  for(i = 0; i < nel; ++i)
246  {
247  eid = fields[0]->GetOffset_Elmt_Id(i);
248  for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
249  {
250  vertId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
251  ->GetGeom2D())->GetVid(j);
252  if(Dofs[0].count(vertId) == 0)
253  {
254  Dofs[0][vertId] = nvel*nz_loc;
255 
256  // Adjust for a Dirichlet boundary condition to give number to be solved
257  if(IsDirVertDof.count(vertId) != 0)
258  {
259  Dofs[0][vertId] -= IsDirVertDof[vertId]*nz_loc;
260  }
261  }
262 
263  edgeId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
264  ->GetGeom2D())->GetEid(j);
265  if(Dofs[1].count(edgeId) == 0)
266  {
267  Dofs[1][edgeId] = nvel*(locExpVector[eid]->GetEdgeNcoeffs(j)-2)*nz_loc;
268  }
269 
270  // Adjust for Dirichlet boundary conditions to give number to be solved
271  if(IsDirEdgeDof.count(edgeId) != 0)
272  {
273  Dofs[1][edgeId] -= IsDirEdgeDof[edgeId]*nz_loc*(locExpVector[eid]->GetEdgeNcoeffs(j)-2);
274  }
275  }
276  }
277 
278  set<int> extraDirVerts, extraDirEdges;
279 
280  CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false,
281  periodicVerts, periodicEdges, periodicFaces,
282  ReorderedGraphVertId, bottomUpGraph, extraDirVerts,
283  extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, 4);
284  /*
285  SetUp2DGraphC0ContMap(*fields[0],
286  bndCondExp,
287  bndConditionsVec,
288  periodicVerts, periodicEdges,
289  Dofs, ReorderedGraphVertId,
290  firstNonDirGraphVertId, nExtraDirichlet,
291  bottomUpGraph, extraDir, false, 4);
292  */
293 
294  /**
295  * STEP 2a: Set the mean pressure modes to edges depending on
296  * type of direct solver technique;
297  */
298 
299  // determine which edge to add mean pressure dof based on
300  // ensuring that at least one pressure dof from an internal
301  // patch is associated with its boundary system
302  if(m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
303  {
304 
305 
306  FindEdgeIdToAddMeanPressure(ReorderedGraphVertId,
307  nel, locExpVector,
308  edgeId, vertId, firstNonDirGraphVertId, IsDirEdgeDof,
309  bottomUpGraph,
310  AddMeanPressureToEdgeId);
311  }
312 
313  // Set unset elmts to non-Dirichlet edges.
314  // special case of singular problem - need to fix one
315  // pressure dof to a dirichlet edge
316  for(i = 0; i < nel; ++i)
317  {
318  eid = fields[0]->GetOffset_Elmt_Id(i);
319  for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
320  {
321  edgeId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
322  ->GetGeom2D())->GetEid(j);
323 
324  if(IsDirEdgeDof.count(edgeId) == 0) // interior edge
325  {
326  // setup AddMeanPressureToEdgeId to decide where to
327  // put pressure
328  if(AddMeanPressureToEdgeId[eid] == -1)
329  {
330  AddMeanPressureToEdgeId[eid] = edgeId;
331  }
332  }
333  }
334  ASSERTL0((AddMeanPressureToEdgeId[eid] != -1),"Did not determine "
335  "an edge to attach mean pressure dof");
336  // Add the mean pressure degree of freedom to this edge
337  Dofs[1][AddMeanPressureToEdgeId[eid]] += nz_loc;
338  }
339 
340  map<int,int> pressureEdgeOffset;
341 
342  /**
343  * STEP 2: Count out the number of Dirichlet vertices and edges first
344  */
345  for(i = 0; i < bndCondExp.num_elements(); i++)
346  {
347  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
348  {
349  bndSegExp = bndCondExp[i]->GetExp(j)
350  ->as<LocalRegions::SegExp>();
351  for(k = 0; k < nvel; ++k)
352  {
353  if(bndConditionsVec[k][i]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
354  {
355  nLocDirBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
356  }
357  nLocBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
358  }
359  }
360  }
361 
362  if(m_systemSingular)
363  {
364  m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet+nz_loc;
365  }
366  else
367  {
368  m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet;
369  }
370 
371  /**
372  * STEP 3: Set up an array which contains the offset information of
373  * the different graph vertices.
374  *
375  * This basically means to identify how many global degrees of
376  * freedom the individual graph vertices correspond. Obviously,
377  * the graph vertices corresponding to the mesh-vertices account
378  * for a single global DOF. However, the graph vertices
379  * corresponding to the element edges correspond to 2*(N-2) global DOF
380  * where N is equal to the number of boundary modes on this edge.
381  */
382  Array<OneD, int> graphVertOffset(nvel*nz_loc*(ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),0);
383  graphVertOffset[0] = 0;
384 
385  m_signChange = false;
386 
387  for(i = 0; i < nel; ++i)
388  {
389  eid = fields[0]->GetOffset_Elmt_Id(i);
390  locExpansion = locExpVector[eid]->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 
556  m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
557  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
558  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int>(nLocBndCondDofs,-1);
559 
560 
561  // Set default sign array.
562  m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
563  m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
564 
565  m_staticCondLevel = staticCondLevel;
566  m_numPatches = nel;
567 
568  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
569  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);
570 
571  for(i = 0; i < nel; ++i)
572  {
573  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) nz_loc*(nvel*locExpVector[fields[0]->GetOffset_Elmt_Id(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;
588  Array<OneD, unsigned int> bmap;
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  eid = fields[0]->GetOffset_Elmt_Id(i);
596  locExpansion = locExpVector[eid]->as<StdRegions::StdExpansion2D>();
597 
598  velnbndry = locExpansion->NumBndryCoeffs();
599 
600  // require an inverse ordering of the bmap system to store
601  // local numbering system which takes matrix these
602  // matrices. Therefore get hold of elemental bmap and set
603  // up an inverse map
604  map<int,int> inv_bmap;
605  locExpansion->GetBoundaryMap(bmap);
606  for(j = 0; j < bmap.num_elements(); ++j)
607  {
608  inv_bmap[bmap[j]] = j;
609  }
610 
611  // Loop over all edges (and vertices) of element i
612  for(j = 0; j < locExpansion->GetNedges(); ++j)
613  {
614  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
615  edgeOrient = (locExpansion->as<LocalRegions::Expansion2D>()
616  ->GetGeom2D())->GetEorient(j);
617  meshEdgeId = (locExpansion->as<LocalRegions::Expansion2D>()
618  ->GetGeom2D())->GetEid(j);
619  meshVertId = (locExpansion->as<LocalRegions::Expansion2D>()
620  ->GetGeom2D())->GetVid(j);
621 
622  locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
623  // Set the global DOF for vertex j of element i
624 
625  for(nv = 0; nv < nvel*nz_loc; ++nv)
626  {
627  m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[locExpansion->GetVertexMap(j)]] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+ nv];
628  // Set the global DOF's for the interior modes of edge j
629  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
630  {
631  m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[edgeInteriorMap[k]]] = graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv]+k;
632  }
633  }
634 
635  // Fill the sign vector if required
636  if(m_signChange)
637  {
638  for(nv = 0; nv < nvel*nz_loc; ++nv)
639  {
640  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
641  {
642  m_localToGlobalSign[cnt+nv*velnbndry + inv_bmap[edgeInteriorMap[k]]] = (NekDouble) edgeInteriorSign[k];
643  }
644  }
645  }
646  }
647 
648  // use difference between two edges of the AddMeanPressureEdgeId to det nEdgeInteriorCoeffs.
649  nEdgeInteriorCoeffs = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc+1] - graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc];
650 
651  int psize = pressure->GetExp(eid)->GetNcoeffs();
652  for(n = 0; n < nz_loc; ++n)
653  {
654  m_localToGlobalMap[cnt + nz_loc*nvel*velnbndry + n*psize] = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]]+1)*nvel*nz_loc-1]+nEdgeInteriorCoeffs + pressureEdgeOffset[AddMeanPressureToEdgeId[eid]];
655 
656  pressureEdgeOffset[AddMeanPressureToEdgeId[eid]] += 1;
657  }
658 
659  cnt += (velnbndry*nvel+ psize)*nz_loc;
660  }
661 
662  // Set up the mapping for the boundary conditions
663  offset = cnt = 0;
664  for(nv = 0; nv < nvel; ++nv)
665  {
666  for(i = 0; i < bndCondExp.num_elements(); i++)
667  {
668  for(n = 0; n < nz_loc; ++n)
669  {
670  int ncoeffcnt = 0;
671  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
672  {
673  bndSegExp = bndCondExp[i]->GetExp(j)
674  ->as<LocalRegions::SegExp>();
675 
676  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
677  for(k = 0; k < 2; k++)
678  {
679  meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
680  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+nv*nz_loc+n];
681  }
682 
683  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
684  bndEdgeCnt = 0;
685  nEdgeCoeffs = bndSegExp->GetNcoeffs();
686  for(k = 0; k < nEdgeCoeffs; k++)
687  {
688  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
689  {
690  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
691  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv*nz_loc+n]+bndEdgeCnt;
692  bndEdgeCnt++;
693  }
694  }
695  ncoeffcnt += nEdgeCoeffs;
696  }
697  // Note: Can not use bndCondExp[i]->GetNcoeffs()
698  // due to homogeneous extension not returning just
699  // the value per plane
700  offset += ncoeffcnt;
701  }
702  }
703  }
704 
705  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
706  m_numGlobalBndCoeffs = globalId;
707 
708  /**
709  * STEP 5: The boundary condition mapping is generated from the
710  * same vertex renumbering and fill in a unique interior map.
711  */
712  cnt=0;
713  for(i = 0; i < m_numLocalCoeffs; ++i)
714  {
715  if(m_localToGlobalMap[i] == -1)
716  {
717  m_localToGlobalMap[i] = globalId++;
718  }
719  else
720  {
721  if(m_signChange)
722  {
723  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
724  }
725  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
726  }
727  }
728  m_numGlobalCoeffs = globalId;
729 
730  // Set up the local to global map for the next level when using
731  // multi-level static condensation
732  if( m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond) )
733  {
734  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
735  {
736  Array<OneD, int> vwgts_perm(
737  Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
738  for(i = 0; i < locExpVector.size(); ++i)
739  {
740  eid = fields[0]->GetOffset_Elmt_Id(i);
741  locExpansion = locExpVector[eid]
743  for(j = 0; j < locExpansion->GetNverts(); ++j)
744  {
745  meshEdgeId = (locExpansion
747  ->GetGeom2D())->GetEid(j);
748  meshVertId = (locExpansion
750  ->GetGeom2D())->GetVid(j);
751 
752  if(ReorderedGraphVertId[0][meshVertId] >=
753  firstNonDirGraphVertId)
754  {
755  vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
756  firstNonDirGraphVertId] =
757  Dofs[0][meshVertId];
758  }
759 
760  if(ReorderedGraphVertId[1][meshEdgeId] >=
761  firstNonDirGraphVertId)
762  {
763  vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
764  firstNonDirGraphVertId] =
765  Dofs[1][meshEdgeId];
766  }
767  }
768  }
769 
770  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
771 
773  AllocateSharedPtr(this,bottomUpGraph);
774  }
775  }
776  }
777 }
778 
779 
780 void CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure(vector<map<int,int> > &ReorderedGraphVertId,
781  int &nel, const LocalRegions::ExpansionVector &locExpVector,
782  int &edgeId, int &vertId, int &firstNonDirGraphVertId, map<int,int> &IsDirEdgeDof,
784  Array<OneD, int> &AddMeanPressureToEdgeId)
785 {
786 
787  int i,j,k;
788 
789  // Make list of homogeneous graph edges to elmt mappings
790  Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(),2,-1);
791  map<int,int> HomGraphEdgeIdToEdgeId;
792 
793  for(i = 0; i < nel; ++i)
794  {
795  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
796  {
797  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
798  ->GetGeom2D())->GetEid(j);
799 
800  // note second condition stops us using mixed boundary condition
801  if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
802  && (IsDirEdgeDof.count(edgeId) == 0))
803  {
804  HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
805 
806  if(EdgeIdToElmts[edgeId][0] == -1)
807  {
808  EdgeIdToElmts[edgeId][0] = i;
809  }
810  else
811  {
812  EdgeIdToElmts[edgeId][1] = i;
813  }
814  }
815  }
816  }
817 
818 
820 
821  // Start at second to last level and find edge on boundary
822  // to attach element
823  int nlevels = bottomUpGraph->GetNlevels();
824 
825  // determine a default edge to attach pressure modes to
826  // which is part of the inner solve;
827  int defedge = -1;
828 
829  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
830  for(i = 0; i < bndgraphs.size(); ++i)
831  {
832  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
833 
834  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
835  {
836  // find edge in graph vert list
837  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
838  {
839  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
840 
841  if(defedge == -1)
842  {
843  defedge = edgeId;
844  break;
845  }
846  }
847  }
848  if(defedge != -1)
849  {
850  break;
851  }
852  }
853 
854  for(int n = 1; n < nlevels; ++n)
855  {
856  // produce a map with a key that is the element id
857  // that contains which next level patch it belongs to
858  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
859 
860  // Fill next level graph of adjacent elements and their level
861  map<int,int> ElmtInBndry;
862 
863  for(i = 0; i < bndgraphs.size(); ++i)
864  {
865  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
866 
867  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
868  {
869  // find edge in graph vert list
870  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
871  {
872  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
873 
874  if(EdgeIdToElmts[edgeId][0] != -1)
875  {
876  ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
877  }
878  if(EdgeIdToElmts[edgeId][1] != -1)
879  {
880  ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
881  }
882  }
883  }
884  }
885 
886  // Now search interior patches in this level for edges
887  // that share the same element as a boundary edge and
888  // assign this elmt that boundary edge
889  vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
890  for(i = 0; i < intgraphs.size(); ++i)
891  {
892  int GlobIdOffset = intgraphs[i]->GetIdOffset();
893  bool SetEdge = false;
894  int elmtid = 0;
895  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
896  {
897  // Check to see if graph vert is an edge
898  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
899  {
900  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
901 
902  for(k = 0; k < 2; ++k)
903  {
904  // relevant edge id
905  elmtid = EdgeIdToElmts[edgeId][k];
906 
907  if(elmtid != -1)
908  {
909  mapIt = ElmtInBndry.find(elmtid);
910 
911  if(mapIt != ElmtInBndry.end())
912  {
913  // now find a edge in the next level boundary graph
914  int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
915  for(int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
916  {
917  // find edge in graph vert list
918  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
919  {
920  //June 2012: commenting this condition apparently
921  //solved the bug caused by the edge reordering procedure
922 
923  //if(AddMeanPressureToEdgeId[elmtid] == -1)
924  //{
925 
926  //AddMeanPressureToEdgeId[elmtid] = HomGraphEdgeIdToEdgeId[GlobIdOffset1+l];
927  AddMeanPressureToEdgeId[elmtid] = defedge;
928 
929  //}
930  SetEdge = true;
931  break;
932  }
933  }
934  }
935  }
936  }
937  }
938  }
939 
940 
941  // if we have failed to find matching edge in next
942  // level patch boundary then set last found elmt
943  // associated to this interior patch to the
944  // default edget value
945  if(SetEdge == false)
946  {
947  if(elmtid == -1) // find an elmtid in patch
948  {
949  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
950  {
951  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
952  {
953  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
954  for(k = 0; k < 2; ++k)
955  {
956  // relevant edge id
957  elmtid = EdgeIdToElmts[edgeId][k];
958  if(elmtid != -1)
959  {
960  break;
961  }
962  }
963  }
964  if(elmtid != -1)
965  {
966  break;
967  }
968  }
969  }
970  if(AddMeanPressureToEdgeId[elmtid] == -1)
971  {
972  AddMeanPressureToEdgeId[elmtid] = defedge;
973  }
974  }
975  }
976  }
977 
978 }
979 
980 
981 
982 
983 
984