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 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 eid, 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  eid = fields[0]->GetOffset_Elmt_Id(i);
251  for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
252  {
253  vertId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
254  ->GetGeom2D())->GetVid(j);
255  if(Dofs[0].count(vertId) == 0)
256  {
257  Dofs[0][vertId] = nvel*nz_loc;
258 
259  // Adjust for a Dirichlet boundary condition to give number to be solved
260  if(IsDirVertDof.count(vertId) != 0)
261  {
262  Dofs[0][vertId] -= IsDirVertDof[vertId]*nz_loc;
263  }
264  }
265 
266  edgeId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
267  ->GetGeom2D())->GetEid(j);
268  if(Dofs[1].count(edgeId) == 0)
269  {
270  Dofs[1][edgeId] = nvel*(locExpVector[eid]->GetEdgeNcoeffs(j)-2)*nz_loc;
271  }
272 
273  // Adjust for Dirichlet boundary conditions to give number to be solved
274  if(IsDirEdgeDof.count(edgeId) != 0)
275  {
276  Dofs[1][edgeId] -= IsDirEdgeDof[edgeId]*nz_loc*(locExpVector[eid]->GetEdgeNcoeffs(j)-2);
277  }
278  }
279  }
280 
281  set<int> extraDirVerts, extraDirEdges;
282 
283  CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false,
284  periodicVerts, periodicEdges, periodicFaces,
285  ReorderedGraphVertId, bottomUpGraph, extraDirVerts,
286  extraDirEdges, firstNonDirGraphVertId, nExtraDirichlet, 4);
287  /*
288  SetUp2DGraphC0ContMap(*fields[0],
289  bndCondExp,
290  bndConditionsVec,
291  periodicVerts, periodicEdges,
292  Dofs, ReorderedGraphVertId,
293  firstNonDirGraphVertId, nExtraDirichlet,
294  bottomUpGraph, extraDir, false, 4);
295  */
296 
297  /**
298  * STEP 2a: Set the mean pressure modes to edges depending on
299  * type of direct solver technique;
300  */
301 
302  // determine which edge to add mean pressure dof based on
303  // ensuring that at least one pressure dof from an internal
304  // patch is associated with its boundary system
305  if(m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond))
306  {
307 
308 
309  FindEdgeIdToAddMeanPressure(ReorderedGraphVertId,
310  nel, locExpVector,
311  edgeId, vertId, firstNonDirGraphVertId, IsDirEdgeDof,
312  bottomUpGraph,
313  AddMeanPressureToEdgeId);
314  }
315 
316  // Set unset elmts to non-Dirichlet edges.
317  // special case of singular problem - need to fix one
318  // pressure dof to a dirichlet edge
319  for(i = 0; i < nel; ++i)
320  {
321  eid = fields[0]->GetOffset_Elmt_Id(i);
322  for(j = 0; j < locExpVector[eid]->GetNverts(); ++j)
323  {
324  edgeId = (locExpVector[eid]->as<LocalRegions::Expansion2D>()
325  ->GetGeom2D())->GetEid(j);
326 
327  if(IsDirEdgeDof.count(edgeId) == 0) // interior edge
328  {
329  // setup AddMeanPressureToEdgeId to decide where to
330  // put pressure
331  if(AddMeanPressureToEdgeId[eid] == -1)
332  {
333  AddMeanPressureToEdgeId[eid] = edgeId;
334  }
335  }
336  }
337  ASSERTL0((AddMeanPressureToEdgeId[eid] != -1),"Did not determine "
338  "an edge to attach mean pressure dof");
339  // Add the mean pressure degree of freedom to this edge
340  Dofs[1][AddMeanPressureToEdgeId[eid]] += nz_loc;
341  }
342 
343  map<int,int> pressureEdgeOffset;
344 
345  /**
346  * STEP 2: Count out the number of Dirichlet vertices and edges first
347  */
348  for(i = 0; i < bndCondExp.num_elements(); i++)
349  {
350  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
351  {
352  bndSegExp = bndCondExp[i]->GetExp(j)
353  ->as<LocalRegions::SegExp>();
354  for(k = 0; k < nvel; ++k)
355  {
356  if(bndConditionsVec[k][i]->GetBoundaryConditionType()==SpatialDomains::eDirichlet)
357  {
358  nLocDirBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
359  }
360  nLocBndCondDofs += bndSegExp->GetNcoeffs()*nz_loc;
361  }
362  }
363  }
364 
365  if(m_systemSingular)
366  {
367  m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet+nz_loc;
368  }
369  else
370  {
371  m_numLocalDirBndCoeffs = nLocDirBndCondDofs+nExtraDirichlet;
372  }
373 
374  /**
375  * STEP 3: Set up an array which contains the offset information of
376  * the different graph vertices.
377  *
378  * This basically means to identify how many global degrees of
379  * freedom the individual graph vertices correspond. Obviously,
380  * the graph vertices corresponding to the mesh-vertices account
381  * for a single global DOF. However, the graph vertices
382  * corresponding to the element edges correspond to 2*(N-2) global DOF
383  * where N is equal to the number of boundary modes on this edge.
384  */
385  Array<OneD, int> graphVertOffset(nvel*nz_loc*(ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),0);
386  graphVertOffset[0] = 0;
387 
388  m_signChange = false;
389 
390  for(i = 0; i < nel; ++i)
391  {
392  eid = fields[0]->GetOffset_Elmt_Id(i);
393  locExpansion = locExpVector[eid]->as<StdRegions::StdExpansion2D>();
394 
395  for(j = 0; j < locExpansion->GetNedges(); ++j)
396  {
397  nEdgeCoeffs = locExpansion->GetEdgeNcoeffs(j);
398  meshEdgeId = (locExpansion->as<LocalRegions::Expansion2D>()
399  ->GetGeom2D())->GetEid(j);
400  meshVertId = (locExpansion->as<LocalRegions::Expansion2D>()
401  ->GetGeom2D())->GetVid(j);
402 
403  for(k = 0; k < nvel*nz_loc; ++k)
404  {
405  graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+k] = 1;
406  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+k] = (nEdgeCoeffs-2);
407  }
408 
409  bType = locExpansion->GetEdgeBasisType(j);
410  // need a sign vector for modal expansions if nEdgeCoeffs >=4
411  if( (nEdgeCoeffs >= 4)&&
412  ( (bType == LibUtilities::eModified_A)||
413  (bType == LibUtilities::eModified_B) ) )
414  {
415  m_signChange = true;
416  }
417  }
418  }
419 
420  // Add mean pressure modes;
421  for(i = 0; i < nel; ++i)
422  {
423  graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]]+1)*nvel*nz_loc-1] += nz_loc;
424  //graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc] += nz_loc;
425  }
426 
427  // Negate the vertices and edges with only a partial
428  // Dirichlet conditon. Essentially we check to see if an edge
429  // has a mixed Dirichlet with Neumann/Robin Condition and if
430  // so negate the offset associated with this vertex.
431 
432  map<int,int> DirVertChk;
433 
434  for(i = 0; i < bndConditionsVec[0].num_elements(); ++i)
435  {
436  cnt = 0;
437  for(j = 0; j < nvel; ++j)
438  {
439  if(bndConditionsVec[j][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
440  {
441  cnt ++;
442  }
443  }
444 
445  // Case where partial Dirichlet boundary condition
446  if((cnt > 0)&&(cnt < nvel))
447  {
448  for(j = 0; j < nvel; ++j)
449  {
450  if(bndConditionsVec[j][i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
451  {
452  //negate graph offsets which should be
453  //Dirichlet conditions
454  for(k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
455  {
456  // vertices with mix condition;
457  id = bndCondExp[i]->GetExp(k)
459  ->GetGeom1D()->GetVid(0);
460  if(DirVertChk.count(id*nvel+j) == 0)
461  {
462  DirVertChk[id*nvel+j] = 1;
463  for(n = 0; n < nz_loc; ++n)
464  {
465  graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc + n] *= -1;
466  }
467  }
468 
469  id = bndCondExp[i]->GetExp(k)
471  ->GetGeom1D()->GetVid(1);
472  if(DirVertChk.count(id*nvel+j) == 0)
473  {
474  DirVertChk[id*nvel+j] = 1;
475  for(n = 0; n < nz_loc; ++n)
476  {
477  graphVertOffset[ReorderedGraphVertId[0][id]*nvel*nz_loc+j*nz_loc+n] *= -1;
478  }
479  }
480 
481  // edges with mixed id;
482  id = bndCondExp[i]->GetExp(k)
484  ->GetGeom1D()->GetEid();
485  for(n = 0; n < nz_loc; ++n)
486  {
487  graphVertOffset[ReorderedGraphVertId[1][id]*nvel*nz_loc+j*nz_loc +n] *= -1;
488  }
489  }
490  }
491  }
492  }
493  }
494 
495 
496  cnt = 0;
497  // assemble accumulative list of full Dirichlet values.
498  for(i = 0; i < firstNonDirGraphVertId*nvel*nz_loc; ++i)
499  {
500  diff = abs(graphVertOffset[i]);
501  graphVertOffset[i] = cnt;
502  cnt += diff;
503  }
504 
505  // set Dirichlet values with negative values to Dirichlet value
506  for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
507  {
508  if(graphVertOffset[i] < 0)
509  {
510  diff = -graphVertOffset[i];
511  graphVertOffset[i] = -cnt;
512  cnt += diff;
513  }
514  }
515 
516  // Accumulate all interior degrees of freedom with positive values
518 
519  // offset values
520  for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
521  {
522  if(graphVertOffset[i] >= 0)
523  {
524  diff = graphVertOffset[i];
525  graphVertOffset[i] = cnt;
526  cnt += diff;
527  }
528  }
529 
530  // Finally set negative entries (corresponding to Dirichlet
531  // values ) to be positive
532  for(i = firstNonDirGraphVertId*nvel*nz_loc; i < graphVertOffset.num_elements(); ++i)
533  {
534  if(graphVertOffset[i] < 0)
535  {
536  graphVertOffset[i] = -graphVertOffset[i];
537  }
538  }
539 
540 
541  // Allocate the proper amount of space for the class-data and fill
542  // information that is already known
543  cnt = 0;
545  m_numLocalCoeffs = 0;
546 
547  for(i = 0; i < nel; ++i)
548  {
549  m_numLocalBndCoeffs += nz_loc*(nvel*locExpVector[i]->NumBndryCoeffs() + 1);
550  // add these coeffs up separately since
551  // pressure->GetNcoeffs can include the coefficient in
552  // multiple planes.
553  m_numLocalCoeffs += (pressure->GetExp(i)->GetNcoeffs()-1)*nz_loc;
554  }
555 
557 
558 
562 
563 
564  // Set default sign array.
567 
568  m_staticCondLevel = staticCondLevel;
569  m_numPatches = nel;
570 
573 
574  for(i = 0; i < nel; ++i)
575  {
576  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) nz_loc*(nvel*locExpVector[fields[0]->GetOffset_Elmt_Id(i)]->NumBndryCoeffs() + 1);
577  m_numLocalIntCoeffsPerPatch[i] = (unsigned int) nz_loc*(pressure->GetExp(i)->GetNcoeffs()-1);
578  }
579 
580  /**
581  * STEP 4: Now, all ingredients are ready to set up the actual
582  * local to global mapping.
583  *
584  * The remainder of the map consists of the element-interior
585  * degrees of freedom. This leads to the block-diagonal submatrix
586  * as each element-interior mode is globally orthogonal to modes
587  * in all other elements.
588  */
589  cnt = 0;
590  int nv,velnbndry;
592 
593 
594  // Loop over all the elements in the domain in shuffled
595  // ordering (element type consistency)
596  for(i = 0; i < nel; ++i)
597  {
598  eid = fields[0]->GetOffset_Elmt_Id(i);
599  locExpansion = locExpVector[eid]->as<StdRegions::StdExpansion2D>();
600 
601  velnbndry = locExpansion->NumBndryCoeffs();
602 
603  // require an inverse ordering of the bmap system to store
604  // local numbering system which takes matrix these
605  // matrices. Therefore get hold of elemental bmap and set
606  // up an inverse map
607  map<int,int> inv_bmap;
608  locExpansion->GetBoundaryMap(bmap);
609  for(j = 0; j < bmap.num_elements(); ++j)
610  {
611  inv_bmap[bmap[j]] = j;
612  }
613 
614  // Loop over all edges (and vertices) of element i
615  for(j = 0; j < locExpansion->GetNedges(); ++j)
616  {
617  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
618  edgeOrient = (locExpansion->as<LocalRegions::Expansion2D>()
619  ->GetGeom2D())->GetEorient(j);
620  meshEdgeId = (locExpansion->as<LocalRegions::Expansion2D>()
621  ->GetGeom2D())->GetEid(j);
622  meshVertId = (locExpansion->as<LocalRegions::Expansion2D>()
623  ->GetGeom2D())->GetVid(j);
624 
625  pIt = periodicEdges.find(meshEdgeId);
626 
627  // See if this edge is periodic. If it is, then we map all
628  // edges to the one with lowest ID, and align all
629  // coefficients to this edge orientation.
630  if (pIt != periodicEdges.end())
631  {
632  pair<int, StdRegions::Orientation> idOrient =
634  meshEdgeId, edgeOrient, pIt->second);
635  edgeOrient = idOrient.second;
636  }
637 
638  locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
639  // Set the global DOF for vertex j of element i
640 
641  for(nv = 0; nv < nvel*nz_loc; ++nv)
642  {
643  m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[locExpansion->GetVertexMap(j)]] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+ nv];
644  // Set the global DOF's for the interior modes of edge j
645  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
646  {
647  m_localToGlobalMap[cnt+nv*velnbndry+inv_bmap[edgeInteriorMap[k]]] = graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv]+k;
648  }
649  }
650 
651  // Fill the sign vector if required
652  if(m_signChange)
653  {
654  for(nv = 0; nv < nvel*nz_loc; ++nv)
655  {
656  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
657  {
658  m_localToGlobalSign[cnt+nv*velnbndry + inv_bmap[edgeInteriorMap[k]]] = (NekDouble) edgeInteriorSign[k];
659  }
660  }
661  }
662  }
663 
664  // use difference between two edges of the AddMeanPressureEdgeId to det nEdgeInteriorCoeffs.
665  nEdgeInteriorCoeffs = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc+1] - graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]])*nvel*nz_loc];
666 
667  int psize = pressure->GetExp(eid)->GetNcoeffs();
668  for(n = 0; n < nz_loc; ++n)
669  {
670  m_localToGlobalMap[cnt + nz_loc*nvel*velnbndry + n*psize] = graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[eid]]+1)*nvel*nz_loc-1]+nEdgeInteriorCoeffs + pressureEdgeOffset[AddMeanPressureToEdgeId[eid]];
671 
672  pressureEdgeOffset[AddMeanPressureToEdgeId[eid]] += 1;
673  }
674 
675  cnt += (velnbndry*nvel+ psize)*nz_loc;
676  }
677 
678  // Set up the mapping for the boundary conditions
679  offset = cnt = 0;
680  for(nv = 0; nv < nvel; ++nv)
681  {
682  for(i = 0; i < bndCondExp.num_elements(); i++)
683  {
684  for(n = 0; n < nz_loc; ++n)
685  {
686  int ncoeffcnt = 0;
687  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
688  {
689  bndSegExp = bndCondExp[i]->GetExp(j)
690  ->as<LocalRegions::SegExp>();
691 
692  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
693  for(k = 0; k < 2; k++)
694  {
695  meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
696  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndSegExp->GetVertexMap(k)] = graphVertOffset[ReorderedGraphVertId[0][meshVertId]*nvel*nz_loc+nv*nz_loc+n];
697  }
698 
699  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
700  bndEdgeCnt = 0;
701  nEdgeCoeffs = bndSegExp->GetNcoeffs();
702  for(k = 0; k < nEdgeCoeffs; k++)
703  {
704  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
705  {
706  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
707  graphVertOffset[ReorderedGraphVertId[1][meshEdgeId]*nvel*nz_loc+nv*nz_loc+n]+bndEdgeCnt;
708  bndEdgeCnt++;
709  }
710  }
711  ncoeffcnt += nEdgeCoeffs;
712  }
713  // Note: Can not use bndCondExp[i]->GetNcoeffs()
714  // due to homogeneous extension not returning just
715  // the value per plane
716  offset += ncoeffcnt;
717  }
718  }
719  }
720 
721  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
722  m_numGlobalBndCoeffs = globalId;
723 
724  /**
725  * STEP 5: The boundary condition mapping is generated from the
726  * same vertex renumbering and fill in a unique interior map.
727  */
728  cnt=0;
729  for(i = 0; i < m_numLocalCoeffs; ++i)
730  {
731  if(m_localToGlobalMap[i] == -1)
732  {
733  m_localToGlobalMap[i] = globalId++;
734  }
735  else
736  {
737  if(m_signChange)
738  {
739  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
740  }
741  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
742  }
743  }
744  m_numGlobalCoeffs = globalId;
745 
746  // Set up the local to global map for the next level when using
747  // multi-level static condensation
748  if( m_session->MatchSolverInfoAsEnum("GlobalSysSoln", MultiRegions::eDirectMultiLevelStaticCond) )
749  {
750  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
751  {
752  Array<OneD, int> vwgts_perm(
753  Dofs[0].size()+Dofs[1].size()-firstNonDirGraphVertId);
754  for(i = 0; i < locExpVector.size(); ++i)
755  {
756  eid = fields[0]->GetOffset_Elmt_Id(i);
757  locExpansion = locExpVector[eid]
759  for(j = 0; j < locExpansion->GetNverts(); ++j)
760  {
761  meshEdgeId = (locExpansion
763  ->GetGeom2D())->GetEid(j);
764  meshVertId = (locExpansion
766  ->GetGeom2D())->GetVid(j);
767 
768  if(ReorderedGraphVertId[0][meshVertId] >=
769  firstNonDirGraphVertId)
770  {
771  vwgts_perm[ReorderedGraphVertId[0][meshVertId]-
772  firstNonDirGraphVertId] =
773  Dofs[0][meshVertId];
774  }
775 
776  if(ReorderedGraphVertId[1][meshEdgeId] >=
777  firstNonDirGraphVertId)
778  {
779  vwgts_perm[ReorderedGraphVertId[1][meshEdgeId]-
780  firstNonDirGraphVertId] =
781  Dofs[1][meshEdgeId];
782  }
783  }
784  }
785 
786  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
787 
789  AllocateSharedPtr(this,bottomUpGraph);
790  }
791  }
792  }
793 
794 
795 
796 void CoupledLocalToGlobalC0ContMap::FindEdgeIdToAddMeanPressure(vector<map<int,int> > &ReorderedGraphVertId,
797  int &nel, const LocalRegions::ExpansionVector &locExpVector,
798  int &edgeId, int &vertId, int &firstNonDirGraphVertId, map<int,int> &IsDirEdgeDof,
800  Array<OneD, int> &AddMeanPressureToEdgeId)
801 {
802 
803  int i,j,k;
804 
805  // Make list of homogeneous graph edges to elmt mappings
806  Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(),2,-1);
807  map<int,int> HomGraphEdgeIdToEdgeId;
808 
809  for(i = 0; i < nel; ++i)
810  {
811  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
812  {
813  edgeId = (locExpVector[i]->as<LocalRegions::Expansion2D>()
814  ->GetGeom2D())->GetEid(j);
815 
816  // note second condition stops us using mixed boundary condition
817  if((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId)
818  && (IsDirEdgeDof.count(edgeId) == 0))
819  {
820  HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId]-firstNonDirGraphVertId] = edgeId;
821 
822  if(EdgeIdToElmts[edgeId][0] == -1)
823  {
824  EdgeIdToElmts[edgeId][0] = i;
825  }
826  else
827  {
828  EdgeIdToElmts[edgeId][1] = i;
829  }
830  }
831  }
832  }
833 
834 
836 
837  // Start at second to last level and find edge on boundary
838  // to attach element
839  int nlevels = bottomUpGraph->GetNlevels();
840 
841  // determine a default edge to attach pressure modes to
842  // which is part of the inner solve;
843  int defedge = -1;
844 
845  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(nlevels);
846  for(i = 0; i < bndgraphs.size(); ++i)
847  {
848  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
849 
850  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
851  {
852  // find edge in graph vert list
853  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
854  {
855  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
856 
857  if(defedge == -1)
858  {
859  defedge = edgeId;
860  break;
861  }
862  }
863  }
864  if(defedge != -1)
865  {
866  break;
867  }
868  }
869 
870  for(int n = 1; n < nlevels; ++n)
871  {
872  // produce a map with a key that is the element id
873  // that contains which next level patch it belongs to
874  vector<MultiRegions::SubGraphSharedPtr> bndgraphs = bottomUpGraph->GetInteriorBlocks(n+1);
875 
876  // Fill next level graph of adjacent elements and their level
877  map<int,int> ElmtInBndry;
878 
879  for(i = 0; i < bndgraphs.size(); ++i)
880  {
881  int GlobIdOffset = bndgraphs[i]->GetIdOffset();
882 
883  for(j = 0; j < bndgraphs[i]->GetNverts(); ++j)
884  {
885  // find edge in graph vert list
886  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
887  {
888  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
889 
890  if(EdgeIdToElmts[edgeId][0] != -1)
891  {
892  ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
893  }
894  if(EdgeIdToElmts[edgeId][1] != -1)
895  {
896  ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
897  }
898  }
899  }
900  }
901 
902  // Now search interior patches in this level for edges
903  // that share the same element as a boundary edge and
904  // assign this elmt that boundary edge
905  vector<MultiRegions::SubGraphSharedPtr> intgraphs = bottomUpGraph->GetInteriorBlocks(n);
906  for(i = 0; i < intgraphs.size(); ++i)
907  {
908  int GlobIdOffset = intgraphs[i]->GetIdOffset();
909  bool SetEdge = false;
910  int elmtid = 0;
911  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
912  {
913  // Check to see if graph vert is an edge
914  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
915  {
916  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
917 
918  for(k = 0; k < 2; ++k)
919  {
920  // relevant edge id
921  elmtid = EdgeIdToElmts[edgeId][k];
922 
923  if(elmtid != -1)
924  {
925  mapIt = ElmtInBndry.find(elmtid);
926 
927  if(mapIt != ElmtInBndry.end())
928  {
929  // now find a edge in the next level boundary graph
930  int GlobIdOffset1 = bndgraphs[mapIt->second]->GetIdOffset();
931  for(int l = 0; l < bndgraphs[mapIt->second]->GetNverts(); ++l)
932  {
933  // find edge in graph vert list
934  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset1+l) != 0)
935  {
936  //June 2012: commenting this condition apparently
937  //solved the bug caused by the edge reordering procedure
938 
939  //if(AddMeanPressureToEdgeId[elmtid] == -1)
940  //{
941 
942  //AddMeanPressureToEdgeId[elmtid] = HomGraphEdgeIdToEdgeId[GlobIdOffset1+l];
943  AddMeanPressureToEdgeId[elmtid] = defedge;
944 
945  //}
946  SetEdge = true;
947  break;
948  }
949  }
950  }
951  }
952  }
953  }
954  }
955 
956 
957  // if we have failed to find matching edge in next
958  // level patch boundary then set last found elmt
959  // associated to this interior patch to the
960  // default edget value
961  if(SetEdge == false)
962  {
963  if(elmtid == -1) // find an elmtid in patch
964  {
965  for(j = 0; j < intgraphs[i]->GetNverts(); ++j)
966  {
967  if(HomGraphEdgeIdToEdgeId.count(GlobIdOffset+j) != 0)
968  {
969  edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset+j];
970  for(k = 0; k < 2; ++k)
971  {
972  // relevant edge id
973  elmtid = EdgeIdToElmts[edgeId][k];
974  if(elmtid != -1)
975  {
976  break;
977  }
978  }
979  }
980  if(elmtid != -1)
981  {
982  break;
983  }
984  }
985  }
986  if(AddMeanPressureToEdgeId[elmtid] == -1)
987  {
988  AddMeanPressureToEdgeId[elmtid] = defedge;
989  }
990  }
991  }
992  }
993 
994 }
995 
996 
997 
998 
999 }
1000 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:320
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:314
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:765
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:331
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:395
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
static const NekDouble kNekZeroTol
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:826
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:388
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
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:390
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
Definition: AssemblyMap.h:348
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:316
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:352
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:312
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:384
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:350
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:271
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:303
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)
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
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:386