Nektar++
CoupledLocalToGlobalC0ContMap.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: CoupledLocalToGlobalC0ContMap.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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Wrapper class around the library
32// LocalToGlobalC0ContMap class for use in the Couplied Linearised NS
33// solver.
34//
35///////////////////////////////////////////////////////////////////////////////
36
40#include <LocalRegions/SegExp.h>
43
44namespace Nektar
45{
46/**
47 * This is an vector extension of
48 * MultiRegions::AssemblyMapCG::SetUp2DExpansionC0ContMap related to the
49 * Linearised Navier Stokes problem
50 */
53 [[maybe_unused]] const SpatialDomains::MeshGraphSharedPtr &graph,
55 &boundaryConditions,
57 const MultiRegions::ExpListSharedPtr &pressure, const size_t nz_loc,
58 const bool CheckforSingularSys)
59 : AssemblyMapCG(pSession, fields[0]->GetComm())
60{
61 size_t i, j, k, n;
62 size_t cnt = 0, offset = 0;
63 size_t meshVertId;
64 size_t meshEdgeId;
65 size_t bndEdgeCnt;
66 size_t globalId;
67 size_t nEdgeCoeffs;
68 size_t nEdgeInteriorCoeffs;
69 int firstNonDirGraphVertId;
70 size_t nLocBndCondDofs = 0;
71 size_t nLocDirBndCondDofs = 0;
72 int nExtraDirichlet = 0;
76 StdRegions::Orientation edgeOrient;
77 Array<OneD, unsigned int> edgeInteriorMap;
78 Array<OneD, int> edgeInteriorSign;
79 size_t nvel = fields.size();
80
81 const LocalRegions::ExpansionVector &locExpVector = *(fields[0]->GetExp());
82 int id, diff;
83 size_t nel = fields[0]->GetNumElmts();
84
85 MultiRegions::PeriodicMap periodicVerts;
86 MultiRegions::PeriodicMap periodicEdges;
87 MultiRegions::PeriodicMap periodicFaces;
88 std::vector<std::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
112 fields[0]->GetBndCondExpansions();
113
115 bndConditionsVec(nvel);
116 for (i = 0; i < nvel; ++i)
117 {
118 bndConditionsVec[i] = fields[i]->GetBndConditions();
119 }
120
121 std::map<int, int> IsDirVertDof;
122 std::map<int, int> IsDirEdgeDof;
123
125 for (j = 0; j < bndCondExp.size(); ++j)
126 {
127 std::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]
132 ->GetExp(k)
134 ->GetGeom1D();
135 BndExpVids[g->GetVid(0)] = g->GetVid(0);
136 BndExpVids[g->GetVid(1)] = g->GetVid(1);
137 }
138
139 for (i = 0; i < nvel; ++i)
140 {
141 if (bndConditionsVec[i][j]->GetBoundaryConditionType() ==
143 {
144 // set number of Dirichlet conditions along edge
145 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
146 {
147 IsDirEdgeDof[bndCondExp[j]
148 ->GetExp(k)
150 ->GetGeom1D()
151 ->GetGlobalID()] += 1;
152 }
153
154 // Set number of Dirichlet conditions at vertices
155 // with a clamp on its maximum value being nvel to
156 // handle corners between expansions
157 for (auto &mapIt : BndExpVids)
158 {
159 id = IsDirVertDof[mapIt.second] + 1;
160 IsDirVertDof[mapIt.second] = (id > (int)nvel) ? nvel : id;
161 }
162 }
163 else
164 {
165 // Check to see that edge normals have non-zero
166 // component in this direction since otherwise
167 // also can be singular.
168 /// @TODO: Fix this so that we can extract normals from edges
169 for (k = 0; k < bndCondExp[j]->GetNumElmts(); ++k)
170 {
173 bndCondExp[j]
174 ->GetExp(k)
176 locnorm =
177 loc_exp->GetLeftAdjacentElementExp()->GetTraceNormal(
178 loc_exp->GetLeftAdjacentElementTrace());
179
180 size_t ndir = locnorm.size();
181 if (i < ndir) // account for Fourier version where n can be
182 // larger then ndir
183 {
184 for (size_t l = 0; l < locnorm[0].size(); ++l)
185 {
186 if (fabs(locnorm[i][l]) > NekConstants::kNekZeroTol)
187 {
188 m_systemSingular = false;
189 break;
190 }
191 }
192 }
193 if (m_systemSingular == false)
194 {
195 break;
196 }
197 }
198 }
199 }
200 }
201
203
204 Array<OneD, int> AddMeanPressureToEdgeId(nel, -1);
205 size_t edgeId, vertId;
206
207 // special case of singular problem - need to fix one pressure
208 // dof to a dirichlet edge. Since we attached pressure dof to
209 // last velocity component of edge need to make sure this
210 // component is Dirichlet
212 {
213 id = -1;
214 for (i = 0; i < bndConditionsVec[0].size(); ++i)
215 {
216 if (bndConditionsVec[nvel - 1][i]->GetBoundaryConditionType() ==
218 {
219 id = bndCondExp[i]
220 ->GetExp(0)
222 ->GetGeom1D()
223 ->GetGlobalID();
224 break;
225 }
226 }
227
228 ASSERTL0(id != -1, " Did not find an edge to attach singular pressure "
229 "degree of freedom");
230
231 // determine element with this edge id. There may be a
232 // more direct way of getting element from spatialDomains
233 for (i = 0; i < nel; ++i)
234 {
235 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
236 {
237 edgeId = (locExpVector[i]
239 ->GetGeom2D())
240 ->GetEid(j);
241
242 if ((int)edgeId == id)
243 {
244 AddMeanPressureToEdgeId[i] = id;
245 break;
246 }
247 }
248
249 if (AddMeanPressureToEdgeId[i] != -1)
250 {
251 break;
252 }
253 }
254 }
255
256 for (i = 0; i < nel; ++i)
257 {
258 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
259 {
260 vertId =
261 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
262 ->GetVid(j);
263 if (Dofs[0].count(vertId) == 0)
264 {
265 Dofs[0][vertId] = nvel * nz_loc;
266
267 // Adjust for a Dirichlet boundary condition to give number to
268 // be solved
269 if (IsDirVertDof.count(vertId) != 0)
270 {
271 Dofs[0][vertId] -= IsDirVertDof[vertId] * nz_loc;
272 }
273 }
274
275 edgeId =
276 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
277 ->GetEid(j);
278 if (Dofs[1].count(edgeId) == 0)
279 {
280 Dofs[1][edgeId] =
281 nvel * (locExpVector[i]->GetTraceNcoeffs(j) - 2) * nz_loc;
282 }
283
284 // Adjust for Dirichlet boundary conditions to give number to be
285 // solved
286 if (IsDirEdgeDof.count(edgeId) != 0)
287 {
288 Dofs[1][edgeId] -= IsDirEdgeDof[edgeId] * nz_loc *
289 (locExpVector[i]->GetTraceNcoeffs(j) - 2);
290 }
291 }
292 }
293
294 std::set<int> extraDirVerts, extraDirEdges;
295
296 CreateGraph(*fields[0], bndCondExp, bndConditionsVec, false, periodicVerts,
297 periodicEdges, periodicFaces, ReorderedGraphVertId,
298 bottomUpGraph, extraDirVerts, extraDirEdges,
299 firstNonDirGraphVertId, nExtraDirichlet, 4);
300 /*
301 SetUp2DGraphC0ContMap(*fields[0],
302 bndCondExp,
303 bndConditionsVec,
304 periodicVerts, periodicEdges,
305 Dofs, ReorderedGraphVertId,
306 firstNonDirGraphVertId, nExtraDirichlet,
307 bottomUpGraph, extraDir, false, 4);
308 */
309
310 /**
311 * STEP 2a: Set the mean pressure modes to edges depending on
312 * type of direct solver technique;
313 */
314
315 // determine which edge to add mean pressure dof based on
316 // ensuring that at least one pressure dof from an internal
317 // patch is associated with its boundary system
318 if (m_session->MatchSolverInfoAsEnum(
320 {
321
322 FindEdgeIdToAddMeanPressure(ReorderedGraphVertId, nel, locExpVector,
323 edgeId, vertId, firstNonDirGraphVertId,
324 IsDirEdgeDof, bottomUpGraph,
325 AddMeanPressureToEdgeId);
326 }
327
328 // Set unset elmts to non-Dirichlet edges.
329 // special case of singular problem - need to fix one
330 // pressure dof to a dirichlet edge
331 for (i = 0; i < nel; ++i)
332 {
333 for (j = 0; j < (size_t)locExpVector[i]->GetNverts(); ++j)
334 {
335 edgeId =
336 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
337 ->GetEid(j);
338
339 if (IsDirEdgeDof.count(edgeId) == 0) // interior edge
340 {
341 // setup AddMeanPressureToEdgeId to decide where to
342 // put pressure
343 if (AddMeanPressureToEdgeId[i] == -1)
344 {
345 AddMeanPressureToEdgeId[i] = edgeId;
346 }
347 }
348 }
349 ASSERTL0((AddMeanPressureToEdgeId[i] != -1),
350 "Did not determine "
351 "an edge to attach mean pressure dof");
352 // Add the mean pressure degree of freedom to this edge
353 Dofs[1][AddMeanPressureToEdgeId[i]] += nz_loc;
354 }
355
356 std::map<int, int> pressureEdgeOffset;
357
358 /**
359 * STEP 2: Count out the number of Dirichlet vertices and edges first
360 */
361 for (i = 0; i < bndCondExp.size(); i++)
362 {
363 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
364 {
365 bndSegExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
366 for (k = 0; k < nvel; ++k)
367 {
368 if (bndConditionsVec[k][i]->GetBoundaryConditionType() ==
370 {
371 nLocDirBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
372 }
373
374 if (bndConditionsVec[k][i]->GetBoundaryConditionType() !=
376 {
377 nLocBndCondDofs += bndSegExp->GetNcoeffs() * nz_loc;
378 }
379 }
380 }
381 }
382
384 {
385 m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet + nz_loc;
386 }
387 else
388 {
389 m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet;
390 }
391
392 /**
393 * STEP 3: Set up an array which contains the offset information of
394 * the different graph vertices.
395 *
396 * This basically means to identify how many global degrees of
397 * freedom the individual graph vertices correspond. Obviously,
398 * the graph vertices corresponding to the mesh-vertices account
399 * for a single global DOF. However, the graph vertices
400 * corresponding to the element edges correspond to 2*(N-2) global DOF
401 * where N is equal to the number of boundary modes on this edge.
402 */
403 Array<OneD, int> graphVertOffset(
404 nvel * nz_loc *
405 (ReorderedGraphVertId[0].size() + ReorderedGraphVertId[1].size()),
406 0);
407 graphVertOffset[0] = 0;
408
409 m_signChange = false;
410
411 for (i = 0; i < nel; ++i)
412 {
413 locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
414
415 for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
416 {
417 nEdgeCoeffs = locExpansion->GetTraceNcoeffs(j);
418 meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
419 meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
420
421 for (k = 0; k < nvel * nz_loc; ++k)
422 {
423 graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
424 nz_loc +
425 k] = 1;
426 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] * nvel *
427 nz_loc +
428 k] = (nEdgeCoeffs - 2);
429 }
430
431 bType = locExpansion->GetBasisType(0);
432 // need a sign vector for modal expansions if nEdgeCoeffs >=4
433 if ((nEdgeCoeffs >= 4) && (bType == LibUtilities::eModified_A))
434 {
435 m_signChange = true;
436 }
437 }
438 }
439
440 // Add mean pressure modes;
441 for (i = 0; i < nel; ++i)
442 {
443 graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] +
444 1) *
445 nvel * nz_loc -
446 1] += nz_loc;
447 // graphVertOffset[(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]])*nvel*nz_loc]
448 // += nz_loc;
449 }
450
451 // Negate the vertices and edges with only a partial
452 // Dirichlet conditon. Essentially we check to see if an edge
453 // has a mixed Dirichlet with Neumann/Robin Condition and if
454 // so negate the offset associated with this vertex.
455
456 std::map<int, int> DirVertChk;
457
458 for (i = 0; i < bndConditionsVec[0].size(); ++i)
459 {
460 cnt = 0;
461 for (j = 0; j < nvel; ++j)
462 {
463 if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
465 {
466 cnt++;
467 }
468 }
469
470 // Case where partial Dirichlet boundary condition
471 if ((cnt > 0) && (cnt < nvel))
472 {
473 for (j = 0; j < nvel; ++j)
474 {
475 if (bndConditionsVec[j][i]->GetBoundaryConditionType() ==
477 {
478 // negate graph offsets which should be
479 // Dirichlet conditions
480 for (k = 0; k < bndCondExp[i]->GetNumElmts(); ++k)
481 {
482 // vertices with mix condition;
483 id = bndCondExp[i]
484 ->GetExp(k)
486 ->GetGeom1D()
487 ->GetVid(0);
488 if (DirVertChk.count(id * nvel + j) == 0)
489 {
490 DirVertChk[id * nvel + j] = 1;
491 for (n = 0; n < nz_loc; ++n)
492 {
493 graphVertOffset[ReorderedGraphVertId[0][id] *
494 nvel * nz_loc +
495 j * nz_loc + n] *= -1;
496 }
497 }
498
499 id = bndCondExp[i]
500 ->GetExp(k)
502 ->GetGeom1D()
503 ->GetVid(1);
504 if (DirVertChk.count(id * nvel + j) == 0)
505 {
506 DirVertChk[id * nvel + j] = 1;
507 for (n = 0; n < nz_loc; ++n)
508 {
509 graphVertOffset[ReorderedGraphVertId[0][id] *
510 nvel * nz_loc +
511 j * nz_loc + n] *= -1;
512 }
513 }
514
515 // edges with mixed id;
516 id = bndCondExp[i]
517 ->GetExp(k)
519 ->GetGeom1D()
520 ->GetGlobalID();
521 for (n = 0; n < nz_loc; ++n)
522 {
523 graphVertOffset[ReorderedGraphVertId[1][id] * nvel *
524 nz_loc +
525 j * nz_loc + n] *= -1;
526 }
527 }
528 }
529 }
530 }
531 }
532
533 cnt = 0;
534 // assemble accumulative list of full Dirichlet values.
535 for (i = 0; i < firstNonDirGraphVertId * nvel * nz_loc; ++i)
536 {
537 diff = abs(graphVertOffset[i]);
538 graphVertOffset[i] = cnt;
539 cnt += diff;
540 }
541
542 // set Dirichlet values with negative values to Dirichlet value
543 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
544 ++i)
545 {
546 if (graphVertOffset[i] < 0)
547 {
548 diff = -graphVertOffset[i];
549 graphVertOffset[i] = -cnt;
550 cnt += diff;
551 }
552 }
553
554 // Accumulate all interior degrees of freedom with positive values
556
557 // offset values
558 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
559 ++i)
560 {
561 if (graphVertOffset[i] >= 0)
562 {
563 diff = graphVertOffset[i];
564 graphVertOffset[i] = cnt;
565 cnt += diff;
566 }
567 }
568
569 // Finally set negative entries (corresponding to Dirichlet
570 // values ) to be positive
571 for (i = firstNonDirGraphVertId * nvel * nz_loc; i < graphVertOffset.size();
572 ++i)
573 {
574 if (graphVertOffset[i] < 0)
575 {
576 graphVertOffset[i] = -graphVertOffset[i];
577 }
578 }
579
580 // Allocate the proper amount of space for the class-data and fill
581 // information that is already known
582 cnt = 0;
585
586 for (i = 0; i < nel; ++i)
587 {
589 nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
590 // add these coeffs up separately since
591 // pressure->GetNcoeffs can include the coefficient in
592 // multiple planes.
593 m_numLocalCoeffs += (pressure->GetExp(i)->GetNcoeffs() - 1) * nz_loc;
594 }
595
597
600 m_bndCondIDToGlobalTraceID = Array<OneD, int>(nLocBndCondDofs, -1);
601
602 // Set default sign array.
605
606 m_staticCondLevel = staticCondLevel;
607 m_numPatches = nel;
608
611
612 for (i = 0; i < nel; ++i)
613 {
615 (unsigned int)nz_loc *
616 (nvel * locExpVector[i]->NumBndryCoeffs() + 1);
618 (unsigned int)nz_loc * (pressure->GetExp(i)->GetNcoeffs() - 1);
619 }
620
621 // Set up local to local bnd and local int maps
625
626 size_t bndcnt = 0;
627 size_t intcnt = 0;
628 cnt = 0;
629 for (i = 0; i < nel; ++i)
630 {
631 for (j = 0; j < nz_loc * (nvel * locExpVector[i]->NumBndryCoeffs());
632 ++j)
633 {
634 m_localToLocalBndMap[bndcnt++] = cnt++;
635 }
636
637 for (n = 0; n < nz_loc; ++n)
638 {
639 m_localToLocalBndMap[bndcnt++] = cnt++;
640 for (j = 1; j < (size_t)pressure->GetExp(i)->GetNcoeffs(); ++j)
641 {
642 m_localToLocalIntMap[intcnt++] = cnt++;
643 }
644 }
645 }
646
647 /**
648 * STEP 4: Now, all ingredients are ready to set up the actual
649 * local to global mapping.
650 *
651 * The remainder of the map consists of the element-interior
652 * degrees of freedom. This leads to the block-diagonal submatrix
653 * as each element-interior mode is globally orthogonal to modes
654 * in all other elements.
655 */
656 cnt = 0;
657 size_t nv, velnbndry;
659
660 // Loop over all the elements in the domain in shuffled
661 // ordering (element type consistency)
662 for (i = 0; i < nel; ++i)
663 {
664 locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
665
666 velnbndry = locExpansion->NumBndryCoeffs();
667
668 // Require an inverse ordering of the bmap system to store
669 // local numbering system. Therefore get hold of elemental
670 // bmap and set up an inverse map
671 std::map<int, int> inv_bmap;
672 locExpansion->GetBoundaryMap(bmap);
673 for (j = 0; j < bmap.size(); ++j)
674 {
675 inv_bmap[bmap[j]] = j;
676 }
677
678 // Loop over all edges (and vertices) of element i
679 for (j = 0; j < (size_t)locExpansion->GetNtraces(); ++j)
680 {
681 nEdgeInteriorCoeffs = locExpansion->GetTraceNcoeffs(j) - 2;
682 edgeOrient = (locExpansion->GetGeom())->GetEorient(j);
683 meshEdgeId = (locExpansion->GetGeom())->GetEid(j);
684 meshVertId = (locExpansion->GetGeom())->GetVid(j);
685
686 auto pIt = periodicEdges.find(meshEdgeId);
687
688 // See if this edge is periodic. If it is, then we map all
689 // edges to the one with lowest ID, and align all
690 // coefficients to this edge orientation.
691 if (pIt != periodicEdges.end())
692 {
693 std::pair<int, StdRegions::Orientation> idOrient =
694 DeterminePeriodicEdgeOrientId(meshEdgeId, edgeOrient,
695 pIt->second);
696 edgeOrient = idOrient.second;
697 }
698
699 locExpansion->GetTraceInteriorToElementMap(
700 j, edgeInteriorMap, edgeInteriorSign, edgeOrient);
701
702 // Set the global DOF for vertex j of element i
703 for (nv = 0; nv < nvel * nz_loc; ++nv)
704 {
705 m_localToGlobalMap[cnt + nv * velnbndry +
706 inv_bmap[locExpansion->GetVertexMap(j)]] =
707 graphVertOffset[ReorderedGraphVertId[0][meshVertId] * nvel *
708 nz_loc +
709 nv];
710
711 // Set the global DOF's for the interior modes of edge j
712 for (k = 0; k < nEdgeInteriorCoeffs; ++k)
713 {
714 m_localToGlobalMap[cnt + nv * velnbndry +
715 inv_bmap[edgeInteriorMap[k]]] =
716 graphVertOffset[ReorderedGraphVertId[1][meshEdgeId] *
717 nvel * nz_loc +
718 nv] +
719 k;
720 }
721 }
722
723 // Fill the sign vector if required
724 if (m_signChange)
725 {
726 for (nv = 0; nv < nvel * nz_loc; ++nv)
727 {
728 for (k = 0; k < nEdgeInteriorCoeffs; ++k)
729 {
730 m_localToGlobalSign[cnt + nv * velnbndry +
731 inv_bmap[edgeInteriorMap[k]]] =
732 (NekDouble)edgeInteriorSign[k];
733 }
734 }
735 }
736 }
737
738 // use difference between two edges of the
739 // AddMeanPressureEdgeId to det nEdgeInteriorCoeffs.
740 nEdgeInteriorCoeffs =
741 graphVertOffset[(ReorderedGraphVertId[1]
742 [AddMeanPressureToEdgeId[i]]) *
743 nvel * nz_loc +
744 1] -
745 graphVertOffset[(ReorderedGraphVertId[1]
746 [AddMeanPressureToEdgeId[i]]) *
747 nvel * nz_loc];
748
749 size_t psize = pressure->GetExp(i)->GetNcoeffs();
750 for (n = 0; n < nz_loc; ++n)
751 {
752 m_localToGlobalMap[cnt + nz_loc * nvel * velnbndry + n * psize] =
753 graphVertOffset
754 [(ReorderedGraphVertId[1][AddMeanPressureToEdgeId[i]] + 1) *
755 nvel * nz_loc -
756 1] +
757 nEdgeInteriorCoeffs +
758 pressureEdgeOffset[AddMeanPressureToEdgeId[i]];
759
760 pressureEdgeOffset[AddMeanPressureToEdgeId[i]] += 1;
761 }
762
763 cnt += (velnbndry * nvel + psize) * nz_loc;
764 }
765
766 // Set up the mapping for the boundary conditions
767 offset = cnt = 0;
768 for (nv = 0; nv < nvel; ++nv)
769 {
770 for (i = 0; i < bndCondExp.size(); i++)
771 {
772 if (bndConditionsVec[nv][i]->GetBoundaryConditionType() ==
774 {
775 continue;
776 }
777
778 for (n = 0; n < nz_loc; ++n)
779 {
780 int ncoeffcnt = 0;
781 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
782 {
783 bndSegExp =
784 bndCondExp[i]->GetExp(j)->as<LocalRegions::SegExp>();
785
786 cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
787 for (k = 0; k < 2; k++)
788 {
789 meshVertId = (bndSegExp->GetGeom1D())->GetVid(k);
791 bndSegExp->GetVertexMap(k)] =
792 graphVertOffset[ReorderedGraphVertId[0]
793 [meshVertId] *
794 nvel * nz_loc +
795 nv * nz_loc + n];
796 }
797
798 meshEdgeId = (bndSegExp->GetGeom1D())->GetGlobalID();
799 bndEdgeCnt = 0;
800 nEdgeCoeffs = bndSegExp->GetNcoeffs();
801 for (k = 0; k < nEdgeCoeffs; k++)
802 {
803 if (m_bndCondIDToGlobalTraceID[cnt + k] == -1)
804 {
806 graphVertOffset
807 [ReorderedGraphVertId[1][meshEdgeId] *
808 nvel * nz_loc +
809 nv * nz_loc + n] +
810 bndEdgeCnt;
811 bndEdgeCnt++;
812 }
813 }
814 ncoeffcnt += nEdgeCoeffs;
815 }
816 // Note: Can not use bndCondExp[i]->GetNcoeffs()
817 // due to homogeneous extension not returning just
818 // the value per plane
819 offset += ncoeffcnt;
820 }
821 }
822 }
823
824 globalId = Vmath::Vmax(m_numLocalCoeffs, &m_localToGlobalMap[0], 1) + 1;
825 m_numGlobalBndCoeffs = globalId;
826
827 /**
828 * STEP 5: The boundary condition mapping is generated from the
829 * same vertex renumbering and fill in a unique interior map.
830 */
831 cnt = 0;
832 for (i = 0; i < (size_t)m_numLocalCoeffs; ++i)
833 {
834 if (m_localToGlobalMap[i] == -1)
835 {
836 m_localToGlobalMap[i] = globalId++;
837 }
838 else
839 {
840 if (m_signChange)
841 {
843 }
845 }
846 }
847 m_numGlobalCoeffs = globalId;
848
849 // Set up the local to global map for the next level when using
850 // multi-level static condensation
851 if (m_session->MatchSolverInfoAsEnum(
853 {
854 if (m_staticCondLevel < (bottomUpGraph->GetNlevels() - 1))
855 {
856 Array<OneD, int> vwgts_perm(Dofs[0].size() + Dofs[1].size() -
857 firstNonDirGraphVertId);
858 for (i = 0; i < locExpVector.size(); ++i)
859 {
860 locExpansion = locExpVector[i]->as<LocalRegions::Expansion2D>();
861 size_t nVert = locExpansion->GetNverts();
862 for (j = 0; j < nVert; ++j)
863 {
864 meshEdgeId = (locExpansion->GetGeom2D())->GetEid(j);
865 meshVertId = (locExpansion->GetGeom2D())->GetVid(j);
866
867 if (ReorderedGraphVertId[0][meshVertId] >=
868 firstNonDirGraphVertId)
869 {
870 vwgts_perm[ReorderedGraphVertId[0][meshVertId] -
871 firstNonDirGraphVertId] =
872 Dofs[0][meshVertId];
873 }
874
875 if (ReorderedGraphVertId[1][meshEdgeId] >=
876 firstNonDirGraphVertId)
877 {
878 vwgts_perm[ReorderedGraphVertId[1][meshEdgeId] -
879 firstNonDirGraphVertId] =
880 Dofs[1][meshEdgeId];
881 }
882 }
883 }
884
885 bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
886
889 bottomUpGraph);
890 }
891 }
892}
893
895 std::vector<std::map<int, int>> &ReorderedGraphVertId, size_t &nel,
896 const LocalRegions::ExpansionVector &locExpVector, size_t &edgeId,
897 [[maybe_unused]] size_t &vertId, int &firstNonDirGraphVertId,
898 std::map<int, int> &IsDirEdgeDof,
900 Array<OneD, int> &AddMeanPressureToEdgeId)
901{
902 size_t i, j, k;
903
904 // Make list of homogeneous graph edges to elmt mappings
905 Array<TwoD, int> EdgeIdToElmts(ReorderedGraphVertId[1].size(), 2, -1);
906 std::map<int, int> HomGraphEdgeIdToEdgeId;
907
908 for (i = 0; i < nel; ++i)
909 {
910 size_t nVert = locExpVector[i]->GetNverts();
911 for (j = 0; j < nVert; ++j)
912 {
913 edgeId =
914 (locExpVector[i]->as<LocalRegions::Expansion2D>()->GetGeom2D())
915 ->GetEid(j);
916
917 // note second condition stops us using mixed boundary condition
918 if ((ReorderedGraphVertId[1][edgeId] >= firstNonDirGraphVertId) &&
919 (IsDirEdgeDof.count(edgeId) == 0))
920 {
921 HomGraphEdgeIdToEdgeId[ReorderedGraphVertId[1][edgeId] -
922 firstNonDirGraphVertId] = edgeId;
923
924 if (EdgeIdToElmts[edgeId][0] == -1)
925 {
926 EdgeIdToElmts[edgeId][0] = i;
927 }
928 else
929 {
930 EdgeIdToElmts[edgeId][1] = i;
931 }
932 }
933 }
934 }
935
936 // Start at second to last level and find edge on boundary
937 // to attach element
938 size_t nlevels = bottomUpGraph->GetNlevels();
939
940 // determine a default edge to attach pressure modes to
941 // which is part of the inner solve;
942 int defedge = -1;
943
944 std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
945 bottomUpGraph->GetInteriorBlocks(nlevels);
946 for (i = 0; i < bndgraphs.size(); ++i)
947 {
948 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
949 size_t nVert = bndgraphs[i]->GetNverts();
950
951 for (j = 0; j < nVert; ++j)
952 {
953 // find edge in graph vert list
954 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
955 {
956 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
957
958 if (defedge == -1)
959 {
960 defedge = edgeId;
961 break;
962 }
963 }
964 }
965 if (defedge != -1)
966 {
967 break;
968 }
969 }
970
971 for (size_t n = 1; n < nlevels; ++n)
972 {
973 // produce a map with a key that is the element id
974 // that contains which next level patch it belongs to
975 std::vector<MultiRegions::SubGraphSharedPtr> bndgraphs =
976 bottomUpGraph->GetInteriorBlocks(n + 1);
977
978 // Fill next level graph of adjacent elements and their level
979 std::map<int, int> ElmtInBndry;
980
981 for (i = 0; i < bndgraphs.size(); ++i)
982 {
983 int GlobIdOffset = bndgraphs[i]->GetIdOffset();
984 size_t nVert = bndgraphs[i]->GetNverts();
985
986 for (j = 0; j < nVert; ++j)
987 {
988 // find edge in graph vert list
989 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
990 {
991 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
992
993 if (EdgeIdToElmts[edgeId][0] != -1)
994 {
995 ElmtInBndry[EdgeIdToElmts[edgeId][0]] = i;
996 }
997 if (EdgeIdToElmts[edgeId][1] != -1)
998 {
999 ElmtInBndry[EdgeIdToElmts[edgeId][1]] = i;
1000 }
1001 }
1002 }
1003 }
1004
1005 // Now search interior patches in this level for edges
1006 // that share the same element as a boundary edge and
1007 // assign this elmt that boundary edge
1008 std::vector<MultiRegions::SubGraphSharedPtr> intgraphs =
1009 bottomUpGraph->GetInteriorBlocks(n);
1010 for (i = 0; i < intgraphs.size(); ++i)
1011 {
1012 int GlobIdOffset = intgraphs[i]->GetIdOffset();
1013 bool SetEdge = false;
1014 int elmtid = 0;
1015 size_t nVert = intgraphs[i]->GetNverts();
1016 for (j = 0; j < nVert; ++j)
1017 {
1018 // Check to see if graph vert is an edge
1019 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
1020 {
1021 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
1022
1023 for (k = 0; k < 2; ++k)
1024 {
1025 // relevant edge id
1026 elmtid = EdgeIdToElmts[edgeId][k];
1027
1028 if (elmtid != -1)
1029 {
1030 auto mapIt = ElmtInBndry.find(elmtid);
1031
1032 if (mapIt != ElmtInBndry.end())
1033 {
1034 // now find a edge in the next level boundary
1035 // graph
1036 int GlobIdOffset1 =
1037 bndgraphs[mapIt->second]->GetIdOffset();
1038 for (int l = 0;
1039 l < bndgraphs[mapIt->second]->GetNverts();
1040 ++l)
1041 {
1042 // find edge in graph vert list
1043 if (HomGraphEdgeIdToEdgeId.count(
1044 GlobIdOffset1 + l) != 0)
1045 {
1046 AddMeanPressureToEdgeId[elmtid] =
1047 defedge;
1048
1049 SetEdge = true;
1050 break;
1051 }
1052 }
1053 }
1054 }
1055 }
1056 }
1057 }
1058
1059 // if we have failed to find matching edge in next
1060 // level patch boundary then set last found elmt
1061 // associated to this interior patch to the
1062 // default edget value
1063 if (SetEdge == false)
1064 {
1065 if (elmtid == -1) // find an elmtid in patch
1066 {
1067 size_t nVert = intgraphs[i]->GetNverts();
1068 for (j = 0; j < nVert; ++j)
1069 {
1070 if (HomGraphEdgeIdToEdgeId.count(GlobIdOffset + j) != 0)
1071 {
1072 edgeId = HomGraphEdgeIdToEdgeId[GlobIdOffset + j];
1073 for (k = 0; k < 2; ++k)
1074 {
1075 // relevant edge id
1076 elmtid = EdgeIdToElmts[edgeId][k];
1077 if (elmtid != -1)
1078 {
1079 break;
1080 }
1081 }
1082 }
1083 if (elmtid != -1)
1084 {
1085 break;
1086 }
1087 }
1088 }
1089 if (AddMeanPressureToEdgeId[elmtid] == -1)
1090 {
1091 AddMeanPressureToEdgeId[elmtid] = defedge;
1092 }
1093 }
1094 }
1095 }
1096}
1097
1098} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
CoupledLocalToGlobalC0ContMap(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &graph, const SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const MultiRegions::ExpListSharedPtr &pressure, const size_t nz_loc, const bool CheeckForSingularSys=true)
void FindEdgeIdToAddMeanPressure(std::vector< std::map< int, int > > &ReorderedGraphVertId, size_t &nel, const LocalRegions::ExpansionVector &locExpVector, size_t &edgeId, size_t &vertId, int &firstNonDirGraphVertId, std::map< int, int > &IsDirEdgeDof, MultiRegions::BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, Array< OneD, int > &AddMeanPressureToEdgeId)
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Constructs mappings for the C0 scalar continuous Galerkin formulation.
Definition: AssemblyMapCG.h:66
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)
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:361
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:375
Array< OneD, int > m_localToLocalIntMap
Integer map of local boundary coeffs to local interior system numbering.
Definition: AssemblyMap.h:384
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:372
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:380
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:427
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:330
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:342
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:434
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:423
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:346
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:348
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:429
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:350
Array< OneD, int > m_localToGlobalBndMap
Integer map of local coeffs to global Boundary Dofs.
Definition: AssemblyMap.h:378
Array< OneD, int > m_localToLocalBndMap
Integer map of local boundary coeffs to local boundary system numbering.
Definition: AssemblyMap.h:382
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:425
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:344
Array< OneD, int > m_bndCondIDToGlobalTraceID
Integer map of bnd cond trace number to global trace number.
Definition: AssemblyMap.h:392
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:688
int GetNverts() const
This function returns the number of vertices of the expansion domain.
Definition: StdExpansion.h:246
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:261
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:253
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:46
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:50
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:68
std::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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...
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
static const NekDouble kNekZeroTol
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:289
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
double NekDouble
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.hpp:644
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289