Nektar++
DisContField.cpp
Go to the documentation of this file.
1//////////////////////////////////////////////////////////////////////////////
2//
3// File: DisContField.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: Field definition for 1D domain with boundary
33// conditions using LDG-H
34//
35///////////////////////////////////////////////////////////////////////////////
36
40#include <LocalRegions/HexExp.h>
42#include <LocalRegions/TetExp.h>
43#include <LocalRegions/TriExp.h>
47
48using namespace std;
49
51{
52/**
53 * @class DisContField
54 * This class augments the list of local expansions inherited from
55 * ExpList with boundary conditions. Inter-element boundaries are
56 * handled using an discontinuous Galerkin scheme.
57 */
58
59/**
60 * Constructs an empty expansion list with no boundary conditions.
61 */
63 : ExpList(), m_bndConditions(), m_bndCondExpansions(),
65{
66}
67
68/**
69 * An expansion list for the boundary expansions is generated first for
70 * the field. These are subsequently evaluated for time zero. The trace
71 * map is then constructed.
72 * @param graph1D A mesh, containing information about the domain
73 * and the spectral/hp element expansions.
74 * @param bcs Information about the enforced boundary
75 * conditions.
76 * @param variable The session variable associated with the
77 * boundary conditions to enforce.
78 * @param solnType Type of global system to use.
79 */
82 const std::string &variable, const bool SetUpJustDG,
83 const bool DeclareCoeffPhysArrays,
85 const std::string bcvariable)
86 : ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
88{
89 std::string bcvar;
90 if (bcvariable == "NotSet")
91 {
92 bcvar = variable;
93 }
94 else
95 {
96 bcvar = bcvariable;
97 }
98
99 if (bcvar.compare("DefaultVar") != 0)
100 {
102
103 GenerateBoundaryConditionExpansion(graph, bcs, bcvar,
104 DeclareCoeffPhysArrays);
105 if (DeclareCoeffPhysArrays)
106 {
107 EvaluateBoundaryConditions(0.0, bcvar);
108 }
110 FindPeriodicTraces(bcs, bcvar);
111 }
112
113 int i, cnt;
114 Array<OneD, int> ElmtID, TraceID;
115 GetBoundaryToElmtMap(ElmtID, TraceID);
116
117 for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
118 {
120 locExpList = m_bndCondExpansions[i];
121
122 for (int e = 0; e < locExpList->GetExpSize(); ++e)
123 {
124 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
125 locExpList->GetExp(e));
126 locExpList->GetExp(e)->SetAdjacentElementExp(
127 TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
128 }
129 cnt += m_bndCondExpansions[i]->GetExpSize();
130 }
131
132 if (SetUpJustDG)
133 {
134 SetUpDG(variable, ImpType);
135 }
136 else
137 {
138 if (m_session->DefinesSolverInfo("PROJECTION"))
139 {
140 std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
141 if ((ProjectStr == "MixedCGDG") ||
142 (ProjectStr == "Mixed_CG_Discontinuous"))
143 {
144 SetUpDG(variable);
145 }
146 else
147 {
149 }
150 }
151 else
152 {
154 }
155 }
156}
157
158/**
159 * @brief Set up all DG member variables and maps.
160 */
161void DisContField::SetUpDG(const std::string variable,
163{
164 // Check for multiple calls
166 {
167 return;
168 }
169
170 // Set up matrix map
172
173 // Set up trace space
176 m_comm, true, "DefaultVar", ImpType);
177
178 PeriodicMap periodicTraces = (m_expType == e1D) ? m_periodicVerts
181
184 m_bndConditions, periodicTraces, variable);
185
186 if (m_session->DefinesCmdLineArgument("verbose"))
187 {
188 m_traceMap->PrintStats(std::cout, variable);
189 }
190
192 m_traceMap->GetElmtToTrace();
193
194 // Scatter trace to elements. For each element, we find
195 // the trace associated to each elemental trace. The element
196 // then retains a pointer to the elemental trace, to
197 // ensure uniqueness of normals when retrieving from two
198 // adjoining elements which do not lie in a plane.
199 for (int i = 0; i < m_exp->size(); ++i)
200 {
201 for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
202 {
203 LocalRegions::ExpansionSharedPtr exp = elmtToTrace[i][j];
204 ;
205
206 exp->SetAdjacentElementExp(j, (*m_exp)[i]);
207 (*m_exp)[i]->SetTraceExp(j, exp);
208 }
209 }
210
211 // Set up physical normals
213
214 // Create interface exchange object
217
218 int cnt, n;
219
220 // Identify boundary trace
221 for (cnt = 0, n = 0; n < m_bndCondExpansions.size(); ++n)
222 {
223 if (m_bndConditions[n]->GetBoundaryConditionType() !=
225 {
226 for (int v = 0; v < m_bndCondExpansions[n]->GetExpSize(); ++v)
227 {
228 m_boundaryTraces.insert(
229 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + v));
230 }
231 cnt += m_bndCondExpansions[n]->GetExpSize();
232 }
233 }
234
235 // Set up information for periodic boundary conditions.
236 std::unordered_map<int, pair<int, int>> perTraceToExpMap;
237 for (cnt = n = 0; n < m_exp->size(); ++n)
238 {
239 for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
240 {
241 auto it = periodicTraces.find((*m_exp)[n]->GetGeom()->GetTid(v));
242
243 if (it != periodicTraces.end())
244 {
245 perTraceToExpMap[it->first] = make_pair(n, v);
246 }
247 }
248 }
249
250 // Set up left-adjacent tracelist.
251 m_leftAdjacentTraces.resize(cnt);
252
253 // count size of trace
254 for (cnt = n = 0; n < m_exp->size(); ++n)
255 {
256 for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
257 {
259 }
260 }
261
262 // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
263 for (cnt = n = 0; n < m_exp->size(); ++n)
264 {
265 for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
266 {
267 int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
268
269 // Check to see if this trace is periodic.
270 auto it = periodicTraces.find(GeomId);
271
272 if (it != periodicTraces.end())
273 {
274 const PeriodicEntity &ent = it->second[0];
275 auto it2 = perTraceToExpMap.find(ent.id);
276
277 if (it2 == perTraceToExpMap.end())
278 {
279 if (m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
280 !ent.isLocal)
281 {
282 continue;
283 }
284 else
285 {
286 ASSERTL1(false, "Periodic trace not found!");
287 }
288 }
289
291 "Periodic trace in non-forward space?");
292
293 int offset = m_trace->GetPhys_Offset(
294 (m_traceMap->GetElmtToTrace())[n][v]->GetElmtId());
295
296 int offset2 = m_trace->GetPhys_Offset(
297 (m_traceMap->GetElmtToTrace())[it2->second.first]
298 [it2->second.second]
299 ->GetElmtId());
300
301 switch (m_expType)
302 {
303 case e1D:
304 {
305 m_periodicFwdCopy.push_back(offset);
306 m_periodicBwdCopy.push_back(offset2);
307 }
308 break;
309 case e2D:
310 {
311 // Calculate relative orientations between trace to
312 // calculate copying map.
313 int nquad = elmtToTrace[n][v]->GetNumPoints(0);
314
315 vector<int> tmpBwd(nquad);
316 vector<int> tmpFwd(nquad);
317
318 if (ent.orient == StdRegions::eForwards)
319 {
320 for (int i = 0; i < nquad; ++i)
321 {
322 tmpBwd[i] = offset2 + i;
323 tmpFwd[i] = offset + i;
324 }
325 }
326 else
327 {
328 for (int i = 0; i < nquad; ++i)
329 {
330 tmpBwd[i] = offset2 + i;
331 tmpFwd[i] = offset + nquad - i - 1;
332 }
333 }
334
335 for (int i = 0; i < nquad; ++i)
336 {
337 m_periodicFwdCopy.push_back(tmpFwd[i]);
338 m_periodicBwdCopy.push_back(tmpBwd[i]);
339 }
340 }
341 break;
342 case e3D:
343 {
344 // Calculate relative orientations between faces to
345 // calculate copying map.
346 int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
347 int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
348
349 vector<int> tmpBwd(nquad1 * nquad2);
350 vector<int> tmpFwd(nquad1 * nquad2);
351
352 if (ent.orient ==
354 ent.orient ==
356 ent.orient ==
359 {
360 for (int i = 0; i < nquad2; ++i)
361 {
362 for (int j = 0; j < nquad1; ++j)
363 {
364 tmpBwd[i * nquad1 + j] =
365 offset2 + i * nquad1 + j;
366 tmpFwd[i * nquad1 + j] =
367 offset + j * nquad2 + i;
368 }
369 }
370 }
371 else
372 {
373 for (int i = 0; i < nquad2; ++i)
374 {
375 for (int j = 0; j < nquad1; ++j)
376 {
377 tmpBwd[i * nquad1 + j] =
378 offset2 + i * nquad1 + j;
379 tmpFwd[i * nquad1 + j] =
380 offset + i * nquad1 + j;
381 }
382 }
383 }
384
385 if (ent.orient ==
387 ent.orient ==
389 ent.orient ==
392 {
393 // Reverse x direction
394 for (int i = 0; i < nquad2; ++i)
395 {
396 for (int j = 0; j < nquad1 / 2; ++j)
397 {
398 swap(tmpFwd[i * nquad1 + j],
399 tmpFwd[i * nquad1 + nquad1 - j - 1]);
400 }
401 }
402 }
403
404 if (ent.orient ==
406 ent.orient ==
408 ent.orient ==
411 {
412 // Reverse y direction
413 for (int j = 0; j < nquad1; ++j)
414 {
415 for (int i = 0; i < nquad2 / 2; ++i)
416 {
417 swap(tmpFwd[i * nquad1 + j],
418 tmpFwd[(nquad2 - i - 1) * nquad1 + j]);
419 }
420 }
421 }
422
423 for (int i = 0; i < nquad1 * nquad2; ++i)
424 {
425 m_periodicFwdCopy.push_back(tmpFwd[i]);
426 m_periodicBwdCopy.push_back(tmpBwd[i]);
427 }
428 }
429 break;
430 default:
431 ASSERTL1(false, "not set up");
432 }
433 }
434 }
435 }
436
438 *this, m_trace, elmtToTrace, m_leftAdjacentTraces);
439}
440
441/**
442 * @brief This routine determines if an element is to the "left"
443 * of the adjacent trace, which arises from the idea there is a local
444 * normal direction between two elements (i.e. on the trace)
445 * and one elements would then be the left.
446 *
447 * This is typically required since we only define one normal
448 * direction along the trace which goes from the left adjacent
449 * element to the right adjacent element. It is also useful
450 * in DG discretisations to set up a local Riemann problem
451 * where we need to have the concept of a local normal
452 * direction which is unique between two elements.
453 *
454 * There are two cases to be checked:
455 *
456 * 1) First is the trace on a boundary condition (perioidic or
457 * otherwise) or on a partition boundary. If a partition
458 * boundary then trace is always considered to be the left
459 * adjacent trace and the normal is pointing outward of the
460 * soltuion domain. We have to perform an additional case for
461 * a periodic boundary where wer chose the element with the
462 * lowest global id. If the trace is on a parallel partition
463 * we use a member of the traceMap where one side is chosen to
464 * contribute to a unique map and have a value which is not
465 * -1 so this is the identifier for the left adjacent side
466 *
467 * 2) If the element is a elemental boundary on one element
468 * the left adjacent element is defined by a link to the left
469 * element from the trace expansion and this is consistent
470 * with the defitiion of the normal which is determined by the
471 * largest id (in contrast to the periodic boundary definition
472 * !). This reversal of convention does not really matter
473 * providing the normals are defined consistently.
474 */
475
476bool DisContField::IsLeftAdjacentTrace(const int n, const int e)
477{
479 m_traceMap->GetElmtToTrace()[n][e];
480
481 PeriodicMap periodicTraces = (m_expType == e1D) ? m_periodicVerts
484
485 bool fwd = true;
486 if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
487 traceEl->GetRightAdjacentElementTrace() == -1)
488 {
489 // Boundary edge (1 connected element). Do nothing in
490 // serial.
491 auto it = m_boundaryTraces.find(traceEl->GetElmtId());
492
493 // If the edge does not have a boundary condition set on
494 // it, then assume it is a partition edge or periodic.
495 if (it == m_boundaryTraces.end())
496 {
497 fwd = true;
498 }
499 }
500 else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
501 traceEl->GetRightAdjacentElementTrace() != -1)
502 {
503 // Non-boundary edge (2 connected elements).
504 fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*m_exp)[n].get();
505 }
506 else
507 {
508 ASSERTL2(false, "Unconnected trace element!");
509 }
510
511 return fwd;
512}
513
514// Given all boundary regions for the whole solution determine
515// which ones (if any) are part of domain and proivde all other
516// conditions are given as UserDefined Dirichlet.
518 const SpatialDomains::CompositeMap &domain,
520 [[maybe_unused]] const std::string &variable)
521{
523
524 returnval =
526
527 map<int, int> GeometryToRegionsMap;
528
530 Allbcs.GetBoundaryRegions();
532 Allbcs.GetBoundaryConditions();
533
534 // Set up a map of all boundary regions
535 for (auto &it : bregions)
536 {
537 for (auto &bregionIt : *it.second)
538 {
539 // can assume that all regions only contain one point in 1D
540 // Really do not need loop above
541 int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
542 GeometryToRegionsMap[id] = it.first;
543 }
544 }
545
546 map<int, SpatialDomains::GeometrySharedPtr> EndOfDomain;
547
548 // Now find out which points in domain have only one vertex
549 for (auto &domIt : domain)
550 {
551 SpatialDomains::CompositeSharedPtr geomvector = domIt.second;
552 for (int i = 0; i < geomvector->m_geomVec.size(); ++i)
553 {
554 for (int j = 0; j < 2; ++j)
555 {
556 int vid = geomvector->m_geomVec[i]->GetVid(j);
557 if (EndOfDomain.count(vid) == 0)
558 {
559 EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
560 }
561 else
562 {
563 EndOfDomain.erase(vid);
564 }
565 }
566 }
567 }
568 ASSERTL1(EndOfDomain.size() == 2, "Did not find two ends of domain");
569
570 for (auto &regIt : EndOfDomain)
571 {
572 if (GeometryToRegionsMap.count(regIt.first) != 0)
573 {
574 // Set up boundary condition up
575 auto iter = GeometryToRegionsMap.find(regIt.first);
576 ASSERTL1(iter != GeometryToRegionsMap.end(),
577 "Failied to find GeometryToRegionMap");
578
579 int regionId = iter->second;
580 auto bregionsIter = bregions.find(regionId);
581 ASSERTL1(bregionsIter != bregions.end(),
582 "Failed to find boundary region");
583
584 SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
585 returnval->AddBoundaryRegions(regionId, breg);
586
587 auto bconditionsIter = bconditions.find(regionId);
588 ASSERTL1(bconditionsIter != bconditions.end(),
589 "Failed to find boundary collection");
590
592 bconditionsIter->second;
593 returnval->AddBoundaryConditions(regionId, bcond);
594 }
595 }
596
597 return returnval;
598}
599
600/**
601 * Constructor for use in multidomain computations where a
602 * domain list can be passed instead of graph1D
603 *
604 * @param domain Subdomain specified in the inputfile from
605 * which the DisContField is set up
606 */
609 const SpatialDomains::CompositeMap &domain,
611 const std::string &variable,
612 const LibUtilities::CommSharedPtr &comm,
613 bool SetToOneSpaceDimension,
615 : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
616 comm, ImpType)
617{
618 if (variable.compare("DefaultVar") != 0)
619 {
621 GetDomainBCs(domain, Allbcs, variable);
622
624 EvaluateBoundaryConditions(0.0, variable);
626 FindPeriodicTraces(*DomBCs, variable);
627 }
628
629 SetUpDG(variable);
630}
631
632/**
633 * Constructs a field as a copy of an existing field.
634 * @param In Existing DisContField object to copy.
635 */
637 const bool DeclareCoeffPhysArrays)
638 : ExpList(In, DeclareCoeffPhysArrays), m_bndConditions(In.m_bndConditions),
639 m_bndCondExpansions(In.m_bndCondExpansions),
640 m_globalBndMat(In.m_globalBndMat), m_traceMap(In.m_traceMap),
641 m_boundaryTraces(In.m_boundaryTraces),
642 m_periodicVerts(In.m_periodicVerts),
643 m_periodicFwdCopy(In.m_periodicFwdCopy),
644 m_periodicBwdCopy(In.m_periodicBwdCopy),
645 m_leftAdjacentTraces(In.m_leftAdjacentTraces),
646 m_locTraceToTraceMap(In.m_locTraceToTraceMap)
647{
648 if (In.m_trace)
649 {
651 *In.m_trace, DeclareCoeffPhysArrays);
652 }
653}
654
655/*
656 * @brief Copy type constructor which declares new boundary conditions
657 * and re-uses mapping info and trace space if possible
658 */
661 const std::string &variable, const bool SetUpJustDG,
662 const bool DeclareCoeffPhysArrays)
663 : ExpList(In, DeclareCoeffPhysArrays)
664{
665
667
668 // Set up boundary conditions for this variable.
669 // Do not set up BCs if default variable
670 if (variable.compare("DefaultVar") != 0)
671 {
673 GenerateBoundaryConditionExpansion(graph, bcs, variable);
674
675 if (DeclareCoeffPhysArrays)
676 {
677 EvaluateBoundaryConditions(0.0, variable);
678 }
679
681 {
682 // Find periodic edges for this variable.
683 FindPeriodicTraces(bcs, variable);
684
685 if (SetUpJustDG)
686 {
687 SetUpDG(variable);
688 }
689 else
690 {
691 // set elmt edges to point to robin bc edges if required
692 int i, cnt = 0;
693 Array<OneD, int> ElmtID, TraceID;
694 GetBoundaryToElmtMap(ElmtID, TraceID);
695
696 for (i = 0; i < m_bndCondExpansions.size(); ++i)
697 {
699
700 int e;
701 locExpList = m_bndCondExpansions[i];
702
703 for (e = 0; e < locExpList->GetExpSize(); ++e)
704 {
705 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
706 TraceID[cnt + e], locExpList->GetExp(e));
707 locExpList->GetExp(e)->SetAdjacentElementExp(
708 TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
709 }
710
711 cnt += m_bndCondExpansions[i]->GetExpSize();
712 }
713
714 if (m_session->DefinesSolverInfo("PROJECTION"))
715 {
716 std::string ProjectStr =
717 m_session->GetSolverInfo("PROJECTION");
718
719 if ((ProjectStr == "MixedCGDG") ||
720 (ProjectStr == "Mixed_CG_Discontinuous"))
721 {
722 SetUpDG(variable);
723 }
724 else
725 {
727 }
728 }
729 else
730 {
732 }
733 }
734 }
735 else
736 {
738 m_trace = In.m_trace;
749
750 if (!SetUpJustDG)
751 {
752 // set elmt edges to point to robin bc edges if required
753 int i, cnt = 0;
754 Array<OneD, int> ElmtID, TraceID;
755 GetBoundaryToElmtMap(ElmtID, TraceID);
756
757 for (i = 0; i < m_bndCondExpansions.size(); ++i)
758 {
760
761 int e;
762 locExpList = m_bndCondExpansions[i];
763
764 for (e = 0; e < locExpList->GetExpSize(); ++e)
765 {
766 (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
767 TraceID[cnt + e], locExpList->GetExp(e));
768 locExpList->GetExp(e)->SetAdjacentElementExp(
769 TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
770 }
771 cnt += m_bndCondExpansions[i]->GetExpSize();
772 }
773
775 }
776 }
777 }
778}
779
780/**
781 * Constructs a field as a copy of an existing explist field.
782 * @param In Existing ExpList object to copy.
783 */
785{
786}
787
788/**
789 *
790 */
792{
793}
794
795/**
796 * \brief This function discretises the boundary conditions by setting
797 * up a list of one-dimensions lower boundary expansions.
798 *
799 * According to their boundary region, the separate boundary
800 * expansions are bundled together in an object of the class
801 *
802 * @param graph A mesh, containing information about the domain
803 * and the spectral/hp element expansions.
804 * @param bcs Information about the enforced boundary
805 * conditions.
806 * @param variable The session variable associated with the
807 * boundary conditions to enforce.
808 * @param DeclareCoeffPhysArrays bool to identify if array
809 * space should be setup.
810 * Default is true.
811 */
814 const SpatialDomains::BoundaryConditions &bcs, const std::string variable,
815 const bool DeclareCoeffPhysArrays)
816{
817 int cnt = 0;
821 bcs.GetBoundaryRegions();
824
829
830 m_bndCondBndWeight = Array<OneD, NekDouble>{bregions.size(), 0.0};
831
832 // count the number of non-periodic boundary points
833 for (auto &it : bregions)
834 {
835 bc = GetBoundaryCondition(bconditions, it.first, variable);
836
838 m_session, *(it.second), graph, DeclareCoeffPhysArrays, variable,
839 false, bc->GetComm());
840
841 m_bndCondExpansions[cnt] = locExpList;
842 m_bndConditions[cnt] = bc;
843
844 std::string type = m_bndConditions[cnt]->GetUserDefined();
845
846 // Set up normals on non-Dirichlet boundary conditions. Second
847 // two conditions ideally should be in local solver setup (when
848 // made into factory)
849 if (bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet ||
850 boost::iequals(type, "I") || boost::iequals(type, "CalcBC"))
851 {
853 }
854 cnt++;
855 }
856}
857
858/**
859 * @brief Determine the periodic faces, edges and vertices for the given
860 * graph.
861 *
862 * @param bcs Information about the boundary conditions.
863 * @param variable Specifies the field.
864 */
866 const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
867{
869 bcs.GetBoundaryRegions();
872
873 LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
874
875 switch (m_expType)
876 {
877 case e1D:
878 {
879 int i, region1ID, region2ID;
880
882 map<int, int> BregionToVertMap;
883
884 // Construct list of all periodic Region and their
885 // global vertex on this process.
886 for (auto &it : bregions)
887 {
888 locBCond =
889 GetBoundaryCondition(bconditions, it.first, variable);
890
891 if (locBCond->GetBoundaryConditionType() !=
893 {
894 continue;
895 }
896 int id =
897 it.second->begin()->second->m_geomVec[0]->GetGlobalID();
898
899 BregionToVertMap[it.first] = id;
900 }
901
902 set<int> islocal;
903
904 int n = vComm->GetSize();
905 int p = vComm->GetRank();
906
907 Array<OneD, int> nregions(n, 0);
908 nregions[p] = BregionToVertMap.size();
909 vComm->AllReduce(nregions, LibUtilities::ReduceSum);
910
911 int totRegions = Vmath::Vsum(n, nregions, 1);
912
913 Array<OneD, int> regOffset(n, 0);
914
915 for (i = 1; i < n; ++i)
916 {
917 regOffset[i] = regOffset[i - 1] + nregions[i - 1];
918 }
919
920 Array<OneD, int> bregmap(totRegions, 0);
921 Array<OneD, int> bregid(totRegions, 0);
922
923 i = regOffset[p];
924 for (auto &iit : BregionToVertMap)
925 {
926 bregid[i] = iit.first;
927 bregmap[i++] = iit.second;
928 islocal.insert(iit.first);
929 }
930
931 vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
932 vComm->AllReduce(bregid, LibUtilities::ReduceSum);
933
934 for (int i = 0; i < totRegions; ++i)
935 {
936 BregionToVertMap[bregid[i]] = bregmap[i];
937 }
938
939 // Construct list of all periodic pairs local to this process.
940 for (auto &it : bregions)
941 {
942 locBCond =
943 GetBoundaryCondition(bconditions, it.first, variable);
944
945 if (locBCond->GetBoundaryConditionType() !=
947 {
948 continue;
949 }
950
951 // Identify periodic boundary region IDs.
952 region1ID = it.first;
953 region2ID =
954 std::static_pointer_cast<
956 ->m_connectedBoundaryRegion;
957
958 ASSERTL0(BregionToVertMap.count(region1ID) != 0,
959 "Cannot determine vertex of region1ID");
960
961 ASSERTL0(BregionToVertMap.count(region2ID) != 0,
962 "Cannot determine vertex of region2ID");
963
964 PeriodicEntity ent(BregionToVertMap[region2ID],
966 islocal.count(region2ID) != 0);
967
968 m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
969 }
970 }
971 break;
972 case e2D:
973 {
974 int region1ID, region2ID, i, j, k, cnt;
976
978 m_graph->GetCompositeOrdering();
980 m_graph->GetBndRegionOrdering();
981 SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
982
983 // Unique collection of pairs of periodic composites
984 // (i.e. if composites 1 and 2 are periodic then this
985 // map will contain either the pair (1,2) or (2,1) but
986 // not both).
987 map<int, int> perComps;
988 map<int, vector<int>> allVerts;
989 set<int> locVerts;
990 map<int, pair<int, StdRegions::Orientation>> allEdges;
991
992 // Set up a set of all local verts and edges.
993 for (i = 0; i < (*m_exp).size(); ++i)
994 {
995 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
996 {
997 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
998 locVerts.insert(id);
999 }
1000 }
1001
1002 // Construct list of all periodic pairs local to this process.
1003 for (auto &it : bregions)
1004 {
1005 locBCond =
1006 GetBoundaryCondition(bconditions, it.first, variable);
1007
1008 if (locBCond->GetBoundaryConditionType() !=
1010 {
1011 continue;
1012 }
1013
1014 // Identify periodic boundary region IDs.
1015 region1ID = it.first;
1016 region2ID =
1017 std::static_pointer_cast<
1019 ->m_connectedBoundaryRegion;
1020
1021 // From this identify composites. Note that in
1022 // serial this will be an empty map.
1023 int cId1, cId2;
1024 if (vComm->GetSize() == 1)
1025 {
1026 cId1 = it.second->begin()->first;
1027 cId2 = bregions.find(region2ID)->second->begin()->first;
1028 }
1029 else
1030 {
1031 cId1 = bndRegOrder.find(region1ID)->second[0];
1032 cId2 = bndRegOrder.find(region2ID)->second[0];
1033 }
1034
1035 ASSERTL0(it.second->size() == 1,
1036 "Boundary region " + std::to_string(region1ID) +
1037 " should only contain 1 composite.");
1038
1039 vector<unsigned int> tmpOrder;
1040
1041 // Construct set containing all periodic edgesd on
1042 // this process
1044 it.second->begin()->second;
1045
1046 for (i = 0; i < c->m_geomVec.size(); ++i)
1047 {
1049 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1050 c->m_geomVec[i]);
1051 ASSERTL0(segGeom, "Unable to cast to shared ptr");
1052
1054 m_graph->GetElementsFromEdge(segGeom);
1055 ASSERTL0(elmt->size() == 1,
1056 "The periodic boundaries belong to "
1057 "more than one element of the mesh");
1058
1060 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1061 elmt->at(0).first);
1062
1063 allEdges[c->m_geomVec[i]->GetGlobalID()] =
1064 make_pair(elmt->at(0).second,
1065 geom->GetEorient(elmt->at(0).second));
1066
1067 // In serial mesh partitioning will not have occurred so
1068 // need to fill composite ordering map manually.
1069 if (vComm->GetSize() == 1)
1070 {
1071 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1072 }
1073
1074 vector<int> vertList(2);
1075 vertList[0] = segGeom->GetVid(0);
1076 vertList[1] = segGeom->GetVid(1);
1077 allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
1078 }
1079
1080 if (vComm->GetSize() == 1)
1081 {
1082 compOrder[it.second->begin()->first] = tmpOrder;
1083 }
1084
1085 // See if we already have either region1 or
1086 // region2 stored in perComps map.
1087 if (perComps.count(cId1) == 0)
1088 {
1089 if (perComps.count(cId2) == 0)
1090 {
1091 perComps[cId1] = cId2;
1092 }
1093 else
1094 {
1095 std::stringstream ss;
1096 ss << "Boundary region " << cId2 << " should be "
1097 << "periodic with " << perComps[cId2] << " but "
1098 << "found " << cId1 << " instead!";
1099 ASSERTL0(perComps[cId2] == cId1, ss.str());
1100 }
1101 }
1102 else
1103 {
1104 std::stringstream ss;
1105 ss << "Boundary region " << cId1 << " should be "
1106 << "periodic with " << perComps[cId1] << " but "
1107 << "found " << cId2 << " instead!";
1108 ASSERTL0(perComps[cId1] == cId1, ss.str());
1109 }
1110 }
1111
1112 // In case of periodic partition being split over many composites
1113 // may not have all periodic matches so check this
1114 int idmax = -1;
1115 for (auto &cIt : perComps)
1116 {
1117 idmax = max(idmax, cIt.first);
1118 idmax = max(idmax, cIt.second);
1119 }
1120 vComm->AllReduce(idmax, LibUtilities::ReduceMax);
1121 idmax++;
1122 Array<OneD, int> perid(idmax, -1);
1123 for (auto &cIt : perComps)
1124 {
1125 perid[cIt.first] = cIt.second;
1126 }
1127 vComm->AllReduce(perid, LibUtilities::ReduceMax);
1128 // update all partitions
1129 for (int i = 0; i < idmax; ++i)
1130 {
1131 if (perid[i] > -1)
1132 {
1133 // skip if complemetary relationship has
1134 // already been speficied in list
1135 if (perComps.count(perid[i]))
1136 {
1137 continue;
1138 }
1139 perComps[i] = perid[i];
1140 }
1141 }
1142
1143 // Process local edge list to obtain relative edge
1144 // orientations.
1145 int n = vComm->GetSize();
1146 int p = vComm->GetRank();
1147 int totEdges;
1148 Array<OneD, int> edgecounts(n, 0);
1149 Array<OneD, int> edgeoffset(n, 0);
1150 Array<OneD, int> vertoffset(n, 0);
1151
1152 edgecounts[p] = allEdges.size();
1153 vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
1154
1155 edgeoffset[0] = 0;
1156 for (i = 1; i < n; ++i)
1157 {
1158 edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
1159 }
1160
1161 totEdges = Vmath::Vsum(n, edgecounts, 1);
1162 Array<OneD, int> edgeIds(totEdges, 0);
1163 Array<OneD, int> edgeIdx(totEdges, 0);
1164 Array<OneD, int> edgeOrient(totEdges, 0);
1165 Array<OneD, int> edgeVerts(totEdges, 0);
1166
1167 auto sIt = allEdges.begin();
1168
1169 for (i = 0; sIt != allEdges.end(); ++sIt)
1170 {
1171 edgeIds[edgeoffset[p] + i] = sIt->first;
1172 edgeIdx[edgeoffset[p] + i] = sIt->second.first;
1173 edgeOrient[edgeoffset[p] + i] = sIt->second.second;
1174 edgeVerts[edgeoffset[p] + i++] = allVerts[sIt->first].size();
1175 }
1176
1177 vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
1178 vComm->AllReduce(edgeIdx, LibUtilities::ReduceSum);
1179 vComm->AllReduce(edgeOrient, LibUtilities::ReduceSum);
1180 vComm->AllReduce(edgeVerts, LibUtilities::ReduceSum);
1181
1182 // Calculate number of vertices on each processor.
1183 Array<OneD, int> procVerts(n, 0);
1184 int nTotVerts;
1185
1186 // Note if there are no periodic edges at all calling Vsum will
1187 // cause a segfault.
1188 if (totEdges > 0)
1189 {
1190 nTotVerts = Vmath::Vsum(totEdges, edgeVerts, 1);
1191 }
1192 else
1193 {
1194 nTotVerts = 0;
1195 }
1196
1197 for (i = 0; i < n; ++i)
1198 {
1199 if (edgecounts[i] > 0)
1200 {
1201 procVerts[i] = Vmath::Vsum(edgecounts[i],
1202 edgeVerts + edgeoffset[i], 1);
1203 }
1204 else
1205 {
1206 procVerts[i] = 0;
1207 }
1208 }
1209 vertoffset[0] = 0;
1210
1211 for (i = 1; i < n; ++i)
1212 {
1213 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1214 }
1215
1216 Array<OneD, int> traceIds(nTotVerts, 0);
1217 for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
1218 {
1219 for (j = 0; j < allVerts[sIt->first].size(); ++j)
1220 {
1221 traceIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
1222 }
1223 }
1224
1225 vComm->AllReduce(traceIds, LibUtilities::ReduceSum);
1226
1227 // For simplicity's sake create a map of edge id -> orientation.
1228 map<int, StdRegions::Orientation> orientMap;
1229 map<int, vector<int>> vertMap;
1230
1231 for (cnt = i = 0; i < totEdges; ++i)
1232 {
1233 ASSERTL0(orientMap.count(edgeIds[i]) == 0,
1234 "Already found edge in orientation map!");
1235
1236 // Work out relative orientations. To avoid having
1237 // to exchange vertex locations, we figure out if
1238 // the edges are backwards or forwards orientated
1239 // with respect to a forwards orientation that is
1240 // CCW. Since our local geometries are
1241 // forwards-orientated with respect to the
1242 // Cartesian axes, we need to invert the
1243 // orientation for the top and left edges of a
1244 // quad and the left edge of a triangle.
1246 (StdRegions::Orientation)edgeOrient[i];
1247
1248 if (edgeIdx[i] > 1)
1249 {
1252 }
1253
1254 orientMap[edgeIds[i]] = o;
1255
1256 vector<int> verts(edgeVerts[i]);
1257
1258 for (j = 0; j < edgeVerts[i]; ++j)
1259 {
1260 verts[j] = traceIds[cnt++];
1261 }
1262 vertMap[edgeIds[i]] = verts;
1263 }
1264
1265 // Go through list of composites and figure out which
1266 // edges are parallel from original ordering in
1267 // session file. This includes composites which are
1268 // not necessarily on this process.
1269 map<int, int> allCompPairs;
1270
1271 // Store temporary map of periodic vertices which will hold all
1272 // periodic vertices on the entire mesh so that doubly periodic
1273 // vertices can be counted properly across partitions. Local
1274 // vertices are copied into m_periodicVerts at the end of the
1275 // function.
1276 PeriodicMap periodicVerts;
1277
1278 for (auto &cIt : perComps)
1279 {
1281 const int id1 = cIt.first;
1282 const int id2 = cIt.second;
1283 std::string id1s = std::to_string(id1);
1284 std::string id2s = std::to_string(id2);
1285
1286 if (compMap.count(id1) > 0)
1287 {
1288 c[0] = compMap[id1];
1289 }
1290
1291 if (compMap.count(id2) > 0)
1292 {
1293 c[1] = compMap[id2];
1294 }
1295
1296 // Loop over composite ordering to construct list of all
1297 // periodic edges regardless of whether they are on this
1298 // process.
1299 map<int, int> compPairs;
1300
1301 ASSERTL0(compOrder.count(id1) > 0,
1302 "Unable to find composite " + id1s + " in order map.");
1303 ASSERTL0(compOrder.count(id2) > 0,
1304 "Unable to find composite " + id2s + " in order map.");
1305 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1306 "Periodic composites " + id1s + " and " + id2s +
1307 " should have the same number of elements.");
1308 ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
1309 id1s + " and " + id2s +
1310 " are empty!");
1311
1312 // TODO: Add more checks.
1313 for (i = 0; i < compOrder[id1].size(); ++i)
1314 {
1315 int eId1 = compOrder[id1][i];
1316 int eId2 = compOrder[id2][i];
1317
1318 ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
1319
1320 if (compPairs.count(eId2) != 0)
1321 {
1322 ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
1323 }
1324 compPairs[eId1] = eId2;
1325 }
1326
1327 // Construct set of all edges that we have
1328 // locally on this processor.
1329 set<int> locEdges;
1330 for (i = 0; i < 2; ++i)
1331 {
1332 if (!c[i])
1333 {
1334 continue;
1335 }
1336
1337 if (c[i]->m_geomVec.size() > 0)
1338 {
1339 for (j = 0; j < c[i]->m_geomVec.size(); ++j)
1340 {
1341 locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
1342 }
1343 }
1344 }
1345
1346 // Loop over all edges in the geometry composite.
1347 for (auto &pIt : compPairs)
1348 {
1349 int ids[2] = {pIt.first, pIt.second};
1350 bool local[2] = {locEdges.count(pIt.first) > 0,
1351 locEdges.count(pIt.second) > 0};
1352
1353 ASSERTL0(orientMap.count(ids[0]) > 0 &&
1354 orientMap.count(ids[1]) > 0,
1355 "Unable to find edge in orientation map.");
1356
1357 allCompPairs[pIt.first] = pIt.second;
1358 allCompPairs[pIt.second] = pIt.first;
1359
1360 for (i = 0; i < 2; ++i)
1361 {
1362 if (!local[i])
1363 {
1364 continue;
1365 }
1366
1367 int other = (i + 1) % 2;
1368
1370 orientMap[ids[i]] == orientMap[ids[other]]
1373
1374 PeriodicEntity ent(ids[other], o, local[other]);
1375 m_periodicEdges[ids[i]].push_back(ent);
1376 }
1377
1378 for (i = 0; i < 2; ++i)
1379 {
1380 int other = (i + 1) % 2;
1381
1383 orientMap[ids[i]] == orientMap[ids[other]]
1386
1387 // Determine periodic vertices.
1388 vector<int> perVerts1 = vertMap[ids[i]];
1389 vector<int> perVerts2 = vertMap[ids[other]];
1390
1391 map<int, pair<int, bool>> tmpMap;
1392 if (o == StdRegions::eForwards)
1393 {
1394 tmpMap[perVerts1[0]] = make_pair(
1395 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1396 tmpMap[perVerts1[1]] = make_pair(
1397 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1398 }
1399 else
1400 {
1401 tmpMap[perVerts1[0]] = make_pair(
1402 perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1403 tmpMap[perVerts1[1]] = make_pair(
1404 perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1405 }
1406
1407 for (auto &mIt : tmpMap)
1408 {
1409 // See if this vertex has been recorded already.
1410 PeriodicEntity ent2(mIt.second.first,
1412 mIt.second.second);
1413 auto perIt = periodicVerts.find(mIt.first);
1414
1415 if (perIt == periodicVerts.end())
1416 {
1417 periodicVerts[mIt.first].push_back(ent2);
1418 perIt = periodicVerts.find(mIt.first);
1419 }
1420 else
1421 {
1422 bool doAdd = true;
1423 for (j = 0; j < perIt->second.size(); ++j)
1424 {
1425 if (perIt->second[j].id == mIt.second.first)
1426 {
1427 doAdd = false;
1428 break;
1429 }
1430 }
1431
1432 if (doAdd)
1433 {
1434 perIt->second.push_back(ent2);
1435 }
1436 }
1437 }
1438 }
1439 }
1440 }
1441
1442 // Search for periodic vertices and edges which are not in
1443 // a periodic composite but lie in this process. First, loop
1444 // over all information we have from other processors.
1445 for (cnt = i = 0; i < totEdges; ++i)
1446 {
1447 int edgeId = edgeIds[i];
1448
1449 ASSERTL0(allCompPairs.count(edgeId) > 0,
1450 "Unable to find matching periodic edge.");
1451
1452 int perEdgeId = allCompPairs[edgeId];
1453
1454 for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1455 {
1456 int vId = traceIds[cnt];
1457
1458 auto perId = periodicVerts.find(vId);
1459
1460 if (perId == periodicVerts.end())
1461 {
1462 // This vertex is not included in the
1463 // map. Figure out which vertex it is
1464 // supposed to be periodic with. perEdgeId
1465 // is the edge ID which is periodic with
1466 // edgeId. The logic is much the same as
1467 // the loop above.
1468 int perVertexId =
1469 orientMap[edgeId] == orientMap[perEdgeId]
1470 ? vertMap[perEdgeId][(j + 1) % 2]
1471 : vertMap[perEdgeId][j];
1472
1473 PeriodicEntity ent(perVertexId,
1475 locVerts.count(perVertexId) > 0);
1476
1477 periodicVerts[vId].push_back(ent);
1478 }
1479 }
1480 }
1481
1482 // Loop over all periodic vertices to complete connectivity
1483 // information.
1484 for (auto &perIt : periodicVerts)
1485 {
1486 // Loop over associated vertices.
1487 for (i = 0; i < perIt.second.size(); ++i)
1488 {
1489 auto perIt2 = periodicVerts.find(perIt.second[i].id);
1490 ASSERTL0(perIt2 != periodicVerts.end(),
1491 "Couldn't find periodic vertex.");
1492
1493 for (j = 0; j < perIt2->second.size(); ++j)
1494 {
1495 if (perIt2->second[j].id == perIt.first)
1496 {
1497 continue;
1498 }
1499
1500 bool doAdd = true;
1501 for (k = 0; k < perIt.second.size(); ++k)
1502 {
1503 if (perIt2->second[j].id == perIt.second[k].id)
1504 {
1505 doAdd = false;
1506 break;
1507 }
1508 }
1509
1510 if (doAdd)
1511 {
1512 perIt.second.push_back(perIt2->second[j]);
1513 }
1514 }
1515 }
1516 }
1517
1518 // Do one final loop over periodic vertices to remove non-local
1519 // vertices from map.
1520 for (auto &perIt : periodicVerts)
1521 {
1522 if (locVerts.count(perIt.first) > 0)
1523 {
1524 m_periodicVerts.insert(perIt);
1525 }
1526 }
1527 }
1528 break;
1529 case e3D:
1530 {
1532 m_graph->GetCompositeOrdering();
1534 m_graph->GetBndRegionOrdering();
1535 SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
1536
1537 // perComps: Stores a unique collection of pairs of periodic
1538 // composites (i.e. if composites 1 and 2 are periodic then this map
1539 // will contain either the pair (1,2) or (2,1) but not both).
1540 //
1541 // The four maps allVerts, allCoord, allEdges and allOrient map a
1542 // periodic face to a vector containing the vertex ids of the face;
1543 // their coordinates; the edge ids of the face; and their
1544 // orientation within that face respectively.
1545 //
1546 // Finally the three sets locVerts, locEdges and locFaces store any
1547 // vertices, edges and faces that belong to a periodic composite and
1548 // lie on this process.
1549 map<int, RotPeriodicInfo> rotComp;
1550 map<int, int> perComps;
1551 map<int, vector<int>> allVerts;
1552 map<int, SpatialDomains::PointGeomVector> allCoord;
1553 map<int, vector<int>> allEdges;
1554 map<int, vector<StdRegions::Orientation>> allOrient;
1555 set<int> locVerts;
1556 set<int> locEdges;
1557 set<int> locFaces;
1558
1559 int region1ID, region2ID, i, j, k, cnt;
1561
1562 // Set up a set of all local verts and edges.
1563 for (i = 0; i < (*m_exp).size(); ++i)
1564 {
1565 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
1566 {
1567 int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1568 locVerts.insert(id);
1569 }
1570
1571 for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
1572 {
1573 int id = (*m_exp)[i]->GetGeom()->GetEid(j);
1574 locEdges.insert(id);
1575 }
1576 }
1577
1578 // Begin by populating the perComps map. We loop over all periodic
1579 // boundary conditions and determine the composite associated with
1580 // it, then fill out the all* maps.
1581 for (auto &it : bregions)
1582 {
1583
1584 locBCond =
1585 GetBoundaryCondition(bconditions, it.first, variable);
1586
1587 if (locBCond->GetBoundaryConditionType() !=
1589 {
1590 continue;
1591 }
1592
1593 // Identify periodic boundary region IDs.
1594 region1ID = it.first;
1595 region2ID =
1596 std::static_pointer_cast<
1598 ->m_connectedBoundaryRegion;
1599
1600 // Check the region only contains a single composite.
1601 ASSERTL0(it.second->size() == 1,
1602 "Boundary region " + std::to_string(region1ID) +
1603 " should only contain 1 composite.");
1604
1605 // From this identify composites by looking at the original
1606 // boundary region ordering. Note that in serial the mesh
1607 // partitioner is not run, so this map will be empty and
1608 // therefore needs to be populated by using the corresponding
1609 // boundary region.
1610 int cId1, cId2;
1611 if (vComm->GetSize() == 1)
1612 {
1613 cId1 = it.second->begin()->first;
1614 cId2 = bregions.find(region2ID)->second->begin()->first;
1615 }
1616 else
1617 {
1618 cId1 = bndRegOrder.find(region1ID)->second[0];
1619 cId2 = bndRegOrder.find(region2ID)->second[0];
1620 }
1621
1622 // check to see if boundary is rotationally aligned
1623 if (boost::icontains(locBCond->GetUserDefined(), "Rotated"))
1624 {
1625 vector<string> tmpstr;
1626
1627 boost::split(tmpstr, locBCond->GetUserDefined(),
1628 boost::is_any_of(":"));
1629
1630 if (boost::iequals(tmpstr[0], "Rotated"))
1631 {
1632 ASSERTL1(tmpstr.size() > 2,
1633 "Expected Rotated user defined string to "
1634 "contain direction and rotation angle "
1635 "and optionally a tolerance, "
1636 "i.e. Rotated:dir:PI/2:1e-6");
1637
1638 ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
1639 (tmpstr[1] == "z"),
1640 "Rotated Dir is "
1641 "not specified as x,y or z");
1642
1643 RotPeriodicInfo RotInfo;
1644 RotInfo.m_dir = (tmpstr[1] == "x") ? 0
1645 : (tmpstr[1] == "y") ? 1
1646 : 2;
1647
1649 int ExprId = strEval.DefineFunction("", tmpstr[2]);
1650 RotInfo.m_angle = strEval.Evaluate(ExprId);
1651
1652 if (tmpstr.size() == 4)
1653 {
1654 try
1655 {
1656 RotInfo.m_tol = std::stod(tmpstr[3]);
1657 }
1658 catch (...)
1659 {
1661 "failed to cast tolerance input "
1662 "to a double value in Rotated"
1663 "boundary information");
1664 }
1665 }
1666 else
1667 {
1668 RotInfo.m_tol = 1e-8;
1669 }
1670 rotComp[cId1] = RotInfo;
1671 }
1672 }
1673
1675 it.second->begin()->second;
1676
1677 vector<unsigned int> tmpOrder;
1678
1679 // store the rotation info of this
1680
1681 // From the composite, we now construct the allVerts, allEdges
1682 // and allCoord map so that they can be transferred across
1683 // processors. We also populate the locFaces set to store a
1684 // record of all faces local to this process.
1685 for (i = 0; i < c->m_geomVec.size(); ++i)
1686 {
1688 std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1689 c->m_geomVec[i]);
1690 ASSERTL1(faceGeom, "Unable to cast to shared ptr");
1691
1692 // Get geometry ID of this face and store in locFaces.
1693 int faceId = c->m_geomVec[i]->GetGlobalID();
1694 locFaces.insert(faceId);
1695
1696 // In serial, mesh partitioning will not have occurred so
1697 // need to fill composite ordering map manually.
1698 if (vComm->GetSize() == 1)
1699 {
1700 tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1701 }
1702
1703 // Loop over vertices and edges of the face to populate
1704 // allVerts, allEdges and allCoord maps.
1705 vector<int> vertList, edgeList;
1707 vector<StdRegions::Orientation> orientVec;
1708 for (j = 0; j < faceGeom->GetNumVerts(); ++j)
1709 {
1710 vertList.push_back(faceGeom->GetVid(j));
1711 edgeList.push_back(faceGeom->GetEid(j));
1712 coordVec.push_back(faceGeom->GetVertex(j));
1713 orientVec.push_back(faceGeom->GetEorient(j));
1714 }
1715
1716 allVerts[faceId] = vertList;
1717 allEdges[faceId] = edgeList;
1718 allCoord[faceId] = coordVec;
1719 allOrient[faceId] = orientVec;
1720 }
1721
1722 // In serial, record the composite ordering in compOrder for
1723 // later in the routine.
1724 if (vComm->GetSize() == 1)
1725 {
1726 compOrder[it.second->begin()->first] = tmpOrder;
1727 }
1728
1729 // See if we already have either region1 or region2 stored in
1730 // perComps map already and do a sanity check to ensure regions
1731 // are mutually periodic.
1732 if (perComps.count(cId1) == 0)
1733 {
1734 if (perComps.count(cId2) == 0)
1735 {
1736 perComps[cId1] = cId2;
1737 }
1738 else
1739 {
1740 std::stringstream ss;
1741 ss << "Boundary region " << cId2 << " should be "
1742 << "periodic with " << perComps[cId2] << " but "
1743 << "found " << cId1 << " instead!";
1744 ASSERTL0(perComps[cId2] == cId1, ss.str());
1745 }
1746 }
1747 else
1748 {
1749 std::stringstream ss;
1750 ss << "Boundary region " << cId1 << " should be "
1751 << "periodic with " << perComps[cId1] << " but "
1752 << "found " << cId2 << " instead!";
1753 ASSERTL0(perComps[cId1] == cId1, ss.str());
1754 }
1755 }
1756
1757 // In case of periodic partition being split over many composites
1758 // may not have all periodic matches so check this
1759 int idmax = -1;
1760 for (auto &cIt : perComps)
1761 {
1762 idmax = max(idmax, cIt.first);
1763 idmax = max(idmax, cIt.second);
1764 }
1765 vComm->AllReduce(idmax, LibUtilities::ReduceMax);
1766 idmax++;
1767 Array<OneD, int> perid(idmax, -1);
1768 for (auto &cIt : perComps)
1769 {
1770 perid[cIt.first] = cIt.second;
1771 }
1772 vComm->AllReduce(perid, LibUtilities::ReduceMax);
1773 // update all partitions
1774 for (int i = 0; i < idmax; ++i)
1775 {
1776 if (perid[i] > -1)
1777 {
1778 // skip if equivlaent relationship has
1779 // already been speficied in list
1780 if (perComps.count(perid[i]))
1781 {
1782 continue;
1783 }
1784 perComps[i] = perid[i];
1785 }
1786 }
1787
1788 // The next routines process local face lists to
1789 // exchange vertices,
1790 // edges and faces.
1791 int n = vComm->GetSize();
1792 int p = vComm->GetRank();
1793 int totFaces;
1794 Array<OneD, int> facecounts(n, 0);
1795 Array<OneD, int> vertcounts(n, 0);
1796 Array<OneD, int> faceoffset(n, 0);
1797 Array<OneD, int> vertoffset(n, 0);
1798
1799 Array<OneD, int> rotcounts(n, 0);
1800 Array<OneD, int> rotoffset(n, 0);
1801
1802 rotcounts[p] = rotComp.size();
1803 vComm->AllReduce(rotcounts, LibUtilities::ReduceSum);
1804 int totrot = Vmath::Vsum(n, rotcounts, 1);
1805
1806 if (totrot)
1807 {
1808 for (i = 1; i < n; ++i)
1809 {
1810 rotoffset[i] = rotoffset[i - 1] + rotcounts[i - 1];
1811 }
1812
1813 Array<OneD, int> compid(totrot, 0);
1814 Array<OneD, int> rotdir(totrot, 0);
1815 Array<OneD, NekDouble> rotangle(totrot, 0.0);
1816 Array<OneD, NekDouble> rottol(totrot, 0.0);
1817
1818 // fill in rotational informaiton
1819 auto rIt = rotComp.begin();
1820
1821 for (i = 0; rIt != rotComp.end(); ++rIt)
1822 {
1823 compid[rotoffset[p] + i] = rIt->first;
1824 rotdir[rotoffset[p] + i] = rIt->second.m_dir;
1825 rotangle[rotoffset[p] + i] = rIt->second.m_angle;
1826 rottol[rotoffset[p] + i++] = rIt->second.m_tol;
1827 }
1828
1829 vComm->AllReduce(compid, LibUtilities::ReduceSum);
1830 vComm->AllReduce(rotdir, LibUtilities::ReduceSum);
1831 vComm->AllReduce(rotangle, LibUtilities::ReduceSum);
1832 vComm->AllReduce(rottol, LibUtilities::ReduceSum);
1833
1834 // Fill in full rotational composite list
1835 for (i = 0; i < totrot; ++i)
1836 {
1837 RotPeriodicInfo rinfo(rotdir[i], rotangle[i], rottol[i]);
1838
1839 rotComp[compid[i]] = rinfo;
1840 }
1841 }
1842
1843 // First exchange the number of faces on each process.
1844 facecounts[p] = locFaces.size();
1845 vComm->AllReduce(facecounts, LibUtilities::ReduceSum);
1846
1847 // Set up an offset map to allow us to distribute face IDs to all
1848 // processors.
1849 faceoffset[0] = 0;
1850 for (i = 1; i < n; ++i)
1851 {
1852 faceoffset[i] = faceoffset[i - 1] + facecounts[i - 1];
1853 }
1854
1855 // Calculate total number of faces.
1856 totFaces = Vmath::Vsum(n, facecounts, 1);
1857
1858 // faceIds holds face IDs for each periodic face. faceVerts holds
1859 // the number of vertices in this face.
1860 Array<OneD, int> faceIds(totFaces, 0);
1861 Array<OneD, int> faceVerts(totFaces, 0);
1862
1863 // Process p writes IDs of its faces into position faceoffset[p] of
1864 // faceIds which allows us to perform an AllReduce to distribute
1865 // information amongst processors.
1866 auto sIt = locFaces.begin();
1867 for (i = 0; sIt != locFaces.end(); ++sIt)
1868 {
1869 faceIds[faceoffset[p] + i] = *sIt;
1870 faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
1871 }
1872
1873 vComm->AllReduce(faceIds, LibUtilities::ReduceSum);
1874 vComm->AllReduce(faceVerts, LibUtilities::ReduceSum);
1875
1876 // procVerts holds number of vertices (and also edges since each
1877 // face is 2D) on each process.
1878 Array<OneD, int> procVerts(n, 0);
1879 int nTotVerts;
1880
1881 // Note if there are no periodic faces at all calling Vsum will
1882 // cause a segfault.
1883 if (totFaces > 0)
1884 {
1885 // Calculate number of vertices on each processor.
1886 nTotVerts = Vmath::Vsum(totFaces, faceVerts, 1);
1887 }
1888 else
1889 {
1890 nTotVerts = 0;
1891 }
1892
1893 for (i = 0; i < n; ++i)
1894 {
1895 if (facecounts[i] > 0)
1896 {
1897 procVerts[i] = Vmath::Vsum(facecounts[i],
1898 faceVerts + faceoffset[i], 1);
1899 }
1900 else
1901 {
1902 procVerts[i] = 0;
1903 }
1904 }
1905
1906 // vertoffset is defined in the same manner as edgeoffset
1907 // beforehand.
1908 vertoffset[0] = 0;
1909 for (i = 1; i < n; ++i)
1910 {
1911 vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1912 }
1913
1914 // At this point we exchange all vertex IDs, edge IDs and vertex
1915 // coordinates for each face. The coordinates are necessary because
1916 // we need to calculate relative face orientations between periodic
1917 // faces to determined edge and vertex connectivity.
1918 Array<OneD, int> vertIds(nTotVerts, 0);
1919 Array<OneD, int> edgeIds(nTotVerts, 0);
1920 Array<OneD, int> edgeOrt(nTotVerts, 0);
1921 Array<OneD, NekDouble> vertX(nTotVerts, 0.0);
1922 Array<OneD, NekDouble> vertY(nTotVerts, 0.0);
1923 Array<OneD, NekDouble> vertZ(nTotVerts, 0.0);
1924
1925 for (cnt = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
1926 {
1927 for (j = 0; j < allVerts[*sIt].size(); ++j)
1928 {
1929 int vertId = allVerts[*sIt][j];
1930 vertIds[vertoffset[p] + cnt] = vertId;
1931 vertX[vertoffset[p] + cnt] = (*allCoord[*sIt][j])(0);
1932 vertY[vertoffset[p] + cnt] = (*allCoord[*sIt][j])(1);
1933 vertZ[vertoffset[p] + cnt] = (*allCoord[*sIt][j])(2);
1934 edgeIds[vertoffset[p] + cnt] = allEdges[*sIt][j];
1935 edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
1936 }
1937 }
1938
1939 vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
1940 vComm->AllReduce(vertX, LibUtilities::ReduceSum);
1941 vComm->AllReduce(vertY, LibUtilities::ReduceSum);
1942 vComm->AllReduce(vertZ, LibUtilities::ReduceSum);
1943 vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
1944 vComm->AllReduce(edgeOrt, LibUtilities::ReduceSum);
1945
1946 // Finally now we have all of this information, we construct maps
1947 // which make accessing the information easier. These are
1948 // conceptually the same as all* maps at the beginning of the
1949 // routine, but now hold information for all periodic vertices.
1950 map<int, vector<int>> vertMap;
1951 map<int, vector<int>> edgeMap;
1952 map<int, SpatialDomains::PointGeomVector> coordMap;
1953
1954 // These final two maps are required for determining the relative
1955 // orientation of periodic edges. vCoMap associates vertex IDs with
1956 // their coordinates, and eIdMap maps an edge ID to the two vertices
1957 // which construct it.
1958 map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
1959 map<int, pair<int, int>> eIdMap;
1960
1961 for (cnt = i = 0; i < totFaces; ++i)
1962 {
1963 vector<int> edges(faceVerts[i]);
1964 vector<int> verts(faceVerts[i]);
1965 SpatialDomains::PointGeomVector coord(faceVerts[i]);
1966
1967 // Keep track of cnt to enable correct edge vertices to be
1968 // inserted into eIdMap.
1969 int tmp = cnt;
1970 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1971 {
1972 edges[j] = edgeIds[cnt];
1973 verts[j] = vertIds[cnt];
1975 AllocateSharedPtr(3, verts[j], vertX[cnt], vertY[cnt],
1976 vertZ[cnt]);
1977 vCoMap[vertIds[cnt]] = coord[j];
1978
1979 // Try to insert edge into the eIdMap to avoid re-inserting.
1980 auto testIns = eIdMap.insert(make_pair(
1981 edgeIds[cnt],
1982 make_pair(vertIds[tmp + j],
1983 vertIds[tmp + ((j + 1) % faceVerts[i])])));
1984
1985 if (testIns.second == false)
1986 {
1987 continue;
1988 }
1989
1990 // If the edge is reversed with respect to the face, then
1991 // swap the edges so that we have the original ordering of
1992 // the edge in the 3D element. This is necessary to properly
1993 // determine edge orientation. Note that the logic relies on
1994 // the fact that all edge forward directions are CCW
1995 // orientated: we use a tensor product ordering for 2D
1996 // elements so need to reverse this for edge IDs 2 and 3.
1997 StdRegions::Orientation edgeOrient =
1998 static_cast<StdRegions::Orientation>(edgeOrt[cnt]);
1999 if (j > 1)
2000 {
2001 edgeOrient = edgeOrient == StdRegions::eForwards
2004 }
2005
2006 if (edgeOrient == StdRegions::eBackwards)
2007 {
2008 swap(testIns.first->second.first,
2009 testIns.first->second.second);
2010 }
2011 }
2012
2013 vertMap[faceIds[i]] = verts;
2014 edgeMap[faceIds[i]] = edges;
2015 coordMap[faceIds[i]] = coord;
2016 }
2017
2018 // Go through list of composites and figure out which edges are
2019 // parallel from original ordering in session file. This includes
2020 // composites which are not necessarily on this process.
2021
2022 // Store temporary map of periodic vertices which will hold all
2023 // periodic vertices on the entire mesh so that doubly periodic
2024 // vertices/edges can be counted properly across partitions. Local
2025 // vertices/edges are copied into m_periodicVerts and
2026 // m_periodicEdges at the end of the function.
2027 PeriodicMap periodicVerts, periodicEdges;
2028
2029 // Construct two maps which determine how vertices and edges of
2030 // faces connect given a specific face orientation. The key of the
2031 // map is the number of vertices in the face, used to determine
2032 // difference between tris and quads.
2033 map<int, map<StdRegions::Orientation, vector<int>>> vmap;
2034 map<int, map<StdRegions::Orientation, vector<int>>> emap;
2035
2036 map<StdRegions::Orientation, vector<int>> quadVertMap;
2037 quadVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
2038 quadVertMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {3, 2, 1, 0};
2039 quadVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 3, 2};
2040 quadVertMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
2041 quadVertMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {0, 3, 2, 1};
2042 quadVertMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
2043 quadVertMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
2044 quadVertMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {2, 1, 0, 3};
2045
2046 map<StdRegions::Orientation, vector<int>> quadEdgeMap;
2047 quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
2048 quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {2, 1, 0, 3};
2049 quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 3, 2, 1};
2050 quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
2051 quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {3, 2, 1, 0};
2052 quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
2053 quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
2054 quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {1, 0, 3, 2};
2055
2056 map<StdRegions::Orientation, vector<int>> triVertMap;
2057 triVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
2058 triVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 2};
2059
2060 map<StdRegions::Orientation, vector<int>> triEdgeMap;
2061 triEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
2062 triEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 2, 1};
2063
2064 vmap[3] = triVertMap;
2065 vmap[4] = quadVertMap;
2066 emap[3] = triEdgeMap;
2067 emap[4] = quadEdgeMap;
2068
2069 map<int, int> allCompPairs;
2070
2071 // Collect composite id's of each periodic face for use if rotation
2072 // is required
2073 map<int, int> fIdToCompId;
2074
2075 // Finally we have enough information to populate the periodic
2076 // vertex, edge and face maps. Begin by looping over all pairs of
2077 // periodic composites to determine pairs of periodic faces.
2078 for (auto &cIt : perComps)
2079 {
2081 const int id1 = cIt.first;
2082 const int id2 = cIt.second;
2083 std::string id1s = std::to_string(id1);
2084 std::string id2s = std::to_string(id2);
2085
2086 if (compMap.count(id1) > 0)
2087 {
2088 c[0] = compMap[id1];
2089 }
2090
2091 if (compMap.count(id2) > 0)
2092 {
2093 c[1] = compMap[id2];
2094 }
2095
2096 // Loop over composite ordering to construct list of all
2097 // periodic faces, regardless of whether they are on this
2098 // process.
2099 map<int, int> compPairs;
2100
2101 ASSERTL0(compOrder.count(id1) > 0,
2102 "Unable to find composite " + id1s + " in order map.");
2103 ASSERTL0(compOrder.count(id2) > 0,
2104 "Unable to find composite " + id2s + " in order map.");
2105 ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
2106 "Periodic composites " + id1s + " and " + id2s +
2107 " should have the same number of elements.");
2108 ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
2109 id1s + " and " + id2s +
2110 " are empty!");
2111
2112 // Look up composite ordering to determine pairs.
2113 for (i = 0; i < compOrder[id1].size(); ++i)
2114 {
2115 int eId1 = compOrder[id1][i];
2116 int eId2 = compOrder[id2][i];
2117
2118 ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
2119
2120 // Sanity check that the faces are mutually periodic.
2121 if (compPairs.count(eId2) != 0)
2122 {
2123 ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
2124 }
2125 compPairs[eId1] = eId2;
2126
2127 // store a map of face ids to composite ids
2128 fIdToCompId[eId1] = id1;
2129 fIdToCompId[eId2] = id2;
2130 }
2131
2132 // Now that we have all pairs of periodic faces, loop over the
2133 // ones local on this process and populate face/edge/vertex
2134 // maps.
2135 for (auto &pIt : compPairs)
2136 {
2137 int ids[2] = {pIt.first, pIt.second};
2138 bool local[2] = {locFaces.count(pIt.first) > 0,
2139 locFaces.count(pIt.second) > 0};
2140
2141 ASSERTL0(coordMap.count(ids[0]) > 0 &&
2142 coordMap.count(ids[1]) > 0,
2143 "Unable to find face in coordinate map");
2144
2145 allCompPairs[pIt.first] = pIt.second;
2146 allCompPairs[pIt.second] = pIt.first;
2147
2148 // Loop up coordinates of the faces, check they have the
2149 // same number of vertices.
2151 coordMap[ids[0]], coordMap[ids[1]]};
2152
2153 ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
2154 "Two periodic faces have different number "
2155 "of vertices!");
2156
2157 // o will store relative orientation of faces. Note that in
2158 // some transpose cases (Dir1FwdDir2_Dir2BwdDir1 and
2159 // Dir1BwdDir1_Dir2FwdDir1) it seems orientation will be
2160 // different going from face1->face2 instead of face2->face1
2161 // (check this).
2163 bool rotbnd = false;
2164 int dir = 0;
2165 NekDouble angle = 0.0;
2166 NekDouble sign = 0.0;
2167 NekDouble tol = 1e-8;
2168
2169 // check to see if perioid boundary is rotated
2170 if (rotComp.count(fIdToCompId[pIt.first]))
2171 {
2172 rotbnd = true;
2173 dir = rotComp[fIdToCompId[pIt.first]].m_dir;
2174 angle = rotComp[fIdToCompId[pIt.first]].m_angle;
2175 tol = rotComp[fIdToCompId[pIt.first]].m_tol;
2176 }
2177
2178 // Record periodic faces.
2179 for (i = 0; i < 2; ++i)
2180 {
2181 if (!local[i])
2182 {
2183 continue;
2184 }
2185
2186 // Reference to the other face.
2187 int other = (i + 1) % 2;
2188
2189 // angle is set up for i = 0 to i = 1
2190 sign = (i == 0) ? 1.0 : -1.0;
2191
2192 // Calculate relative face orientation.
2193 if (tmpVec[0].size() == 3)
2194 {
2196 tmpVec[i], tmpVec[other], rotbnd, dir,
2197 sign * angle, tol);
2198 }
2199 else
2200 {
2202 tmpVec[i], tmpVec[other], rotbnd, dir,
2203 sign * angle, tol);
2204 }
2205
2206 // Record face ID, orientation and whether other face is
2207 // local.
2208 PeriodicEntity ent(ids[other], o, local[other]);
2209 m_periodicFaces[ids[i]].push_back(ent);
2210 }
2211
2212 int nFaceVerts = vertMap[ids[0]].size();
2213
2214 // Determine periodic vertices.
2215 for (i = 0; i < 2; ++i)
2216 {
2217 int other = (i + 1) % 2;
2218
2219 // angle is set up for i = 0 to i = 1
2220 sign = (i == 0) ? 1.0 : -1.0;
2221
2222 // Calculate relative face orientation.
2223 if (tmpVec[0].size() == 3)
2224 {
2226 tmpVec[i], tmpVec[other], rotbnd, dir,
2227 sign * angle, tol);
2228 }
2229 else
2230 {
2232 tmpVec[i], tmpVec[other], rotbnd, dir,
2233 sign * angle, tol);
2234 }
2235
2236 if (nFaceVerts == 3)
2237 {
2238 ASSERTL0(
2241 "Unsupported face orientation for face " +
2242 std::to_string(ids[i]));
2243 }
2244
2245 // Look up vertices for this face.
2246 vector<int> per1 = vertMap[ids[i]];
2247 vector<int> per2 = vertMap[ids[other]];
2248
2249 // tmpMap will hold the pairs of vertices which are
2250 // periodic.
2251 map<int, pair<int, bool>> tmpMap;
2252
2253 // Use vmap to determine which vertices connect given
2254 // the orientation o.
2255 for (j = 0; j < nFaceVerts; ++j)
2256 {
2257 int v = vmap[nFaceVerts][o][j];
2258 tmpMap[per1[j]] =
2259 make_pair(per2[v], locVerts.count(per2[v]) > 0);
2260 }
2261
2262 // Now loop over tmpMap to associate periodic vertices.
2263 for (auto &mIt : tmpMap)
2264 {
2265 PeriodicEntity ent2(mIt.second.first,
2267 mIt.second.second);
2268
2269 // See if this vertex has been recorded already.
2270 auto perIt = periodicVerts.find(mIt.first);
2271
2272 if (perIt == periodicVerts.end())
2273 {
2274 // Vertex is new - just record this entity as
2275 // usual.
2276 periodicVerts[mIt.first].push_back(ent2);
2277 perIt = periodicVerts.find(mIt.first);
2278 }
2279 else
2280 {
2281 // Vertex is known - loop over the vertices
2282 // inside the record and potentially add vertex
2283 // mIt.second to the list.
2284 for (k = 0; k < perIt->second.size(); ++k)
2285 {
2286 if (perIt->second[k].id == mIt.second.first)
2287 {
2288 break;
2289 }
2290 }
2291
2292 if (k == perIt->second.size())
2293 {
2294 perIt->second.push_back(ent2);
2295 }
2296 }
2297 }
2298 }
2299
2300 // Determine periodic edges. Logic is the same as above,
2301 // and perhaps should be condensed to avoid replication.
2302 for (i = 0; i < 2; ++i)
2303 {
2304 int other = (i + 1) % 2;
2305
2306 // angle is set up for i = 0 to i = 1
2307 sign = (i == 0) ? 1.0 : -1.0;
2308
2309 if (tmpVec[0].size() == 3)
2310 {
2312 tmpVec[i], tmpVec[other], rotbnd, dir,
2313 sign * angle, tol);
2314 }
2315 else
2316 {
2318 tmpVec[i], tmpVec[other], rotbnd, dir,
2319 sign * angle, tol);
2320 }
2321
2322 vector<int> per1 = edgeMap[ids[i]];
2323 vector<int> per2 = edgeMap[ids[other]];
2324
2325 map<int, pair<int, bool>> tmpMap;
2326
2327 for (j = 0; j < nFaceVerts; ++j)
2328 {
2329 int e = emap[nFaceVerts][o][j];
2330 tmpMap[per1[j]] =
2331 make_pair(per2[e], locEdges.count(per2[e]) > 0);
2332 }
2333
2334 for (auto &mIt : tmpMap)
2335 {
2336 // Note we assume orientation of edges is forwards -
2337 // this may be reversed later.
2338 PeriodicEntity ent2(mIt.second.first,
2340 mIt.second.second);
2341 auto perIt = periodicEdges.find(mIt.first);
2342
2343 if (perIt == periodicEdges.end())
2344 {
2345 periodicEdges[mIt.first].push_back(ent2);
2346 perIt = periodicEdges.find(mIt.first);
2347 }
2348 else
2349 {
2350 for (k = 0; k < perIt->second.size(); ++k)
2351 {
2352 if (perIt->second[k].id == mIt.second.first)
2353 {
2354 break;
2355 }
2356 }
2357
2358 if (k == perIt->second.size())
2359 {
2360 perIt->second.push_back(ent2);
2361 }
2362 }
2363 }
2364 }
2365 }
2366 }
2367
2368 // Make sure that the nubmer of face pairs and the
2369 // face Id to composite Id map match in size
2370 ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
2371 "At this point the size of allCompPairs "
2372 "should have been the same as fIdToCompId");
2373
2374 // also will need an edge id to composite id at
2375 // end of routine
2376 map<int, int> eIdToCompId;
2377
2378 // Search for periodic vertices and edges which are not
2379 // in a periodic composite but lie in this process. First,
2380 // loop over all information we have from other
2381 // processors.
2382 for (cnt = i = 0; i < totFaces; ++i)
2383 {
2384 bool rotbnd = false;
2385 int dir = 0;
2386 NekDouble angle = 0.0;
2387 NekDouble tol = 1e-8;
2388
2389 int faceId = faceIds[i];
2390
2391 ASSERTL0(allCompPairs.count(faceId) > 0,
2392 "Unable to find matching periodic face. faceId = " +
2393 std::to_string(faceId));
2394
2395 int perFaceId = allCompPairs[faceId];
2396
2397 // check to see if periodic boundary is rotated
2398 ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
2399 "Face " + std::to_string(faceId) +
2400 " not found in fIdtoCompId map");
2401 if (rotComp.count(fIdToCompId[faceId]))
2402 {
2403 rotbnd = true;
2404 dir = rotComp[fIdToCompId[faceId]].m_dir;
2405 angle = rotComp[fIdToCompId[faceId]].m_angle;
2406 tol = rotComp[fIdToCompId[faceId]].m_tol;
2407 }
2408
2409 for (j = 0; j < faceVerts[i]; ++j, ++cnt)
2410 {
2411 int vId = vertIds[cnt];
2412
2413 auto perId = periodicVerts.find(vId);
2414
2415 if (perId == periodicVerts.end())
2416 {
2417
2418 // This vertex is not included in the
2419 // map. Figure out which vertex it is supposed
2420 // to be periodic with. perFaceId is the face
2421 // ID which is periodic with faceId. The logic
2422 // is much the same as the loop above.
2424 coordMap[faceId], coordMap[perFaceId]};
2425
2426 int nFaceVerts = tmpVec[0].size();
2428 nFaceVerts == 3
2430 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2431 tol)
2433 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2434 tol);
2435
2436 // Use vmap to determine which vertex of the other face
2437 // should be periodic with this one.
2438 int perVertexId =
2439 vertMap[perFaceId][vmap[nFaceVerts][o][j]];
2440
2441 PeriodicEntity ent(perVertexId,
2443 locVerts.count(perVertexId) > 0);
2444
2445 periodicVerts[vId].push_back(ent);
2446 }
2447
2448 int eId = edgeIds[cnt];
2449
2450 perId = periodicEdges.find(eId);
2451
2452 // this map is required at very end to determine rotation of
2453 // edges.
2454 if (rotbnd)
2455 {
2456 eIdToCompId[eId] = fIdToCompId[faceId];
2457 }
2458
2459 if (perId == periodicEdges.end())
2460 {
2461 // This edge is not included in the map. Figure
2462 // out which edge it is supposed to be periodic
2463 // with. perFaceId is the face ID which is
2464 // periodic with faceId. The logic is much the
2465 // same as the loop above.
2467 coordMap[faceId], coordMap[perFaceId]};
2468
2469 int nFaceEdges = tmpVec[0].size();
2471 nFaceEdges == 3
2473 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2474 tol)
2476 tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2477 tol);
2478
2479 // Use emap to determine which edge of the other
2480 // face should be periodic with this one.
2481 int perEdgeId =
2482 edgeMap[perFaceId][emap[nFaceEdges][o][j]];
2483
2485 locEdges.count(perEdgeId) > 0);
2486
2487 periodicEdges[eId].push_back(ent);
2488
2489 // this map is required at very end to
2490 // determine rotation of edges.
2491 if (rotbnd)
2492 {
2493 eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
2494 }
2495 }
2496 }
2497 }
2498
2499 // Finally, we must loop over the periodicVerts and periodicEdges
2500 // map to complete connectivity information.
2501 for (auto &perIt : periodicVerts)
2502 {
2503 // For each vertex that is periodic with this one...
2504 for (i = 0; i < perIt.second.size(); ++i)
2505 {
2506 // Find the vertex in the periodicVerts map...
2507 auto perIt2 = periodicVerts.find(perIt.second[i].id);
2508 ASSERTL0(perIt2 != periodicVerts.end(),
2509 "Couldn't find periodic vertex.");
2510
2511 // Now search through this vertex's list and make sure that
2512 // we have a record of any vertices which aren't in the
2513 // original list.
2514 for (j = 0; j < perIt2->second.size(); ++j)
2515 {
2516 if (perIt2->second[j].id == perIt.first)
2517 {
2518 continue;
2519 }
2520
2521 for (k = 0; k < perIt.second.size(); ++k)
2522 {
2523 if (perIt2->second[j].id == perIt.second[k].id)
2524 {
2525 break;
2526 }
2527 }
2528
2529 if (k == perIt.second.size())
2530 {
2531 perIt.second.push_back(perIt2->second[j]);
2532 }
2533 }
2534 }
2535 }
2536
2537 for (auto &perIt : periodicEdges)
2538 {
2539 for (i = 0; i < perIt.second.size(); ++i)
2540 {
2541 auto perIt2 = periodicEdges.find(perIt.second[i].id);
2542 ASSERTL0(perIt2 != periodicEdges.end(),
2543 "Couldn't find periodic edge.");
2544
2545 for (j = 0; j < perIt2->second.size(); ++j)
2546 {
2547 if (perIt2->second[j].id == perIt.first)
2548 {
2549 continue;
2550 }
2551
2552 for (k = 0; k < perIt.second.size(); ++k)
2553 {
2554 if (perIt2->second[j].id == perIt.second[k].id)
2555 {
2556 break;
2557 }
2558 }
2559
2560 if (k == perIt.second.size())
2561 {
2562 perIt.second.push_back(perIt2->second[j]);
2563 }
2564 }
2565 }
2566 }
2567
2568 // Loop over periodic edges to determine relative edge orientations.
2569 for (auto &perIt : periodicEdges)
2570 {
2571 bool rotbnd = false;
2572 int dir = 0;
2573 NekDouble angle = 0.0;
2574 NekDouble tol = 1e-8;
2575
2576 // Find edge coordinates
2577 auto eIt = eIdMap.find(perIt.first);
2578 SpatialDomains::PointGeom v[2] = {*vCoMap[eIt->second.first],
2579 *vCoMap[eIt->second.second]};
2580
2581 // check to see if perioid boundary is rotated
2582 if (rotComp.count(eIdToCompId[perIt.first]))
2583 {
2584 rotbnd = true;
2585 dir = rotComp[eIdToCompId[perIt.first]].m_dir;
2586 angle = rotComp[eIdToCompId[perIt.first]].m_angle;
2587 tol = rotComp[eIdToCompId[perIt.first]].m_tol;
2588 }
2589
2590 // Loop over each edge, and construct a vector that takes us
2591 // from one vertex to another. Use this to figure out which
2592 // vertex maps to which.
2593 for (i = 0; i < perIt.second.size(); ++i)
2594 {
2595 eIt = eIdMap.find(perIt.second[i].id);
2596
2598 *vCoMap[eIt->second.first],
2599 *vCoMap[eIt->second.second]};
2600
2601 int vMap[2] = {-1, -1};
2602 if (rotbnd)
2603 {
2604
2606
2607 r.Rotate(v[0], dir, angle);
2608
2609 if (r.dist(w[0]) < tol)
2610 {
2611 vMap[0] = 0;
2612 }
2613 else
2614 {
2615 r.Rotate(v[1], dir, angle);
2616 if (r.dist(w[0]) < tol)
2617 {
2618 vMap[0] = 1;
2619 }
2620 else
2621 {
2623 "Unable to align rotationally "
2624 "periodic edge vertex");
2625 }
2626 }
2627 }
2628 else // translation test
2629 {
2630 NekDouble cx =
2631 0.5 * (w[0](0) - v[0](0) + w[1](0) - v[1](0));
2632 NekDouble cy =
2633 0.5 * (w[0](1) - v[0](1) + w[1](1) - v[1](1));
2634 NekDouble cz =
2635 0.5 * (w[0](2) - v[0](2) + w[1](2) - v[1](2));
2636
2637 for (j = 0; j < 2; ++j)
2638 {
2639 NekDouble x = v[j](0);
2640 NekDouble y = v[j](1);
2641 NekDouble z = v[j](2);
2642 for (k = 0; k < 2; ++k)
2643 {
2644 NekDouble x1 = w[k](0) - cx;
2645 NekDouble y1 = w[k](1) - cy;
2646 NekDouble z1 = w[k](2) - cz;
2647
2648 if (sqrt((x1 - x) * (x1 - x) +
2649 (y1 - y) * (y1 - y) +
2650 (z1 - z) * (z1 - z)) < 1e-8)
2651 {
2652 vMap[k] = j;
2653 break;
2654 }
2655 }
2656 }
2657
2658 // Sanity check the map.
2659 ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
2660 "Unable to align periodic edge vertex.");
2661 ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
2662 (vMap[1] == 0 || vMap[1] == 1) &&
2663 (vMap[0] != vMap[1]),
2664 "Unable to align periodic edge vertex.");
2665 }
2666
2667 // If 0 -> 0 then edges are aligned already; otherwise
2668 // reverse the orientation.
2669 if (vMap[0] != 0)
2670 {
2671 perIt.second[i].orient = StdRegions::eBackwards;
2672 }
2673 }
2674 }
2675
2676 // Do one final loop over periodic vertices/edges to remove
2677 // non-local vertices/edges from map.
2678 for (auto &perIt : periodicVerts)
2679 {
2680 if (locVerts.count(perIt.first) > 0)
2681 {
2682 m_periodicVerts.insert(perIt);
2683 }
2684 }
2685
2686 for (auto &perIt : periodicEdges)
2687 {
2688 if (locEdges.count(perIt.first) > 0)
2689 {
2690 m_periodicEdges.insert(perIt);
2691 }
2692 }
2693 }
2694 break;
2695 default:
2696 ASSERTL1(false, "Not setup for this expansion");
2697 break;
2698 }
2699}
2700
2701/**
2702 *
2703 */
2705 const GlobalLinSysKey &mkey)
2706{
2708 "Routine currently only tested for HybridDGHelmholtz");
2709
2711 "Full matrix global systems are not supported for HDG "
2712 "expansions");
2713
2714 ASSERTL1(mkey.GetGlobalSysSolnType() == m_traceMap->GetGlobalSysSolnType(),
2715 "The local to global map is not set up for the requested "
2716 "solution type");
2717
2718 GlobalLinSysSharedPtr glo_matrix;
2719 auto matrixIter = m_globalBndMat->find(mkey);
2720
2721 if (matrixIter == m_globalBndMat->end())
2722 {
2723 glo_matrix = GenGlobalBndLinSys(mkey, m_traceMap);
2724 (*m_globalBndMat)[mkey] = glo_matrix;
2725 }
2726 else
2727 {
2728 glo_matrix = matrixIter->second;
2729 }
2730
2731 return glo_matrix;
2732}
2733
2734/**
2735 * For each boundary region, checks that the types and number of
2736 * boundary expansions in that region match.
2737 * @param In Field to compare with.
2738 * @return True if boundary conditions match.
2739 */
2741{
2742 int i;
2743 bool returnval = true;
2744
2745 for (i = 0; i < m_bndConditions.size(); ++i)
2746 {
2747 // check to see if boundary condition type is the same
2748 // and there are the same number of boundary
2749 // conditions in the boundary definition.
2750 if ((m_bndConditions[i]->GetBoundaryConditionType() !=
2751 In.m_bndConditions[i]->GetBoundaryConditionType()) ||
2752 (m_bndCondExpansions[i]->GetExpSize() !=
2753 In.m_bndCondExpansions[i]->GetExpSize()))
2754 {
2755 returnval = false;
2756 break;
2757 }
2758 }
2759
2760 // Compare with all other processes. Return true only if all
2761 // processes report having the same boundary conditions.
2762 int vSame = (returnval ? 1 : 0);
2763 m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
2764
2765 return (vSame == 1);
2766}
2767
2769{
2770 if (m_negatedFluxNormal.size() == 0)
2771 {
2773 &elmtToTrace = m_traceMap->GetElmtToTrace();
2774
2775 m_negatedFluxNormal.resize(2 * GetExpSize());
2776
2777 for (int i = 0; i < GetExpSize(); ++i)
2778 {
2779
2780 for (int v = 0; v < 2; ++v)
2781 {
2783 elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
2784
2785 if (vertExp->GetLeftAdjacentElementExp()
2786 ->GetGeom()
2787 ->GetGlobalID() !=
2788 (*m_exp)[i]->GetGeom()->GetGlobalID())
2789 {
2790 m_negatedFluxNormal[2 * i + v] = true;
2791 }
2792 else
2793 {
2794 m_negatedFluxNormal[2 * i + v] = false;
2795 }
2796 }
2797 }
2798 }
2799
2800 return m_negatedFluxNormal;
2801}
2802
2804{
2806 {
2807 SetUpDG();
2808 }
2809
2810 return m_trace;
2811}
2812
2814{
2815 return m_traceMap;
2816}
2817
2819{
2820 return m_interfaceMap;
2821}
2822
2824 void) const
2825{
2826 return m_locTraceToTraceMap;
2827}
2828
2829// Return the internal vector which identifieds if trace
2830// is left adjacent definiing which trace the normal
2831// points otwards from
2833{
2834 return m_leftAdjacentTraces;
2835}
2836
2839{
2840 return m_bndCondExpansions;
2841}
2842
2845{
2846 return m_bndConditions;
2847}
2848
2850{
2851 return m_bndCondExpansions[i];
2852}
2853
2856{
2857 return m_bndConditions;
2858}
2859
2862{
2863 for (int n = 0; n < m_periodicFwdCopy.size(); ++n)
2864 {
2865 Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
2866 }
2867}
2868
2871{
2872 v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
2873}
2874
2875/**
2876 * @brief Obtain a copy of the periodic edges and vertices for this
2877 * field.
2878 */
2880 PeriodicMap &periodicEdges,
2881 PeriodicMap &periodicFaces)
2882{
2883 periodicVerts = m_periodicVerts;
2884 periodicEdges = m_periodicEdges;
2885 periodicFaces = m_periodicFaces;
2886}
2887/**
2888 * \brief This method extracts the "forward" and "backward" trace
2889 * data from the array \a field and puts the data into output
2890 * vectors \a Fwd and \a Bwd.
2891 *
2892 * We first define the convention which defines "forwards" and
2893 * "backwards". First an association is made between the vertex/edge/face of
2894 * each element and its corresponding vertex/edge/face in the trace space
2895 * using the mapping #m_traceMap. The element can either be
2896 * left-adjacent or right-adjacent to this trace face (see
2897 * Expansion2D::GetLeftAdjacentElementExp). Boundary faces are
2898 * always left-adjacent since left-adjacency is populated first.
2899 *
2900 * If the element is left-adjacent we extract the trace data
2901 * from \a field into the forward trace space \a Fwd; otherwise,
2902 * we place it in the backwards trace space \a Bwd. In this way,
2903 * we form a unique set of trace normals since these are always
2904 * extracted from left-adjacent elements.
2905 *
2906 * \param field is a NekDouble array which contains the fielddata
2907 * from which we wish to extract the backward and forward
2908 * orientated trace/face arrays.
2909 *
2910 * \return Updates a NekDouble array \a Fwd and \a Bwd
2911 */
2916 const Array<OneD, const ExpListSharedPtr> &BndCondExp)
2917{
2918 ASSERTL1((*m_exp)[0]->GetBasis(0)->GetBasisType() !=
2920 "Routine not set up for Gauss Lagrange points distribution");
2921
2922 // blocked routine
2923 Array<OneD, NekDouble> tracevals(m_locTraceToTraceMap->GetNLocTracePts());
2924
2925 m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2926 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2927
2928 Array<OneD, NekDouble> invals =
2929 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2930 m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2931
2933
2934 FillBwdWithBoundCond(Fwd, Bwd, BndCond, BndCondExp, false);
2935
2936 // Do parallel exchange for forwards/backwards spaces.
2937 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2938
2939 // Do exchange of interface traces (local and parallel)
2940 // We may have to split this out into separate local and
2941 // parallel for IP method???
2942 m_interfaceMap->ExchangeTrace(Fwd, Bwd);
2943}
2944
2948 bool FillBnd, // these should be template params
2949 bool PutFwdInBwdOnBCs, bool DoExchange)
2950{
2951 // Is this zeroing necessary?
2952 // Zero forward/backward vectors.
2953 Vmath::Zero(Fwd.size(), Fwd, 1);
2954 Vmath::Zero(Bwd.size(), Bwd, 1);
2955
2956 // Basis definition on each element
2957 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
2958 if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
2959 {
2960 // blocked routine
2961 Array<OneD, NekDouble> tracevals(
2962 m_locTraceToTraceMap->GetNLocTracePts());
2963
2964 m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2965 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2966
2967 Array<OneD, NekDouble> invals =
2968 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2969 m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2970 }
2971 else
2972 {
2973 // Loop over elements and collect forward expansion
2974 auto nexp = GetExpSize();
2977
2979 &elmtToTrace = m_traceMap->GetElmtToTrace();
2980
2981 for (int n = 0, cnt = 0; n < nexp; ++n)
2982 {
2983 exp = (*m_exp)[n];
2984 auto phys_offset = GetPhys_Offset(n);
2985
2986 for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2987 {
2988 auto offset =
2989 m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2990
2991 e_tmp =
2992 (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
2993
2994 exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
2995 e_tmp);
2996 }
2997 }
2998 }
2999
3001
3002 if (FillBnd)
3003 {
3004 v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
3005 }
3006
3007 if (DoExchange)
3008 {
3009 // Do parallel exchange for forwards/backwards spaces.
3010 m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
3011
3012 // Do exchange of interface traces (local and parallel)
3013 // We may have to split this out into separate local and
3014 // parallel for IP method???
3015 m_interfaceMap->ExchangeTrace(Fwd, Bwd);
3016 }
3017}
3018
3021 bool PutFwdInBwdOnBCs)
3022{
3024 PutFwdInBwdOnBCs);
3025}
3026
3027/** function allowing different BCs to be passed to routine
3028 */
3032 &bndConditions,
3033 const Array<OneD, const ExpListSharedPtr> &bndCondExpansions,
3034 bool PutFwdInBwdOnBCs)
3035{
3036
3037 // Fill boundary conditions into missing elements
3038 if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
3039 {
3040 // Fill boundary conditions into missing elements
3041 for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3042 {
3043 if (bndConditions[n]->GetBoundaryConditionType() ==
3045 {
3046 auto ne = bndCondExpansions[n]->GetExpSize();
3047 for (int e = 0; e < ne; ++e)
3048 {
3049 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3050 auto id2 = m_trace->GetPhys_Offset(
3051 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3052 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3053 }
3054
3055 cnt += ne;
3056 }
3057 else if (bndConditions[n]->GetBoundaryConditionType() ==
3059 bndConditions[n]->GetBoundaryConditionType() ==
3061 {
3062 auto ne = bndCondExpansions[n]->GetExpSize();
3063 for (int e = 0; e < ne; ++e)
3064 {
3065 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3066 auto id2 = m_trace->GetPhys_Offset(
3067 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3068 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3069 }
3070 cnt += ne;
3071 }
3072 else if (bndConditions[n]->GetBoundaryConditionType() !=
3074 {
3075 ASSERTL0(false, "Method not set up for this "
3076 "boundary condition.");
3077 }
3078 }
3079 }
3080 else
3081 {
3082 for (int n = 0, cnt = 0; n < bndCondExpansions.size(); ++n)
3083 {
3084 if (bndConditions[n]->GetBoundaryConditionType() ==
3086 {
3087 auto ne = bndCondExpansions[n]->GetExpSize();
3088 for (int e = 0; e < ne; ++e)
3089 {
3090 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3091 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3092 auto id2 = m_trace->GetPhys_Offset(
3093 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3094
3095 Vmath::Vcopy(npts, &(bndCondExpansions[n]->GetPhys())[id1],
3096 1, &Bwd[id2], 1);
3097 }
3098 cnt += ne;
3099 }
3100 else if (bndConditions[n]->GetBoundaryConditionType() ==
3102 bndConditions[n]->GetBoundaryConditionType() ==
3104 {
3105 auto ne = m_bndCondExpansions[n]->GetExpSize();
3106 for (int e = 0; e < ne; ++e)
3107 {
3108 auto npts = bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3109 auto id1 = bndCondExpansions[n]->GetPhys_Offset(e);
3110
3111 ASSERTL0((bndCondExpansions[n]->GetPhys())[id1] == 0.0,
3112 "Method not set up for non-zero "
3113 "Neumann boundary condition");
3114 auto id2 = m_trace->GetPhys_Offset(
3115 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3116
3117 Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
3118 }
3119
3120 cnt += ne;
3121 }
3122 else if (bndConditions[n]->GetBoundaryConditionType() ==
3124 {
3125 // Do nothing
3126 }
3127 else if (bndConditions[n]->GetBoundaryConditionType() !=
3129 {
3131 "Method not set up for this boundary "
3132 "condition.");
3133 }
3134 }
3135 }
3136}
3137
3139{
3140 return m_bndCondBndWeight;
3141}
3142
3143void DisContField::v_SetBndCondBwdWeight(const int index, const NekDouble value)
3144{
3145 m_bndCondBndWeight[index] = value;
3146}
3147
3151{
3152 // Basis definition on each element
3153 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3154 if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
3155 {
3156 Array<OneD, NekDouble> tracevals(
3157 m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
3158
3159 Array<OneD, NekDouble> invals =
3160 tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
3161
3162 m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
3163
3164 m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
3165
3166 m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
3167 }
3168 else
3169 {
3171 "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3172 }
3173}
3174
3176{
3177 ASSERTL1(m_physState == true, "local physical space is not true ");
3178 v_ExtractTracePhys(m_phys, outarray);
3179}
3180/**
3181 * @brief This method extracts the trace (verts in 1D) from
3182 * the field @a inarray and puts the values in @a outarray.
3183 *
3184 * It assumes the field is C0 continuous so that it can
3185 * overwrite the edge data when visited by the two adjacent
3186 * elements.
3187 *
3188 * @param inarray An array containing the 1D data from which we wish
3189 * to extract the edge data.
3190 * @param outarray The resulting edge information.
3191 *
3192 * This will not work for non-boundary expansions
3193 */
3195 const Array<OneD, const NekDouble> &inarray,
3196 Array<OneD, NekDouble> &outarray)
3197{
3198 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3199 if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3200 {
3201 Vmath::Zero(outarray.size(), outarray, 1);
3202 Array<OneD, NekDouble> tracevals(
3203 m_locTraceToTraceMap->GetNFwdLocTracePts());
3204 m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
3205 m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
3206 m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3207 }
3208 else
3209 {
3210 // Loop over elemente and collect forward expansion
3211 int nexp = GetExpSize();
3212 int n, p, offset, phys_offset;
3215 &elmtToTrace = m_traceMap->GetElmtToTrace();
3216 ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
3217 "input array is of insufficient length");
3218 for (n = 0; n < nexp; ++n)
3219 {
3220 phys_offset = GetPhys_Offset(n);
3221 for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
3222 {
3223 offset =
3224 m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
3225 (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
3226 inarray + phys_offset,
3227 t_tmp = outarray + offset);
3228 }
3229 }
3230 }
3231}
3232/**
3233 * @brief Add trace contributions into elemental coefficient spaces.
3234 *
3235 * Given some quantity \f$ \vec{Fn} \f$, which conatins this
3236 * routine calculates the integral
3237 *
3238 * \f[
3239 * \int_{\Omega^e} \vec{Fn}, \mathrm{d}S
3240 * \f]
3241 *
3242 * and adds this to the coefficient space provided by outarray.
3243 *
3244 * @see Expansion3D::AddFaceNormBoundaryInt
3245 *
3246 * @param Fn The trace quantities.
3247 * @param outarray Resulting 3D coefficient space.
3248 *
3249 */
3251 Array<OneD, NekDouble> &outarray)
3252{
3253 if (m_expType == e1D)
3254 {
3255 int n, offset, t_offset;
3257 &elmtToTrace = m_traceMap->GetElmtToTrace();
3258 vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
3259 for (n = 0; n < GetExpSize(); ++n)
3260 {
3261 // Number of coefficients on each element
3262 int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3263 offset = GetCoeff_Offset(n);
3264 // Implementation for every points except Gauss points
3265 if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
3267 {
3268 t_offset =
3269 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3270 if (negatedFluxNormal[2 * n])
3271 {
3272 outarray[offset] -= Fn[t_offset];
3273 }
3274 else
3275 {
3276 outarray[offset] += Fn[t_offset];
3277 }
3278 t_offset =
3279 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3280 if (negatedFluxNormal[2 * n + 1])
3281 {
3282 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3283 Fn[t_offset];
3284 }
3285 else
3286 {
3287 outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3288 Fn[t_offset];
3289 }
3290 }
3291 else
3292 {
3293 int j;
3294 static DNekMatSharedPtr m_Ixm, m_Ixp;
3295 static int sav_ncoeffs = 0;
3296 if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3297 {
3299 const LibUtilities::PointsKey BS_p(
3301 const LibUtilities::BasisKey BS_k(
3302 LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
3303 BASE = LibUtilities::BasisManager()[BS_k];
3304 Array<OneD, NekDouble> coords(1, 0.0);
3305 coords[0] = -1.0;
3306 m_Ixm = BASE->GetI(coords);
3307 coords[0] = 1.0;
3308 m_Ixp = BASE->GetI(coords);
3309 sav_ncoeffs = e_ncoeffs;
3310 }
3311 t_offset =
3312 GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3313 if (negatedFluxNormal[2 * n])
3314 {
3315 for (j = 0; j < e_ncoeffs; j++)
3316 {
3317 outarray[offset + j] -=
3318 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3319 }
3320 }
3321 else
3322 {
3323 for (j = 0; j < e_ncoeffs; j++)
3324 {
3325 outarray[offset + j] +=
3326 (m_Ixm->GetPtr())[j] * Fn[t_offset];
3327 }
3328 }
3329 t_offset =
3330 GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3331 if (negatedFluxNormal[2 * n + 1])
3332 {
3333 for (j = 0; j < e_ncoeffs; j++)
3334 {
3335 outarray[offset + j] -=
3336 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3337 }
3338 }
3339 else
3340 {
3341 for (j = 0; j < e_ncoeffs; j++)
3342 {
3343 outarray[offset + j] +=
3344 (m_Ixp->GetPtr())[j] * Fn[t_offset];
3345 }
3346 }
3347 }
3348 }
3349 }
3350 else // other dimensions
3351 {
3352 // Basis definition on each element
3353 LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3354 if ((m_expType != e1D) &&
3355 (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3356 {
3357 Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
3358 m_trace->IProductWRTBase(Fn, Fcoeffs);
3359 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
3360 outarray);
3361 }
3362 else
3363 {
3364 int e, n, offset, t_offset;
3365 Array<OneD, NekDouble> e_outarray;
3367 &elmtToTrace = m_traceMap->GetElmtToTrace();
3368 if (m_expType == e2D)
3369 {
3370 for (n = 0; n < GetExpSize(); ++n)
3371 {
3372 offset = GetCoeff_Offset(n);
3373 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3374 {
3375 t_offset = GetTrace()->GetPhys_Offset(
3376 elmtToTrace[n][e]->GetElmtId());
3377 (*m_exp)[n]->AddEdgeNormBoundaryInt(
3378 e, elmtToTrace[n][e], Fn + t_offset,
3379 e_outarray = outarray + offset);
3380 }
3381 }
3382 }
3383 else if (m_expType == e3D)
3384 {
3385 for (n = 0; n < GetExpSize(); ++n)
3386 {
3387 offset = GetCoeff_Offset(n);
3388 for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3389 {
3390 t_offset = GetTrace()->GetPhys_Offset(
3391 elmtToTrace[n][e]->GetElmtId());
3392 (*m_exp)[n]->AddFaceNormBoundaryInt(
3393 e, elmtToTrace[n][e], Fn + t_offset,
3394 e_outarray = outarray + offset);
3395 }
3396 }
3397 }
3398 }
3399 }
3400}
3401/**
3402 * @brief Add trace contributions into elemental coefficient spaces.
3403 *
3404 * Given some quantity \f$ \vec{q} \f$, calculate the elemental integral
3405 *
3406 * \f[
3407 * \int_{\Omega^e} \vec{q}, \mathrm{d}S
3408 * \f]
3409 *
3410 * and adds this to the coefficient space provided by
3411 * outarray. The value of q is determined from the routine
3412 * IsLeftAdjacentTrace() which if true we use Fwd else we use
3413 * Bwd
3414 *
3415 * @param Fwd The trace quantities associated with left (fwd)
3416 * adjancent elmt.
3417 * @param Bwd The trace quantities associated with right (bwd)
3418 * adjacent elet.
3419 * @param outarray Resulting coefficient space.
3420 */
3424{
3425 ASSERTL0(m_expType != e1D, "This method is not setup or "
3426 "tested for 1D expansion");
3427 Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
3428 m_trace->IProductWRTBase(Fwd, Coeffs);
3429 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, Coeffs, outarray);
3430 m_trace->IProductWRTBase(Bwd, Coeffs);
3431 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, Coeffs, outarray);
3432}
3433
3435 const Array<OneD, const NekDouble> &inarray,
3437 const StdRegions::VarCoeffMap &varcoeff,
3438 [[maybe_unused]] const MultiRegions::VarFactorsMap &varfactors,
3439 [[maybe_unused]] const Array<OneD, const NekDouble> &dirForcing,
3440 const bool PhysSpaceForcing)
3441{
3442 int i, n, cnt, nbndry;
3443 int nexp = GetExpSize();
3445 DNekVec F(m_ncoeffs, f, eWrapper);
3446 Array<OneD, NekDouble> e_f, e_l;
3447 //----------------------------------
3448 // Setup RHS Inner product if required
3449 //----------------------------------
3450 if (PhysSpaceForcing)
3451 {
3452 IProductWRTBase(inarray, f);
3453 Vmath::Neg(m_ncoeffs, f, 1);
3454 }
3455 else
3456 {
3457 Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
3458 }
3459 //----------------------------------
3460 // Solve continuous Boundary System
3461 //----------------------------------
3462 int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
3463 int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
3464 int e_ncoeffs;
3467 const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
3468 // Retrieve number of local trace space coefficients N_{\lambda},
3469 // and set up local elemental trace solution \lambda^e.
3470 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
3471 Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
3472 Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
3473 DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
3474 //----------------------------------
3475 // Evaluate Trace Forcing
3476 // Kirby et al, 2010, P23, Step 5.
3477 //----------------------------------
3478 // Determing <u_lam,f> terms using HDGLamToU matrix
3479 for (cnt = n = 0; n < nexp; ++n)
3480 {
3481 nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3482 e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3483 e_f = f + m_coeff_offset[n];
3484 e_l = bndrhs + cnt;
3485 // use outarray as tmp space
3486 DNekVec Floc(nbndry, e_l, eWrapper);
3487 DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
3488 Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3489 cnt += nbndry;
3490 }
3491 Array<OneD, const int> bndCondMap =
3492 m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3493 Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
3494 // Copy Dirichlet boundary conditions and weak forcing
3495 // into trace space
3496 int locid;
3497 cnt = 0;
3498 for (i = 0; i < m_bndCondExpansions.size(); ++i)
3499 {
3500 const Array<OneD, const NekDouble> bndcoeffs =
3501 m_bndCondExpansions[i]->GetCoeffs();
3502
3503 if (m_bndConditions[i]->GetBoundaryConditionType() ==
3505 {
3506 for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3507 {
3508 locid = bndCondMap[cnt + j];
3509 loclambda[locid] = Sign[locid] * bndcoeffs[j];
3510 }
3511 }
3512 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3514 m_bndConditions[i]->GetBoundaryConditionType() ==
3516 {
3517 // Add weak boundary condition to trace forcing
3518 for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3519 {
3520 locid = bndCondMap[cnt + j];
3521 bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3522 }
3523 }
3524 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3526 {
3527 // skip increment of cnt if ePeriodic
3528 // because bndCondMap does not include ePeriodic
3529 continue;
3530 }
3531 cnt += (m_bndCondExpansions[i])->GetNcoeffs();
3532 }
3533
3534 //----------------------------------
3535 // Solve trace problem: \Lambda = K^{-1} F
3536 // K is the HybridDGHelmBndLam matrix.
3537 //----------------------------------
3538 if (GloBndDofs - NumDirBCs > 0)
3539 {
3541 factors, varcoeff);
3543 LinSys->Solve(bndrhs, loclambda, m_traceMap);
3544
3545 // For consistency with previous version put global
3546 // solution into m_trace->m_coeffs
3547 m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
3548 }
3549
3550 //----------------------------------
3551 // Internal element solves
3552 //----------------------------------
3555
3556 const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
3557 DNekVec out(m_ncoeffs, outarray, eWrapper);
3558 Vmath::Zero(m_ncoeffs, outarray, 1);
3559
3560 // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
3561 out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3562
3563 // Return empty GlobalLinSysKey
3564 return NullGlobalLinSysKey;
3565}
3566
3567/* \brief This function evaluates the boundary conditions at a certain
3568 * time-level.
3569 *
3570 * Based on the boundary condition \f$g(\boldsymbol{x},t)\f$ evaluated
3571 * at a given time-level \a t, this function transforms the boundary
3572 * conditions onto the coefficients of the (one-dimensional) boundary
3573 * expansion. Depending on the type of boundary conditions, these
3574 * expansion coefficients are calculated in different ways:
3575 * - <b>Dirichlet boundary conditions</b><BR>
3576 * In order to ensure global \f$C^0\f$ continuity of the spectral/hp
3577 * approximation, the Dirichlet boundary conditions are projected onto
3578 * the boundary expansion by means of a modified \f$C^0\f$ continuous
3579 * Galerkin projection. This projection can be viewed as a collocation
3580 * projection at the vertices, followed by an \f$L^2\f$ projection on
3581 * the interior modes of the edges. The resulting coefficients
3582 * \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ will be stored for the
3583 * boundary expansion.
3584 * - <b>Neumann boundary conditions</b>
3585 * In the discrete Galerkin formulation of the problem to be solved,
3586 * the Neumann boundary conditions appear as the set of surface
3587 * integrals: \f[\boldsymbol{\hat{g}}=\int_{\Gamma}
3588 * \phi^e_n(\boldsymbol{x})g(\boldsymbol{x})d(\boldsymbol{x})\quad
3589 * \forall n \f]
3590 * As a result, it are the coefficients \f$\boldsymbol{\hat{g}}\f$
3591 * that will be stored in the boundary expansion
3592 *
3593 * @param time The time at which the boundary conditions
3594 * should be evaluated.
3595 * @param bndCondExpansions List of boundary conditions.
3596 * @param bndConditions Information about the boundary conditions.
3597 *
3598 * This will only be undertaken for time dependent
3599 * boundary conditions unless time == 0.0 which is the
3600 * case when the method is called from the constructor.
3601 */
3603 const std::string varName,
3604 const NekDouble x2_in,
3605 const NekDouble x3_in)
3606{
3607 int i;
3608 int npoints;
3609
3611
3612 for (i = 0; i < m_bndCondExpansions.size(); ++i)
3613 {
3614 if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
3615 {
3616 m_bndCondBndWeight[i] = 1.0;
3617 locExpList = m_bndCondExpansions[i];
3618
3619 npoints = locExpList->GetNpoints();
3620 Array<OneD, NekDouble> x0(npoints, 0.0);
3621 Array<OneD, NekDouble> x1(npoints, 0.0);
3622 Array<OneD, NekDouble> x2(npoints, 0.0);
3623
3624 locExpList->GetCoords(x0, x1, x2);
3625
3626 if (x2_in != NekConstants::kNekUnsetDouble &&
3628 {
3629 Vmath::Fill(npoints, x2_in, x1, 1);
3630 Vmath::Fill(npoints, x3_in, x2, 1);
3631 }
3632 else if (x2_in != NekConstants::kNekUnsetDouble)
3633 {
3634 Vmath::Fill(npoints, x2_in, x2, 1);
3635 }
3636
3637 // treat 1D expansions separately since we only
3638 // require an evaluation at a point rather than
3639 // any projections or inner products that are not
3640 // available in a PointExp
3641 if (m_expType == e1D)
3642 {
3643 if (m_bndConditions[i]->GetBoundaryConditionType() ==
3645 {
3646
3647 m_bndCondExpansions[i]->SetCoeff(
3648 0, (std::static_pointer_cast<
3650 m_bndConditions[i])
3651 ->m_dirichletCondition)
3652 .Evaluate(x0[0], x1[0], x2[0], time));
3653 m_bndCondExpansions[i]->SetPhys(
3654 0, m_bndCondExpansions[i]->GetCoeff(0));
3655 }
3656 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3658 {
3659 m_bndCondExpansions[i]->SetCoeff(
3660 0, (std::static_pointer_cast<
3662 m_bndConditions[i])
3663 ->m_neumannCondition)
3664 .Evaluate(x0[0], x1[0], x2[0], time));
3665 }
3666 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3668 {
3669 m_bndCondExpansions[i]->SetCoeff(
3670 0, (std::static_pointer_cast<
3672 m_bndConditions[i])
3673 ->m_robinFunction)
3674 .Evaluate(x0[0], x1[0], x2[0], time));
3675 }
3676 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3678 {
3679 continue;
3680 }
3681 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3683 {
3684 }
3685 else
3686 {
3688 "This type of BC not implemented yet");
3689 }
3690 }
3691 else // 2D and 3D versions
3692 {
3693 if (m_bndConditions[i]->GetBoundaryConditionType() ==
3695 {
3697 std::static_pointer_cast<
3699 m_bndConditions[i]);
3700
3701 Array<OneD, NekDouble> valuesFile(npoints, 1.0),
3702 valuesExp(npoints, 1.0);
3703
3704 string bcfilename = bcPtr->m_filename;
3705 string exprbcs = bcPtr->m_expr;
3706
3707 if (bcfilename != "")
3708 {
3709 locExpList->ExtractCoeffsFromFile(
3710 bcfilename, bcPtr->GetComm(), varName,
3711 locExpList->UpdateCoeffs());
3712 locExpList->BwdTrans(locExpList->GetCoeffs(),
3713 locExpList->UpdatePhys());
3714 valuesFile = locExpList->GetPhys();
3715 }
3716
3717 if (exprbcs != "")
3718 {
3719 LibUtilities::Equation condition =
3720 std::static_pointer_cast<
3722 m_bndConditions[i])
3723 ->m_dirichletCondition;
3724
3725 condition.Evaluate(x0, x1, x2, time, valuesExp);
3726 }
3727
3728 Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
3729 locExpList->UpdatePhys(), 1);
3730 locExpList->FwdTransBndConstrained(
3731 locExpList->GetPhys(), locExpList->UpdateCoeffs());
3732 }
3733 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3735 {
3737 std::static_pointer_cast<
3739 m_bndConditions[i]);
3740 string bcfilename = bcPtr->m_filename;
3741
3742 if (bcfilename != "")
3743 {
3744 locExpList->ExtractCoeffsFromFile(
3745 bcfilename, bcPtr->GetComm(), varName,
3746 locExpList->UpdateCoeffs());
3747 locExpList->BwdTrans(locExpList->GetCoeffs(),
3748 locExpList->UpdatePhys());
3749 }
3750 else
3751 {
3752 LibUtilities::Equation condition =
3753 std::static_pointer_cast<
3755 m_bndConditions[i])
3756 ->m_neumannCondition;
3757 condition.Evaluate(x0, x1, x2, time,
3758 locExpList->UpdatePhys());
3759 }
3760
3761 locExpList->IProductWRTBase(locExpList->GetPhys(),
3762 locExpList->UpdateCoeffs());
3763 }
3764 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3766 {
3768 std::static_pointer_cast<
3770 m_bndConditions[i]);
3771
3772 string bcfilename = bcPtr->m_filename;
3773
3774 if (bcfilename != "")
3775 {
3776 locExpList->ExtractCoeffsFromFile(
3777 bcfilename, bcPtr->GetComm(), varName,
3778 locExpList->UpdateCoeffs());
3779 locExpList->BwdTrans(locExpList->GetCoeffs(),
3780 locExpList->UpdatePhys());
3781 }
3782 else
3783 {
3784 LibUtilities::Equation condition =
3785 std::static_pointer_cast<
3787 m_bndConditions[i])
3788 ->m_robinFunction;
3789 condition.Evaluate(x0, x1, x2, time,
3790 locExpList->UpdatePhys());
3791 }
3792
3793 locExpList->IProductWRTBase(locExpList->GetPhys(),
3794 locExpList->UpdateCoeffs());
3795 }
3796 else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3798 {
3799 continue;
3800 }
3801 else
3802 {
3804 "This type of BC not implemented yet");
3805 }
3806 }
3807 }
3808 }
3809}
3810
3811/**
3812 * @brief Fill the weight with m_bndCondBndWeight.
3813 */
3815 Array<OneD, NekDouble> &weightjmp)
3816{
3817 int cnt;
3818 int npts;
3819 int e = 0;
3820
3821 // Fill boundary conditions into missing elements
3822 int id2 = 0;
3823
3824 for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
3825 {
3826
3827 if (m_bndConditions[n]->GetBoundaryConditionType() ==
3829 {
3830 for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3831 {
3832 npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3833 id2 = m_trace->GetPhys_Offset(
3834 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3835 Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3836 Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3837 }
3838
3839 cnt += e;
3840 }
3841 else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3843 m_bndConditions[n]->GetBoundaryConditionType() ==
3845 {
3846 for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3847 {
3848 npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3849 id2 = m_trace->GetPhys_Offset(
3850 m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3851 Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3852 Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3853 }
3854
3855 cnt += e;
3856 }
3857 else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3859 {
3861 "Method not set up for this boundary condition.");
3862 }
3863 }
3864}
3865
3866// Set up a list of element ids and trace ids that link to the
3867// boundary conditions
3869 Array<OneD, int> &TraceID)
3870{
3871
3872 if (m_BCtoElmMap.size() == 0)
3873 {
3874 switch (m_expType)
3875 {
3876 case e1D:
3877 {
3878 map<int, int> VertGID;
3879 int i, n, id;
3880 int bid, cnt, Vid;
3881 int nbcs = 0;
3882 // Determine number of boundary condition expansions.
3883 for (i = 0; i < m_bndConditions.size(); ++i)
3884 {
3885 nbcs += m_bndCondExpansions[i]->GetExpSize();
3886 }
3887
3888 // make sure arrays are of sufficient length
3889 m_BCtoElmMap = Array<OneD, int>(nbcs, -1);
3891
3892 // setup map of all global ids along boundary
3893 cnt = 0;
3894 for (n = 0; n < m_bndCondExpansions.size(); ++n)
3895 {
3896 for (i = 0; i < m_bndCondExpansions[n]->GetExpSize(); ++i)
3897 {
3898 Vid = m_bndCondExpansions[n]
3899 ->GetExp(i)
3900 ->GetGeom()
3901 ->GetVertex(0)
3902 ->GetVid();
3903 VertGID[Vid] = cnt++;
3904 }
3905 }
3906
3907 // Loop over elements and find verts that match;
3909 for (cnt = n = 0; n < GetExpSize(); ++n)
3910 {
3911 exp = (*m_exp)[n];
3912 for (i = 0; i < exp->GetNverts(); ++i)
3913 {
3914 id = exp->GetGeom()->GetVid(i);
3915
3916 if (VertGID.count(id) > 0)
3917 {
3918 bid = VertGID.find(id)->second;
3919 ASSERTL1(m_BCtoElmMap[bid] == -1,
3920 "Edge already set");
3921 m_BCtoElmMap[bid] = n;
3922 m_BCtoTraceMap[bid] = i;
3923 cnt++;
3924 }
3925 }
3926 }
3927 ASSERTL1(cnt == nbcs,
3928 "Failed to visit all boundary condtiions");
3929 }
3930 break;
3931 case e2D:
3932 {
3933 map<int, int> globalIdMap;
3934 int i, n;
3935 int cnt;
3936 int nbcs = 0;
3937
3938 // Populate global ID map (takes global geometry
3939 // ID to local expansion list ID).
3940 for (i = 0; i < GetExpSize(); ++i)
3941 {
3942 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3943 }
3944
3945 // Determine number of boundary condition expansions.
3946 for (i = 0; i < m_bndConditions.size(); ++i)
3947 {
3948 nbcs += m_bndCondExpansions[i]->GetExpSize();
3949 }
3950
3951 // Initialize arrays
3954
3956 cnt = 0;
3957 for (n = 0; n < m_bndCondExpansions.size(); ++n)
3958 {
3959 for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3960 ++i, ++cnt)
3961 {
3962 exp1d = m_bndCondExpansions[n]
3963 ->GetExp(i)
3965
3966 // Use edge to element map from MeshGraph.
3968 m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3969
3970 m_BCtoElmMap[cnt] =
3971 globalIdMap[(*tmp)[0].first->GetGlobalID()];
3972 m_BCtoTraceMap[cnt] = (*tmp)[0].second;
3973 }
3974 }
3975 }
3976 break;
3977 case e3D:
3978 {
3979 map<int, int> globalIdMap;
3980 int i, n;
3981 int cnt;
3982 int nbcs = 0;
3983
3984 // Populate global ID map (takes global geometry ID to local
3985 // expansion list ID).
3987 for (i = 0; i < GetExpSize(); ++i)
3988 {
3989 globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3990 }
3991
3992 // Determine number of boundary condition expansions.
3993 for (i = 0; i < m_bndConditions.size(); ++i)
3994 {
3995 nbcs += m_bndCondExpansions[i]->GetExpSize();
3996 }
3997
3998 // Initialize arrays
4001
4003 for (cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
4004 {
4005 for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
4006 ++i, ++cnt)
4007 {
4008 exp2d = m_bndCondExpansions[n]
4009 ->GetExp(i)
4011
4013 m_graph->GetElementsFromFace(exp2d->GetGeom2D());
4014 m_BCtoElmMap[cnt] =
4015 globalIdMap[tmp->at(0).first->GetGlobalID()];
4016 m_BCtoTraceMap[cnt] = tmp->at(0).second;
4017 }
4018 }
4019 }
4020 break;
4021 default:
4022 ASSERTL1(false, "Not setup for this expansion");
4023 break;
4024 }
4025 }
4026
4027 ElmtID = m_BCtoElmMap;
4028 TraceID = m_BCtoTraceMap;
4029}
4030
4032 std::shared_ptr<ExpList> &result,
4033 const bool DeclareCoeffPhysArrays)
4034{
4035 int n, cnt;
4036 std::vector<unsigned int> eIDs;
4037
4038 Array<OneD, int> ElmtID, TraceID;
4039 GetBoundaryToElmtMap(ElmtID, TraceID);
4040
4041 // Skip other boundary regions
4042 for (cnt = n = 0; n < i; ++n)
4043 {
4044 cnt += m_bndCondExpansions[n]->GetExpSize();
4045 }
4046
4047 // Populate eIDs with information from BoundaryToElmtMap
4048 for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
4049 {
4050 eIDs.push_back(ElmtID[cnt + n]);
4051 }
4052
4053 // Create expansion list
4055 *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
4056 // Copy phys and coeffs to new explist
4057 if (DeclareCoeffPhysArrays)
4058 {
4059 int nq;
4060 int offsetOld, offsetNew;
4061 Array<OneD, NekDouble> tmp1, tmp2;
4062 for (n = 0; n < result->GetExpSize(); ++n)
4063 {
4064 nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
4065 offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
4066 offsetNew = result->GetPhys_Offset(n);
4067 Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
4068 tmp2 = result->UpdatePhys() + offsetNew, 1);
4069 nq = GetExp(ElmtID[cnt + n])->GetNcoeffs();
4070 offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
4071 offsetNew = result->GetCoeff_Offset(n);
4072 Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
4073 tmp2 = result->UpdateCoeffs() + offsetNew, 1);
4074 }
4075 }
4076}
4077
4078/**
4079 * @brief Reset this field, so that geometry information can be updated.
4080 */
4082{
4084
4085 // Reset boundary condition expansions.
4086 for (int n = 0; n < m_bndCondExpansions.size(); ++n)
4087 {
4088 m_bndCondExpansions[n]->Reset();
4089 }
4090}
4091
4092/**
4093 * Search through the edge expansions and identify which ones
4094 * have Robin/Mixed type boundary conditions. If find a Robin
4095 * boundary then store the edge id of the boundary condition
4096 * and the array of points of the physical space boundary
4097 * condition which are hold the boundary condition primitive
4098 * variable coefficient at the quatrature points
4099 *
4100 * \return std map containing the robin boundary condition
4101 * info using a key of the element id
4102 *
4103 * There is a next member to allow for more than one Robin
4104 * boundary condition per element
4105 */
4106map<int, RobinBCInfoSharedPtr> DisContField::v_GetRobinBCInfo(void)
4107{
4108 int i, cnt;
4109 map<int, RobinBCInfoSharedPtr> returnval;
4110 Array<OneD, int> ElmtID, TraceID;
4111 GetBoundaryToElmtMap(ElmtID, TraceID);
4112
4113 for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
4114 {
4116
4117 if (m_bndConditions[i]->GetBoundaryConditionType() ==
4119 {
4120 int e, elmtid;
4121 Array<OneD, NekDouble> Array_tmp;
4122
4123 locExpList = m_bndCondExpansions[i];
4124
4125 int npoints = locExpList->GetNpoints();
4126 Array<OneD, NekDouble> x0(npoints, 0.0);
4127 Array<OneD, NekDouble> x1(npoints, 0.0);
4128 Array<OneD, NekDouble> x2(npoints, 0.0);
4129 Array<OneD, NekDouble> coeffphys(npoints);
4130
4131 locExpList->GetCoords(x0, x1, x2);
4132
4133 LibUtilities::Equation coeffeqn =
4134 std::static_pointer_cast<
4136 ->m_robinPrimitiveCoeff;
4137
4138 // evalaute coefficient
4139 coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
4140
4141 for (e = 0; e < locExpList->GetExpSize(); ++e)
4142 {
4143 RobinBCInfoSharedPtr rInfo =
4145 TraceID[cnt + e],
4146 Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4147
4148 elmtid = ElmtID[cnt + e];
4149 // make link list if necessary
4150 if (returnval.count(elmtid) != 0)
4151 {
4152 rInfo->next = returnval.find(elmtid)->second;
4153 }
4154 returnval[elmtid] = rInfo;
4155 }
4156 }
4157 cnt += m_bndCondExpansions[i]->GetExpSize();
4158 }
4159
4160 return returnval;
4161}
4162
4163/**
4164 * @brief Calculate the \f$ L^2 \f$ error of the \f$ Q_{\rm dir} \f$
4165 * derivative using the consistent DG evaluation of \f$ Q_{\rm dir} \f$.
4166 *
4167 * The solution provided is of the primative variation at the quadrature
4168 * points and the derivative is compared to the discrete derivative at
4169 * these points, which is likely to be undesirable unless using a much
4170 * higher number of quadrature points than the polynomial order used to
4171 * evaluate \f$ Q_{\rm dir} \f$.
4172 */
4174 const Array<OneD, const NekDouble> &coeffs,
4175 const Array<OneD, const NekDouble> &soln)
4176{
4177 int i, e, ncoeff_edge;
4179 Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
4180
4182 m_traceMap->GetElmtToTrace();
4183
4185
4186 int cnt;
4187 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4188 Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
4189 m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4190
4191 edge_lambda = loc_lambda;
4192
4193 // Calculate Q using standard DG formulation.
4194 for (i = cnt = 0; i < GetExpSize(); ++i)
4195 {
4196 // Probably a better way of setting up lambda than this.
4197 // Note cannot use PutCoeffsInToElmts since lambda space
4198 // is mapped during the solve.
4199 int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4200 Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
4201
4202 for (e = 0; e < nEdges; ++e)
4203 {
4204 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4205 ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4206 edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
4207 Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4208 elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4209 edgeCoeffs[e]);
4210 edge_lambda = edge_lambda + ncoeff_edge;
4211 }
4212
4213 (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = coeffs + m_coeff_offset[i],
4214 elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4215 cnt += (*m_exp)[i]->GetNcoeffs();
4216 }
4217
4219 BwdTrans(out_d, phys);
4220 Vmath::Vsub(m_npoints, phys, 1, soln, 1, phys, 1);
4221 return L2(phys);
4222}
4223
4224/**
4225 * @brief Evaluate HDG post-processing to increase polynomial order of
4226 * solution.
4227 *
4228 * This function takes the solution (assumed to be one order lower) in
4229 * physical space, and postprocesses at the current polynomial order by
4230 * solving the system:
4231 *
4232 * \f[
4233 * \begin{aligned}
4234 * (\nabla w, \nabla u^*) &= (\nabla w, u), \\
4235 * \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle
4236 * \end{aligned}
4237 * \f]
4238 *
4239 * where \f$ u \f$ corresponds with the current solution as stored
4240 * inside #m_coeffs.
4241 *
4242 * @param outarray The resulting field \f$ u^* \f$.
4243 */
4245 const Array<OneD, const NekDouble> &coeffs,
4246 Array<OneD, NekDouble> &outarray)
4247{
4248 int i, cnt, e, ncoeff_trace;
4249 Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4251 m_traceMap->GetElmtToTrace();
4252
4254
4255 int nq_elmt, nm_elmt;
4256 int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4257 Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4258 Array<OneD, NekDouble> tmp_coeffs;
4259 m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4260
4261 trace_lambda = loc_lambda;
4262
4263 int dim = (m_expType == e2D) ? 2 : 3;
4264
4265 int num_points[] = {0, 0, 0};
4266 int num_modes[] = {0, 0, 0};
4267
4268 // Calculate Q using standard DG formulation.
4269 for (i = cnt = 0; i < GetExpSize(); ++i)
4270 {
4271 nq_elmt = (*m_exp)[i]->GetTotPoints();
4272 nm_elmt = (*m_exp)[i]->GetNcoeffs();
4273 qrhs = Array<OneD, NekDouble>(nq_elmt);
4274 qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4275 force = Array<OneD, NekDouble>(2 * nm_elmt);
4276 out_tmp = force + nm_elmt;
4278
4279 for (int j = 0; j < dim; ++j)
4280 {
4281 num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4282 num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4283 }
4284
4285 // Probably a better way of setting up lambda than this. Note
4286 // cannot use PutCoeffsInToElmts since lambda space is mapped
4287 // during the solve.
4288 int nTraces = (*m_exp)[i]->GetNtraces();
4289 Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
4290
4291 for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4292 {
4293 edgedir = (*m_exp)[i]->GetTraceOrient(e);
4294 ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4295 traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4296 Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4297 if (dim == 2)
4298 {
4299 elmtToTrace[i][e]->SetCoeffsToOrientation(
4300 edgedir, traceCoeffs[e], traceCoeffs[e]);
4301 }
4302 else
4303 {
4304 (*m_exp)[i]
4306 ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4307 }
4308 trace_lambda = trace_lambda + ncoeff_trace;
4309 }
4310
4311 // creating orthogonal expansion (checking if we have quads or
4312 // triangles)
4313 LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4314 switch (shape)
4315 {
4317 {
4318 const LibUtilities::PointsKey PkeyQ1(
4320 const LibUtilities::PointsKey PkeyQ2(
4323 num_modes[0], PkeyQ1);
4325 num_modes[1], PkeyQ2);
4327 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4328 (*m_exp)[i]->GetGeom());
4330 BkeyQ1, BkeyQ2, qGeom);
4331 }
4332 break;
4334 {
4335 const LibUtilities::PointsKey PkeyT1(
4337 const LibUtilities::PointsKey PkeyT2(
4338 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4340 num_modes[0], PkeyT1);
4342 num_modes[1], PkeyT2);
4344 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4345 (*m_exp)[i]->GetGeom());
4347 BkeyT1, BkeyT2, tGeom);
4348 }
4349 break;
4351 {
4352 const LibUtilities::PointsKey PkeyH1(
4354 const LibUtilities::PointsKey PkeyH2(
4356 const LibUtilities::PointsKey PkeyH3(
4359 num_modes[0], PkeyH1);
4361 num_modes[1], PkeyH2);
4363 num_modes[2], PkeyH3);
4365 std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4366 (*m_exp)[i]->GetGeom());
4368 BkeyH1, BkeyH2, BkeyH3, hGeom);
4369 }
4370 break;
4372 {
4373 const LibUtilities::PointsKey PkeyT1(
4375 const LibUtilities::PointsKey PkeyT2(
4376 num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4377 const LibUtilities::PointsKey PkeyT3(
4378 num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4380 num_modes[0], PkeyT1);
4382 num_modes[1], PkeyT2);
4384 num_modes[2], PkeyT3);
4386 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4387 (*m_exp)[i]->GetGeom());
4389 BkeyT1, BkeyT2, BkeyT3, tGeom);
4390 }
4391 break;
4392 default:
4394 "Wrong shape type, HDG postprocessing is not "
4395 "implemented");
4396 };
4397
4398 // DGDeriv
4399 // (d/dx w, d/dx q_0)
4400 (*m_exp)[i]->DGDeriv(0, tmp_coeffs = coeffs + m_coeff_offset[i],
4401 elmtToTrace[i], traceCoeffs, out_tmp);
4402 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4403 ppExp->IProductWRTDerivBase(0, qrhs, force);
4404
4405 // + (d/dy w, d/dy q_1)
4406 (*m_exp)[i]->DGDeriv(1, tmp_coeffs = coeffs + m_coeff_offset[i],
4407 elmtToTrace[i], traceCoeffs, out_tmp);
4408
4409 (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4410 ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4411
4412 Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4413
4414 // determine force[0] = (1,u)
4415 (*m_exp)[i]->BwdTrans(tmp_coeffs = coeffs + m_coeff_offset[i], qrhs);
4416 force[0] = (*m_exp)[i]->Integral(qrhs);
4417
4418 // multiply by inverse Laplacian matrix
4419 // get matrix inverse
4421 ppExp->DetShapeType(), *ppExp);
4422 DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4423
4424 NekVector<NekDouble> in(nm_elmt, force, eWrapper);
4425 NekVector<NekDouble> out(nm_elmt);
4426
4427 out = (*lapsys) * in;
4428
4429 // Transforming back to modified basis
4430 Array<OneD, NekDouble> work(nq_elmt);
4431 ppExp->BwdTrans(out.GetPtr(), work);
4432 (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4433 }
4434}
4435
4439 Array<OneD, NekDouble> &locTraceFwd, Array<OneD, NekDouble> &locTraceBwd)
4440{
4441 if (locTraceFwd != NullNekDouble1DArray)
4442 {
4443 m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
4444 locTraceFwd);
4445 }
4446
4447 if (locTraceBwd != NullNekDouble1DArray)
4448 {
4449 m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
4450 locTraceBwd);
4451 }
4452}
4453
4455 const Array<OneD, const NekDouble> &FwdFlux,
4456 const Array<OneD, const NekDouble> &BwdFlux,
4457 Array<OneD, NekDouble> &outarray)
4458{
4459 Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
4460
4461 m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4462 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
4463 m_trace->IProductWRTBase(BwdFlux, FCoeffs);
4464 m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
4465}
4466} // namespace Nektar::MultiRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
Describes the specification for a Basis.
Definition: Basis.h:45
Interpreter class for the evaluation of mathematical expressions.
Definition: Interpreter.h:76
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
Defines a specification for a set of points.
Definition: Points.h:50
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
This class is the abstractio n of a global discontinuous two- dimensional spectral/hp element expansi...
Definition: DisContField.h:56
~DisContField() override
Destructor.
Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions() override
std::vector< int > m_periodicBwdCopy
Definition: DisContField.h:198
void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray) override
Add trace contributions into elemental coefficient spaces.
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
Definition: DisContField.h:185
void SetUpDG(const std::string="DefaultVar", const Collections::ImplementationType ImpType=Collections::eNoImpType)
Set up all DG member variables and maps.
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
Definition: DisContField.h:190
std::set< int > m_boundaryTraces
A set storing the global IDs of any boundary Verts.
Definition: DisContField.h:175
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Definition: DisContField.h:210
std::vector< bool > m_negatedFluxNormal
Definition: DisContField.h:338
SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs(const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
AssemblyMapDGSharedPtr m_traceMap
Local to global DG mapping for trace space.
Definition: DisContField.h:170
Array< OneD, int > m_BCtoTraceMap
Definition: DisContField.h:59
void EvaluateHDGPostProcessing(const Array< OneD, const NekDouble > &coeffs, Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces) override
Obtain a copy of the periodic edges and vertices for this field.
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
For a given key, returns the associated global linear system.
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions() override
MultiRegions::ExpListSharedPtr & v_UpdateBndCondExpansion(int i) override
void v_SetBndCondBwdWeight(const int index, const NekDouble value) override
void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID) override
void v_ExtractTracePhys(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
This method extracts the trace (verts in 1D) from the field inarray and puts the values in outarray.
InterfaceMapDGSharedPtr & v_GetInterfaceMap(void) override
InterfaceMapDGSharedPtr m_interfaceMap
Interfaces mapping for trace space.
Definition: DisContField.h:161
void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays) override
void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp) override
Fill the weight with m_bndCondBndWeight.
const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight() override
const Array< OneD, const MultiRegions::ExpListSharedPtr > & v_GetBndCondExpansions() override
void FindPeriodicTraces(const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Generate a associative map of periodic vertices in a mesh.
Array< OneD, NekDouble > m_bndCondBndWeight
Definition: DisContField.h:158
bool SameTypeOfBoundaryConditions(const DisContField &In)
Check to see if expansion has the same BCs as In.
ExpListSharedPtr & v_GetTrace() override
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition structure definition on the diff...
Definition: DisContField.h:143
void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble) override
Evaluate all boundary conditions at a given time..
void FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndConditions, const Array< OneD, const ExpListSharedPtr > &BndCondExpansions, bool PutFwdInBwdOnBCs)
void GetFwdBwdTracePhys(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > &bndCond, const Array< OneD, const ExpListSharedPtr > &BndCondExp)
This method extracts the "forward" and "backward" trace data from the array field and puts the data i...
void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray) override
Add trace contributions into elemental coefficient spaces.
const LocTraceToTraceMapSharedPtr & v_GetLocTraceToTraceMap(void) const override
void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray) override
void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs) override
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
Definition: DisContField.h:164
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
Definition: DisContField.h:180
ExpListSharedPtr m_trace
Trace space storage for points between elements.
Definition: DisContField.h:167
void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field) override
void v_Reset() override
Reset this field, so that geometry information can be updated.
std::vector< bool > & GetNegatedFluxNormal(void)
DisContField()
Default constructor.
std::vector< bool > & v_GetLeftAdjacentTraces(void) override
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
Definition: DisContField.h:156
GlobalLinSysKey v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing) override
Solve the Helmholtz equation.
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
Definition: DisContField.h:197
bool IsLeftAdjacentTrace(const int n, const int e)
This routine determines if an element is to the "left" of the adjacent trace, which arises from the i...
std::vector< bool > m_leftAdjacentTraces
Definition: DisContField.h:204
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &coeffs, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd) override
AssemblyMapDGSharedPtr & v_GetTraceMap(void) override
std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo() override
void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd) override
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, const bool DeclareCoeffPhysArrays=true)
Discretises the boundary conditions.
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:99
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:3339
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
Definition: ExpList.h:1650
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2153
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1512
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2656
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2261
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:5772
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
Definition: ExpList.h:1944
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
Definition: ExpList.h:2030
int GetCoeff_Offset(int n) const
Get the start offset position for a local contiguous list of coeffs correspoinding to element n.
Definition: ExpList.h:2078
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2038
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1123
std::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:3093
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:2010
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1107
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1056
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1118
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1063
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2250
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2070
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:1060
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global This function calculates the error with respect to...
Definition: ExpList.h:514
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:3327
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1058
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
Definition: ExpList.h:1717
int GetPhys_Offset(int n) const
Get the start offset position for a local contiguous list of quadrature points in a full array corres...
Definition: ExpList.h:2085
ExpansionType m_expType
Expansion type.
Definition: ExpList.h:1051
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1099
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Describes a matrix with ordering defined by a local to global map.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:217
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:235
const BoundaryConditionCollection & GetBoundaryConditions(void) const
Definition: Conditions.h:245
void Rotate(PointGeom &a, int dir, NekDouble angle)
_this = rotation of a by angle 'angle' around axis dir
Definition: PointGeom.cpp:137
NekDouble dist(PointGeom &a)
return distance between this and input a
Definition: PointGeom.cpp:175
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2, bool doRot=false, int dir=0, NekDouble angle=0.0, NekDouble tol=1e-8)
Get the orientation of face1.
Definition: QuadGeom.cpp:129
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2, bool doRot, int dir, NekDouble angle, NekDouble tol)
Definition: TriGeom.cpp:109
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:46
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:46
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:46
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
Definition: Expansion0D.h:47
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:50
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:47
std::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:46
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:51
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1499
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
std::shared_ptr< InterfaceMapDG > InterfaceMapDGSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static GlobalLinSysKey NullGlobalLinSysKey(StdRegions::eNoMatrixType)
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
static const NekDouble kNekUnsetDouble
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:211
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshGraph.h:107
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:45
std::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:209
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:289
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:213
std::shared_ptr< std::vector< std::pair< GeometrySharedPtr, int > > > GeometryLinkSharedPtr
Definition: MeshGraph.h:166
std::map< int, std::vector< unsigned int > > BndRegionOrdering
Definition: MeshGraph.h:108
std::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
Definition: Conditions.h:214
std::shared_ptr< Composite > CompositeSharedPtr
Definition: MeshGraph.h:135
std::vector< PointGeomSharedPtr > PointGeomVector
Definition: Geometry2D.h:60
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:220
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:58
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Definition: Conditions.h:215
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:219
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:56
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:216
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:136
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
std::vector< double > w(NPUPPER)
std::vector< double > z(NPUPPER)
StdRegions::ConstFactorMap factors
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
static Array< OneD, NekDouble > NullNekDouble1DArray
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
int id
Geometry ID of entity.
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
bool isLocal
Flag specifying if this entity is local to this partition.
NekDouble m_tol
Tolerance to rotation is considered identical.
int m_dir
Axis of rotation. 0 = 'x', 1 = 'y', 2 = 'z'.
NekDouble m_angle
Angle of rotation in radians.