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