Nektar++
ExpList.cpp
Go to the documentation of this file.
1//////////////////////////////////////////////////////////////////////////////
2//
3// File: ExpList.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Expansion list definition
32//
33///////////////////////////////////////////////////////////////////////////////
34
40#include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
47#include <LocalRegions/HexExp.h>
48#include <LocalRegions/MatrixKey.h> // for MatrixKey
52#include <LocalRegions/PyrExp.h>
54#include <LocalRegions/SegExp.h>
55#include <LocalRegions/TetExp.h>
56#include <LocalRegions/TriExp.h>
57#include <MultiRegions/AssemblyMap/AssemblyMapCG.h> // for AssemblyMapCG, etc
58#include <MultiRegions/AssemblyMap/AssemblyMapDG.h> // for AssemblyMapDG, etc
59#include <MultiRegions/AssemblyMap/InterfaceMapDG.h> // for InterfaceMapDG, etc
62#include <MultiRegions/GlobalLinSysKey.h> // for GlobalLinSysKey
63#include <MultiRegions/GlobalMatrix.h> // for GlobalMatrix, etc
64#include <MultiRegions/GlobalMatrixKey.h> // for GlobalMatrixKey
65#include <iomanip>
66
67using namespace std;
68
70{
71/**
72 * @class ExpList
73 * All multi-elemental expansions \f$u^{\delta}(\boldsymbol{x})\f$ can
74 * be considered as the assembly of the various elemental contributions.
75 * On a discrete level, this yields,
76 * \f[u^{\delta}(\boldsymbol{x}_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
77 * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(\boldsymbol{x}_i).\f]
78 * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
79 * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
80 * As it is the lowest level class, it contains the definition of the
81 * common data and common routines to all multi-elemental expansions.
82 *
83 * The class stores a vector of expansions, \a m_exp, (each derived from
84 * StdRegions#StdExpansion) which define the constituent components of
85 * the domain.
86 */
87
88/**
89 * Creates an empty expansion list.
90 */
92 : m_expType(type), m_ncoeffs(0), m_npoints(0), m_physState(false),
93 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
94 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
95 m_WaveSpace(false)
96{
97}
98
99/*----------------------------------------------------------------*/
100/* Copy Construtor */
101/*-----------------------------------------------------------------*/
102
103/**
104 * Copies an existing expansion list.
105 * @param in Source expansion list.
106 */
107ExpList::ExpList(const ExpList &in, const bool DeclareCoeffPhysArrays)
108 : std::enable_shared_from_this<ExpList>(in), m_expType(in.m_expType),
109
110 m_comm(in.m_comm), m_session(in.m_session), m_graph(in.m_graph),
111 m_ncoeffs(in.m_ncoeffs), m_npoints(in.m_npoints), m_physState(false),
112 m_exp(in.m_exp), m_collections(in.m_collections),
113 m_collectionsDoInit(in.m_collectionsDoInit),
114 m_coeff_offset(in.m_coeff_offset), m_phys_offset(in.m_phys_offset),
115 m_blockMat(in.m_blockMat), m_WaveSpace(false),
116 m_elmtToExpId(in.m_elmtToExpId)
117{
118 // Set up m_coeffs, m_phys and offset arrays.
119 // use this to keep memory declaration in one place
120 SetupCoeffPhys(DeclareCoeffPhysArrays, false);
121}
122
123/**
124 * Copies the eIds elements from an existing expansion list.
125 * @param in Source expansion list.
126 * @param in elements that will be in the new exp list.
127 */
128ExpList::ExpList(const ExpList &in, const std::vector<unsigned int> &eIDs,
129 const bool DeclareCoeffPhysArrays,
131 : m_expType(in.m_expType), m_comm(in.m_comm), m_session(in.m_session),
132 m_graph(in.m_graph), m_physState(false),
133 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
134 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
135 m_WaveSpace(false)
136{
137 for (int i = 0; i < eIDs.size(); ++i)
138 {
139 (*m_exp).push_back((*(in.m_exp))[eIDs[i]]);
140 }
141
142 // Set up m_coeffs, m_phys and offset arrays.
143 SetupCoeffPhys(DeclareCoeffPhysArrays);
144
145 // set up collections
146 CreateCollections(ImpType);
147}
148
149/**
150 * Given a meshgraph \a graph, containing information about
151 * the domain and the spectral/hp element expansion, this
152 * constructor fills the list of local expansions
153 * \texttt{m_exp} with the proper expansions, calculates the
154 * total number of quadrature points \f$x_i\f$ and local
155 * expansion coefficients \f$\hat{u}^e_n\f$ and
156 *
157 * @param pSession A session within information about expansion
158 *
159 * @param graph A meshgraph, containing information about the
160 * domain and the spectral/hp element expansion.
161 *
162 * @param DeclareCoeffPhysArrays Declare the coefficient and
163 * phys space arrays
164 *
165 * @param ImpType Detail about the implementation type to use
166 * in operators
167 */
170 const bool DeclareCoeffPhysArrays, const std::string &var,
172 : m_comm(pSession->GetComm()), m_session(pSession), m_graph(graph),
173 m_physState(false),
174 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
175 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
176 m_WaveSpace(false)
177{
178 // Retrieve the list of expansions
179 const SpatialDomains::ExpansionInfoMap &expansions =
180 graph->GetExpansionInfo(var);
181
182 // Initialise Expansion Vector
183 InitialiseExpVector(expansions);
184
185 // Setup phys coeff space
186 SetupCoeffPhys(DeclareCoeffPhysArrays);
187
188 // Initialise collection
189 CreateCollections(ImpType);
190}
191
192/**
193 * Given an expansion vector \a expansions, containing
194 * information about the domain and the spectral/hp element
195 * expansion, this constructor fills the list of local
196 * expansions \texttt{m_exp} with the proper expansions,
197 * calculates the total number of quadrature points
198 * \f$\boldsymbol{x}_i\f$ and local expansion coefficients
199 * \f$\hat{u}^e_n\f$.
200 *
201 * @param pSession A session within information about expansion
202 * @param expansions A vector containing information about the
203 * domain and the spectral/hp element
204 * expansion.
205 * @param DeclareCoeffPhysArrays Declare the coefficient and
206 * phys space arrays
207 * @param ImpType Detail about the implementation type to use
208 * in operators
209 */
211 const SpatialDomains::ExpansionInfoMap &expansions,
212 const bool DeclareCoeffPhysArrays,
214 : m_comm(pSession->GetComm()), m_session(pSession), m_physState(false),
215 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
216 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
217 m_WaveSpace(false)
218{
219 // Initialise expansion vector
220 InitialiseExpVector(expansions);
221
222 // Set up m_coeffs, m_phys and offset arrays.
223 SetupCoeffPhys(DeclareCoeffPhysArrays);
224
225 // Setup Collection
226 CreateCollections(ImpType);
227}
228
229//----------------------------------------------------------------------
230// 0D Expansion Constructors
231//----------------------------------------------------------------------
233 : m_expType(e0D), m_ncoeffs(1), m_npoints(1), m_physState(false),
234 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
235 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
236 m_WaveSpace(false)
237{
240 (*m_exp).push_back(Point);
241
243}
244
245/**
246 * Store expansions for the trace space expansions used in
247 * DisContField
248 *
249 * @param pSession A session within information about expansion
250 * @param bndConstraint Array of ExpList1D objects each containing a
251 * 1D spectral/hp element expansion on a single
252 * boundary region.
253 * @param bndCond Array of BoundaryCondition objects which contain
254 * information about the boundary conditions on the
255 * different boundary regions.
256 * @param locexp Complete domain expansion list.
257 * @param graph mesh corresponding to the expansion list.
258 * @param DeclareCoeffPhysArrays Declare the coefficient and
259 * phys space arrays
260 * @param variable The variable name associated with the expansion
261 * @param ImpType Detail about the implementation type to use
262 * in operators
263 *
264 * By default the routine will try and order the expansions in a
265 * manner which is optimal for collection type operations. This can be
266 * disabled by the command line option --no-exp-opt
267 */
270 const Array<OneD, const ExpListSharedPtr> &bndConstraint,
272 const LocalRegions::ExpansionVector &locexp,
274 const LibUtilities::CommSharedPtr &comm, const bool DeclareCoeffPhysArrays,
275 [[maybe_unused]] const std::string variable,
276 [[maybe_unused]] const Collections::ImplementationType ImpType)
277 : m_comm(comm), m_session(pSession), m_graph(graph), m_physState(false),
278 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
279 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
280 m_WaveSpace(false)
281{
282 int i, j, id, elmtid = 0;
283 set<int> tracesDone;
284
291
297
298 map<int, vector<SpatialDomains::ExpansionInfoShPtr>> ExpOrder;
300
301 bool DoOptOnCollection =
302 m_session->DefinesCmdLineArgument("no-exp-opt") ? false : true;
303 int cnt = 0;
304 for (i = 0; i < bndCond.size(); ++i)
305 {
306 if (bndCond[i]->GetBoundaryConditionType() ==
308 {
309 bool IsNot0D = true; // Cehck for 0D expansion
310 for (j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
311 {
315 bndConstraint[i]->GetExp(j)->GetGeom(), PtBvec);
316
317 if ((exp1D =
318 std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
319 bndConstraint[i]->GetExp(j))))
320 {
322 exp1D->GetBasis(0)->GetBasisKey();
323 eInfo->m_basisKeyVector.push_back(bkey);
324 }
325 else if ((exp2D = std::dynamic_pointer_cast<
327 bndConstraint[i]->GetExp(j))))
328 {
330 exp2D->GetBasis(0)->GetBasisKey();
332 exp2D->GetBasis(1)->GetBasisKey();
333
334 eInfo->m_basisKeyVector.push_back(bkey0);
335 eInfo->m_basisKeyVector.push_back(bkey1);
336 }
337 else // 0D case no optimisation necessary
338 {
339 IsNot0D = false;
340 }
341
342 if (DoOptOnCollection && IsNot0D)
343 {
344 int i;
345 for (i = 0; i < cnt; ++i)
346 {
347 if ((eInfo->m_basisKeyVector ==
348 ExpOrder[i][0]->m_basisKeyVector) &&
349 (eInfo->m_geomShPtr->GetGeomFactors()->GetGtype() ==
350 ExpOrder[i][0]
351 ->m_geomShPtr->GetGeomFactors()
352 ->GetGtype()))
353 {
354 ExpOrder[i].push_back(eInfo);
355 break;
356 }
357 }
358
359 if (i == cnt)
360 {
361 ExpOrder[cnt++].push_back(eInfo);
362 }
363 }
364 else
365 {
366 ExpOrder[0].push_back(eInfo);
367 }
368 }
369 }
370 }
371
372 // decalare expansions in provided order using geom and basis info
373 for (auto &ordIt : ExpOrder)
374 {
375 for (auto &eit : ordIt.second)
376 {
377
378 if ((PointGeom =
379 std::dynamic_pointer_cast<SpatialDomains::PointGeom>(
380 eit->m_geomShPtr)))
381 {
382 m_expType = e0D;
383
385 PointGeom);
386 tracesDone.insert(PointGeom->GetVid());
387 }
388 else if ((segGeom =
389 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
390 eit->m_geomShPtr)))
391 {
392 m_expType = e1D;
393
395 eit->m_basisKeyVector[0], segGeom);
396 tracesDone.insert(segGeom->GetGlobalID());
397 }
398 else if ((TriGeom =
399 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
400 eit->m_geomShPtr)))
401 {
402 m_expType = e2D;
403
405 eit->m_basisKeyVector[0], eit->m_basisKeyVector[1],
406 TriGeom);
407
408 tracesDone.insert(TriGeom->GetGlobalID());
409 }
410 else if ((QuadGeom =
411 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
412 eit->m_geomShPtr)))
413 {
414 m_expType = e2D;
416 eit->m_basisKeyVector[0], eit->m_basisKeyVector[1],
417 QuadGeom);
418
419 tracesDone.insert(QuadGeom->GetGlobalID());
420 }
421
422 // Assign next id
423 exp->SetElmtId(elmtid++);
424
425 // Add the expansion
426 (*m_exp).push_back(exp);
427 }
428 }
429
430 map<int, pair<SpatialDomains::Geometry1DSharedPtr, LibUtilities::BasisKey>>
431 edgeOrders;
432
434 pair<LibUtilities::BasisKey, LibUtilities::BasisKey>>>
435 faceOrders;
436
437 for (i = 0; i < locexp.size(); ++i)
438 {
439 if ((exp1D = std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
440 locexp[i])))
441 {
442 m_expType = e0D;
443
444 for (j = 0; j < 2; ++j)
445 {
446 PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
447 id = PointGeom->GetVid();
448
449 // Ignore Dirichlet edges
450 if (tracesDone.count(id) != 0)
451 {
452 continue;
453 }
454
456 PointGeom);
457 tracesDone.insert(id);
458 exp->SetElmtId(elmtid++);
459 (*m_exp).push_back(exp);
460 }
461 }
462 else if ((exp2D = std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
463 locexp[i])))
464 {
465 m_expType = e1D;
466 for (j = 0; j < locexp[i]->GetNtraces(); ++j)
467 {
468 segGeom = exp2D->GetGeom2D()->GetEdge(j);
469 id = segGeom->GetGlobalID();
470 // Ignore Dirichlet edges
471 if (tracesDone.count(id) != 0)
472 {
473 continue;
474 }
475
476 auto it = edgeOrders.find(id);
477
478 if (it == edgeOrders.end()) // not exist, insert new one
479 {
480 edgeOrders.insert(std::make_pair(
481 id, std::make_pair(segGeom,
482 locexp[i]->GetTraceBasisKey(j))));
483 }
484 else // already exist, check if we need to update it
485 {
486 LibUtilities::BasisKey edge =
487 locexp[i]->GetTraceBasisKey(j);
488 LibUtilities::BasisKey existing = it->second.second;
489
490 int np1 = edge.GetNumPoints();
491 int np2 = existing.GetNumPoints();
492 int nm1 = edge.GetNumModes();
493 int nm2 = existing.GetNumModes();
494
495 // if the existing edge has less points/modes than the
496 // present edge, then we update the existing edge with
497 // present one (trace should always have highest order)
498
499 // The pointsType is always GLL for edges
500 // So we can directly compare them.
501 if (np2 >= np1 && nm2 >= nm1)
502 {
503 continue;
504 }
505 else if (np2 <= np1 && nm2 <= nm1)
506 {
507 it->second.second = edge;
508 }
509 else
510 {
512 "inappropriate number of points/modes (max"
513 "num of points is not set with max order)");
514 }
515 }
516 }
517 }
518 else if ((exp3D = dynamic_pointer_cast<LocalRegions::Expansion3D>(
519 locexp[i])))
520 {
521 m_expType = e2D;
522 for (j = 0; j < exp3D->GetNtraces(); ++j)
523 {
524 FaceGeom = exp3D->GetGeom3D()->GetFace(j);
525 id = FaceGeom->GetGlobalID();
526 // Ignore Dirichlet edges
527 if (tracesDone.count(id) != 0)
528 {
529 continue;
530 }
531 auto it = faceOrders.find(id);
532
533 if (it == faceOrders.end()) // not exist, insert new one
534 {
535 // GetTraceBasisKey does not take into account
536 // the orientaion. It is w.r.t local face axes.
537 LibUtilities::BasisKey face_dir0 =
538 locexp[i]->GetTraceBasisKey(j, 0);
539 LibUtilities::BasisKey face_dir1 =
540 locexp[i]->GetTraceBasisKey(j, 1);
541
542 // If the axes 1, 2 of the local face correspond to
543 // the axes 2, 1 of global face, respectively,
544 // then swap face0 and face1
545 if (locexp[i]->GetTraceOrient(j) >= 9)
546 {
547 std::swap(face_dir0, face_dir1);
548 }
549
550 faceOrders.insert(std::make_pair(
551 id,
552 std::make_pair(FaceGeom,
553 std::make_pair(face_dir0, face_dir1))));
554 }
555 else // already exist, check if we need to update it
556 {
557 LibUtilities::BasisKey face0 =
558 locexp[i]->GetTraceBasisKey(j, 0);
559 LibUtilities::BasisKey face1 =
560 locexp[i]->GetTraceBasisKey(j, 1);
561 LibUtilities::BasisKey existing0 = it->second.second.first;
562 LibUtilities::BasisKey existing1 = it->second.second.second;
563
564 // np -- number of points; nm -- number of modes;
565 // np_I_J --- I=1 current; I=2 existing; J=1 dir0; J=2 dir1;
566 int np11 = face0.GetNumPoints();
567 int np12 = face1.GetNumPoints();
568 int np21 = existing0.GetNumPoints();
569 int np22 = existing1.GetNumPoints();
570 int nm11 = face0.GetNumModes();
571 int nm12 = face1.GetNumModes();
572 int nm21 = existing0.GetNumModes();
573 int nm22 = existing1.GetNumModes();
574
575 // If the axes 1, 2 of the current face correspond to
576 // the axes 2, 1 of existing face, respectively,
577 // then swap face0 and face1
578 // eDir1FwdDir2_Dir2FwdDir1 = 9
579 // eDir1FwdDir2_Dir2BwdDir1 = 10
580 // eDir1BwdDir2_Dir2FwdDir1 = 11
581 // eDir1BwdDir2_Dir2BwdDir1 = 12
582 if (locexp[i]->GetTraceOrient(j) >= 9)
583 {
584 std::swap(np11, np12);
585 std::swap(nm11, nm12);
586 std::swap(face0, face1);
587 }
588
589 // The baiskey return by GetTraceBasisKey should always
590 // have GLL for eModified_A and GR for eModified_B.
591 // But we still use GetPointsType to check this.
592 if (existing1.GetPointsType() ==
593 LibUtilities::eGaussRadauMAlpha1Beta0)
594 {
595 if (face1.GetPointsType() ==
597 {
598 np12--; // make np12 comparable to np22
599 }
600 }
601 else
602 {
603 if (face1.GetPointsType() ==
604 LibUtilities::eGaussRadauMAlpha1Beta0)
605 {
606 np12++; // make np12 comparable to np22
607 }
608 }
609
610 if (existing0.GetPointsType() ==
611 LibUtilities::eGaussRadauMAlpha1Beta0)
612 {
613 if (face0.GetPointsType() ==
615 {
616 np11--; // make np11 comparable to np21
617 }
618 }
619 else
620 {
621 if (face0.GetPointsType() ==
622 LibUtilities::eGaussRadauMAlpha1Beta0)
623 {
624 np11++; // make np11 comparable to np21
625 }
626 }
627
628 // if the existing face_i has less points/modes than the
629 // present face_i, then we update the existing face_i with
630 // present one (trace should always have highest order)
631 if (np22 >= np12 && nm22 >= nm12)
632 {
633 // keep existing face_i and do nothing
634 }
635 else if (np22 <= np12 && nm22 <= nm12)
636 {
637 // Instead of using face0 directly, We create new
638 // basiskey with original Type but higher order.
639 LibUtilities::BasisKey newbkey(
640 existing1.GetBasisType(), nm12,
641 LibUtilities::PointsKey(np12,
642 existing1.GetPointsType()));
643 it->second.second.second = newbkey;
644 }
645 else // np22 > np12 but nm22 < nm12
646 {
648 "inappropriate number of points/modes (max"
649 "num of points is not set with max order)");
650 }
651
652 if (np21 >= np11 && nm21 >= nm11)
653 {
654 // keep existing face_i and do nothing
655 }
656 else if (np21 <= np11 && nm21 <= nm11)
657 {
658 // Instead of using face0 directly, We create new
659 // basiskey with original Type but higher order.
660 LibUtilities::PointsKey newpkey(
661 np11, existing0.GetPointsType());
662 LibUtilities::BasisKey newbkey(existing0.GetBasisType(),
663 nm11, newpkey);
664 it->second.second.first = newbkey;
665 }
666 else // np21 > np11 but nm21 < nm11
667 {
669 "inappropriate number of points/modes (max"
670 "num of points is not set with max order)");
671 }
672 }
673 }
674 }
675 }
676
677 int nproc = m_comm->GetRowComm()->GetSize(); // number of processors
678 int tracepr = m_comm->GetRowComm()->GetRank(); // ID processor
679
680 if (nproc > 1)
681 {
682 int tCnt = 0;
683
684 // Count the number of traces on each partition
685 for (i = 0; i < locexp.size(); ++i)
686 {
687 tCnt += locexp[i]->GetNtraces();
688 }
689
690 // Record the number traces on each partition,
691 // then reduce this across processors.
692 Array<OneD, int> tracesCnt(nproc, 0);
693 tracesCnt[tracepr] = tCnt;
694 m_comm->GetRowComm()->AllReduce(tracesCnt, LibUtilities::ReduceSum);
695
696 // Set up offset array.
697 int totTraceCnt = Vmath::Vsum(nproc, tracesCnt, 1);
698 Array<OneD, int> tTotOffsets(nproc, 0);
699
700 for (i = 1; i < nproc; ++i)
701 {
702 tTotOffsets[i] = tTotOffsets[i - 1] + tracesCnt[i - 1];
703 }
704
705 // Set up arrays that are large enough to hold all traces
706 // in the domain, not just those in local partition.
707 Array<OneD, int> TracesTotID(totTraceCnt, 0);
708 Array<OneD, int> TracesTotNm0(totTraceCnt, 0);
709 Array<OneD, int> TracesTotNm1(totTraceCnt, 0);
710 Array<OneD, int> TracesTotPnts0(totTraceCnt, 0);
711 Array<OneD, int> TracesTotPnts1(totTraceCnt, 0);
712 Array<OneD, int> TracesPointsType0(totTraceCnt, 0);
713 Array<OneD, int> TracesPointsType1(totTraceCnt, 0);
714 // convert enum PointsType to integer, so that MPI can handle it
715
716 // TODO: Change this design to save memory and improve performance.
717 //
718 // In the current design, After AllReduce, each processor will
719 // store a full copy these large arrays, which might be too
720 // expensive when nproc and mesh are large. So basically, this
721 // design is just a naive implementation and we need to redesign it
722 // in the future.
723
724 int cntr = tTotOffsets[tracepr];
725
726 for (i = 0; i < locexp.size(); ++i)
727 {
728 if ((exp2D = locexp[i]->as<LocalRegions::Expansion2D>()))
729 {
730 int nedges = locexp[i]->GetNtraces();
731
732 for (j = 0; j < nedges; ++j, ++cntr)
733 {
734 LibUtilities::BasisKey bkeyEdge =
735 locexp[i]->GetTraceBasisKey(j);
736 TracesTotID[cntr] = exp2D->GetGeom2D()->GetEid(j);
737 TracesTotNm0[cntr] = bkeyEdge.GetNumModes();
738 TracesTotPnts0[cntr] = bkeyEdge.GetNumPoints();
739 TracesPointsType0[cntr] =
740 static_cast<int>(bkeyEdge.GetPointsType());
741 // Althought for edges, we only have GLL points, we still
742 // keep this design for backup.
743 }
744 }
745 else if ((exp3D = locexp[i]->as<LocalRegions::Expansion3D>()))
746 {
747 int nfaces = locexp[i]->GetNtraces();
748
749 for (j = 0; j < nfaces; ++j, ++cntr)
750 {
751 LibUtilities::BasisKey face_dir0 =
752 locexp[i]->GetTraceBasisKey(j, 0);
753 LibUtilities::BasisKey face_dir1 =
754 locexp[i]->GetTraceBasisKey(j, 1);
755
756 // If the axes 1, 2 of the local face correspond to
757 // the axes 2, 1 of global face, respectively,
758 // then swap face0 and face1
759 if (locexp[i]->GetTraceOrient(j) >= 9)
760 {
761 std::swap(face_dir0, face_dir1);
762 }
763
764 TracesTotID[cntr] = exp3D->GetGeom3D()->GetFid(j);
765 TracesTotNm0[cntr] = face_dir0.GetNumModes();
766 TracesTotNm1[cntr] = face_dir1.GetNumModes();
767 TracesTotPnts0[cntr] = face_dir0.GetNumPoints();
768 TracesTotPnts1[cntr] = face_dir1.GetNumPoints();
769 TracesPointsType0[cntr] =
770 static_cast<int>(face_dir0.GetPointsType());
771 TracesPointsType1[cntr] =
772 static_cast<int>(face_dir1.GetPointsType());
773 }
774 }
775 }
776
777 m_comm->GetRowComm()->AllReduce(TracesTotID, LibUtilities::ReduceSum);
778 m_comm->GetRowComm()->AllReduce(TracesTotNm0, LibUtilities::ReduceSum);
779 m_comm->GetRowComm()->AllReduce(TracesTotPnts0,
781 m_comm->GetRowComm()->AllReduce(TracesPointsType0,
783 if (m_expType == e2D) // 2D face
784 {
785 m_comm->GetRowComm()->AllReduce(TracesTotNm1,
787 m_comm->GetRowComm()->AllReduce(TracesTotPnts1,
789 m_comm->GetRowComm()->AllReduce(TracesPointsType1,
791 }
792 // TracesTotXXX has collected traces info of entire domain
793 // Now compare them with current traces info - edgeOrders
794 // and update traces basiskey. This step is to ensure the
795 // m_trace has highest order even on parition interfaces.
796
797 if (edgeOrders.size())
798 {
799 for (i = 0; i < totTraceCnt; ++i)
800 {
801 auto it = edgeOrders.find(TracesTotID[i]);
802
803 if (it == edgeOrders.end())
804 {
805 continue;
806 }
807
808 LibUtilities::BasisKey existing = it->second.second;
809
810 int np1 = TracesTotPnts0[i];
811 int np2 = existing.GetNumPoints();
812 int nm1 = TracesTotNm0[i];
813 int nm2 = existing.GetNumModes();
814
815 // The pointsType is always GLL for edges
816 // So we can directly compare them.
817 if (np2 >= np1 && nm2 >= nm1)
818 {
819 continue;
820 }
821 else if (np2 <= np1 && nm2 <= nm1)
822 {
823 // We have to rebuild the basis key locally because
824 // MPI::AllReduce does not support BasisKey directly.
825 LibUtilities::BasisKey newbkey(
826 existing.GetBasisType(), nm1,
827 LibUtilities::PointsKey(np1, existing.GetPointsType()));
828 it->second.second = newbkey;
829 }
830 else
831 {
833 "inappropriate number of points/modes (max "
834 "num of points is not set with max order)");
835 }
836 }
837 }
838 else if (faceOrders.size())
839 {
840 for (i = 0; i < totTraceCnt; ++i)
841 {
842 auto it = faceOrders.find(TracesTotID[i]);
843
844 if (it == faceOrders.end())
845 {
846 continue;
847 }
848
849 LibUtilities::BasisKey existing0 = it->second.second.first;
850 LibUtilities::BasisKey existing1 = it->second.second.second;
851
852 // np -- number of points; nm -- number of modes;
853 // np_I_J --- I=1 current; I=2 existing; J=1 dir0; J=2 dir1;
854 int np11 = TracesTotPnts0[i];
855 int np12 = TracesTotPnts1[i];
856 int np21 = existing0.GetNumPoints();
857 int np22 = existing1.GetNumPoints();
858 int nm11 = TracesTotNm0[i];
859 int nm12 = TracesTotNm1[i];
860 int nm21 = existing0.GetNumModes();
861 int nm22 = existing1.GetNumModes();
862
863 // The orientation is already aligned
864 // Here we only need to compare pointsType
865 // and adjust np
866 if (existing1.GetPointsType() ==
867 LibUtilities::eGaussRadauMAlpha1Beta0)
868 {
869 if (static_cast<LibUtilities::PointsType>(
870 TracesPointsType1[i]) ==
872 {
873 np12--; // make np12 comparable to np22
874 }
875 }
876 else
877 {
878 if (static_cast<LibUtilities::PointsType>(
879 TracesPointsType1[i]) ==
880 LibUtilities::eGaussRadauMAlpha1Beta0)
881 {
882 np12++; // make np12 comparable to np22
883 }
884 }
885
886 if (existing0.GetPointsType() ==
887 LibUtilities::eGaussRadauMAlpha1Beta0)
888 {
889 if (static_cast<LibUtilities::PointsType>(
890 TracesPointsType0[i]) ==
892 {
893 np11--; // make np11 comparable to np21
894 }
895 }
896 else
897 {
898 if (static_cast<LibUtilities::PointsType>(
899 TracesPointsType0[i]) ==
900 LibUtilities::eGaussRadauMAlpha1Beta0)
901 {
902 np11++; // make np11 comparable to np21
903 }
904 }
905
906 // if the existing face_i has less points/modes than the
907 // present face_i, then we update the existing face_i with
908 // present one (trace should always have highest order)
909 if (np22 >= np12 && nm22 >= nm12)
910 {
911 // keep existing face_i and do nothing
912 }
913 else if (np22 <= np12 && nm22 <= nm12)
914 {
915 LibUtilities::BasisKey newbkey(
916 existing1.GetBasisType(), nm12,
917 LibUtilities::PointsKey(np12,
918 existing1.GetPointsType()));
919 it->second.second.second = newbkey;
920 }
921 else // np22 > np12 but nm22 < nm12
922 {
924 "inappropriate number of points/modes (max "
925 "num of points is not set with max order)");
926 }
927
928 if (np21 >= np11 && nm21 >= nm11)
929 {
930 // keep existing face_i and do nothing
931 }
932 else if (np21 <= np11 && nm21 <= nm11)
933 {
934 LibUtilities::BasisKey newbkey(
935 existing0.GetBasisType(), nm11,
936 LibUtilities::PointsKey(np11,
937 existing0.GetPointsType()));
938 it->second.second.first = newbkey;
939 }
940 else // np21 > np11 but nm21 < nm11
941 {
943 "inappropriate number of points/modes (max"
944 "num of points is not set with max order)");
945 }
946 }
947 }
948 }
949
950 if (edgeOrders.size())
951 {
952 map<int, vector<int>> opt;
953 int cnt = 0;
954
955 if (DoOptOnCollection)
956 {
957 for (auto &it : edgeOrders)
958 {
959 int i;
960 for (i = 0; i < cnt; ++i)
961 {
962 auto it1 = edgeOrders.find(opt[i][0]);
963
964 if ((it.second.second == it1->second.second) &&
965 (it.second.first->GetGeomFactors()->GetGtype() ==
966 it1->second.first->GetGeomFactors()->GetGtype()))
967 {
968 opt[i].push_back(it.first);
969 break;
970 }
971 }
972
973 if (i == cnt)
974 {
975 opt[cnt++].push_back(it.first);
976 }
977 }
978 }
979 else
980 {
981 for (auto &it : edgeOrders)
982 {
983 opt[0].push_back(it.first);
984 }
985 }
986
987 for (int i = 0; i < opt.size(); ++i)
988 {
989
990 for (int j = 0; j < opt[i].size(); ++j)
991 {
992 auto it = edgeOrders.find(opt[i][j]);
993
995 it->second.second, it->second.first);
996
997 exp->SetElmtId(elmtid++);
998 (*m_exp).push_back(exp);
999 }
1000 }
1001 }
1002 else
1003 {
1004 map<int, vector<int>> opt;
1005 int cnt = 0;
1006
1007 if (DoOptOnCollection)
1008 {
1009 for (auto &it : faceOrders)
1010 {
1011 int i;
1012 for (i = 0; i < cnt; ++i)
1013 {
1014 auto it1 = faceOrders.find(opt[i][0]);
1015
1016 if ((it.second.second.first == it1->second.second.first) &&
1017 (it.second.second.second ==
1018 it1->second.second.second) &&
1019 (it.second.first->GetGeomFactors()->GetGtype() ==
1020 it1->second.first->GetGeomFactors()->GetGtype()))
1021 {
1022 opt[i].push_back(it.first);
1023 break;
1024 }
1025 }
1026
1027 if (i == cnt)
1028 {
1029 opt[cnt++].push_back(it.first);
1030 }
1031 }
1032 }
1033 else
1034 {
1035 for (auto &it : faceOrders)
1036 {
1037 opt[0].push_back(it.first);
1038 }
1039 }
1040
1041 for (int i = 0; i < opt.size(); ++i)
1042 {
1043
1044 for (int j = 0; j < opt[i].size(); ++j)
1045 {
1046 auto it = faceOrders.find(opt[i][j]);
1047
1048 FaceGeom = it->second.first;
1049
1050 if ((QuadGeom =
1051 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
1052 FaceGeom)))
1053 {
1054 exp =
1056 it->second.second.first, it->second.second.second,
1057 QuadGeom);
1058 }
1059 else if ((TriGeom = std::dynamic_pointer_cast<
1060 SpatialDomains::TriGeom>(FaceGeom)))
1061 {
1062 exp =
1064 it->second.second.first, it->second.second.second,
1065 TriGeom);
1066 }
1067 exp->SetElmtId(elmtid++);
1068 (*m_exp).push_back(exp);
1069 }
1070 }
1071 }
1072
1073 // Set up m_coeffs, m_phys and offset arrays.
1074 SetupCoeffPhys(DeclareCoeffPhysArrays);
1075
1076 // Set up collections
1077 if (m_expType != e0D)
1078 {
1079 CreateCollections(ImpType);
1080 }
1081
1082 // Setup element to expansion ID maps for the trace elements
1083 // Loop in reverse order so that in case where using a
1084 // Homogeneous expansion it sets geometry ids to first part of
1085 // m_exp list. Otherwise will set to second (complex) expansion
1086 for (int i = (*m_exp).size() - 1; i >= 0; --i)
1087 {
1088 m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1089 }
1090}
1091
1092/**
1093 * Set expansions for localtrace space expansions used in
1094 * DisContField as part of Gradient Jump Penalisation
1095 *
1096 * @param pSession A session within information about expansion
1097 * @param locexp Complete domain expansion list.
1098 * @param graph mesh corresponding to the expansion list.
1099 * @param DeclareCoeffPhysArrays Declare the coefficient and
1100 * phys space arrays
1101 * @param variable The variable name associated with the expansion
1102 * @param ImpType Detail about the implementation type to use
1103 * in operators
1104 */
1106 const LocalRegions::ExpansionVector &locexp,
1108 const bool DeclareCoeffPhysArrays,
1109 [[maybe_unused]] const std::string variable,
1110 [[maybe_unused]] const Collections::ImplementationType ImpType)
1111 : m_comm(pSession->GetComm()), m_session(pSession), m_graph(graph),
1112 m_physState(false),
1113 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
1114 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
1115 m_WaveSpace(false)
1116{
1117 int i, j, elmtid = 0;
1118
1125
1131
1132 for (i = 0; i < locexp.size(); ++i)
1133 {
1134 if ((exp1D = std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
1135 locexp[i])))
1136 {
1137 m_expType = e0D;
1138
1139 for (j = 0; j < 2; ++j)
1140 {
1141 PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
1142
1144 PointGeom);
1145 exp->SetElmtId(elmtid++);
1146 (*m_exp).push_back(exp);
1147 }
1148 }
1149 else if ((exp2D = std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
1150 locexp[i])))
1151 {
1152 m_expType = e1D;
1153 LibUtilities::BasisKey edgeKey0 =
1154 locexp[i]->GetBasis(0)->GetBasisKey();
1155
1156 for (j = 0; j < locexp[i]->GetNtraces(); ++j)
1157 {
1158 segGeom = exp2D->GetGeom2D()->GetEdge(j);
1159
1160 int dir = exp2D->GetGeom2D()->GetDir(j);
1161
1162 if (locexp[i]->GetNtraces() == 3)
1163 {
1164 LibUtilities::BasisKey edgeKey =
1165 locexp[i]->GetBasis(dir)->GetBasisKey();
1166
1167 LibUtilities::BasisKey nEdgeKey(edgeKey0.GetBasisType(),
1168 edgeKey.GetNumModes(),
1169 edgeKey.GetPointsKey());
1170
1171 exp =
1173 nEdgeKey, segGeom);
1174 }
1175 else
1176 {
1177 exp =
1179 locexp[i]->GetBasis(dir)->GetBasisKey(), segGeom);
1180 }
1181
1182 exp->SetElmtId(elmtid++);
1183 (*m_exp).push_back(exp);
1184 }
1185 }
1186 else if ((exp3D = dynamic_pointer_cast<LocalRegions::Expansion3D>(
1187 locexp[i])))
1188 {
1189 m_expType = e2D;
1190
1191 LibUtilities::BasisKey face0_dir0 =
1192 locexp[i]->GetBasis(0)->GetBasisKey();
1193 LibUtilities::BasisKey face0_dir1 =
1194 locexp[i]->GetBasis(1)->GetBasisKey();
1195
1196 for (j = 0; j < exp3D->GetNtraces(); ++j)
1197 {
1198 FaceGeom = exp3D->GetGeom3D()->GetFace(j);
1199
1200 int dir0 = exp3D->GetGeom3D()->GetDir(j, 0);
1201 int dir1 = exp3D->GetGeom3D()->GetDir(j, 1);
1202
1203 LibUtilities::BasisKey face_dir0 =
1204 locexp[i]->GetBasis(dir0)->GetBasisKey();
1205 LibUtilities::BasisKey face_dir1 =
1206 locexp[i]->GetBasis(dir1)->GetBasisKey();
1207
1208 if ((QuadGeom =
1209 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
1210 FaceGeom)))
1211 {
1212 exp =
1214 face_dir0, face_dir1, QuadGeom);
1215 }
1216 else if ((TriGeom = std::dynamic_pointer_cast<
1217 SpatialDomains::TriGeom>(FaceGeom)))
1218 {
1219 LibUtilities::BasisKey nface_dir0(face0_dir0.GetBasisType(),
1220 face_dir0.GetNumModes(),
1221 face_dir0.GetPointsKey());
1222 LibUtilities::BasisKey nface_dir1(face0_dir1.GetBasisType(),
1223 face_dir1.GetNumModes(),
1224 face_dir1.GetPointsKey());
1225 exp =
1227 nface_dir0, nface_dir1, TriGeom);
1228 }
1229 exp->SetElmtId(elmtid++);
1230 (*m_exp).push_back(exp);
1231 }
1232 }
1233 }
1234
1235 // Set up m_coeffs, m_phys and offset arrays.
1236 SetupCoeffPhys(DeclareCoeffPhysArrays);
1237
1238 // Set up collections
1239 if (m_expType != e0D)
1240 {
1241 CreateCollections(ImpType);
1242 }
1243}
1244
1245/**
1246 * Fills the list of local expansions with the trace from the mesh
1247 * specified by \a domain. This CompositeMap contains a list of
1248 * Composites which define the boundary. It is also used to set up
1249 * expansion domains in the 1D Pulse Wave solver.
1250 *
1251 * @param pSession A session within information about expansion
1252 * @param domain A domain, comprising of one or more composite
1253 * regions,
1254 * @param graph A mesh, containing information about the
1255 * domain and the spectral/hp element expansion.
1256 * @param DeclareCoeffPhysArrays Declare the coefficient and
1257 * phys space arrays. Default is true.
1258 * @param variable The variable name associated with the expansion
1259 * @param SetToOneSpaceDimension Reduce to one space dimension expansion
1260 * @param comm An optional communicator that can be used with the
1261 * boundary expansion in case of more global
1262 * parallel operations. Default to a Null Communicator
1263 * @param ImpType Detail about the implementation type to use
1264 * in operators. Default is eNoImpType.
1265 *
1266 * By default the routine will try and order the expansions in a
1267 * manner which is optimal for collection type operations. This can be
1268 * disabled by the command line option --no-exp-opt
1269 */
1271 const SpatialDomains::CompositeMap &domain,
1273 const bool DeclareCoeffPhysArrays, const std::string variable,
1274 bool SetToOneSpaceDimension,
1275 const LibUtilities::CommSharedPtr comm,
1276 const Collections::ImplementationType ImpType)
1277 : m_comm(comm), m_session(pSession), m_graph(graph), m_physState(false),
1278 m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
1279 m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
1280 m_WaveSpace(false)
1281{
1282 int j, elmtid = 0;
1287
1289
1291
1292 int meshdim = graph->GetMeshDimension();
1293
1294 // Retrieve the list of expansions (element exp)
1295 const SpatialDomains::ExpansionInfoMap &expansions =
1296 graph->GetExpansionInfo(variable);
1297 map<int, vector<SpatialDomains::ExpansionInfoShPtr>> ExpOrder;
1299
1300 bool UseGLLOnTri = false;
1301 // use GLL in all directions on Triangles if Continuous Expansion
1302 pSession->MatchSolverInfo("Projection", "Continuous", UseGLLOnTri, false);
1303
1304 bool DoOptOnCollection =
1305 m_session->DefinesCmdLineArgument("no-exp-opt") ? false : true;
1306 int cnt = 0;
1307 for (auto &compIt : domain)
1308 {
1309 bool IsNot0D = true; // Cehck for 0D expansion
1310 // Process each expansion in the region.
1311 for (j = 0; j < compIt.second->m_geomVec.size(); ++j)
1312 {
1316 compIt.second->m_geomVec[j], PtBvec);
1317
1318 if ((SegGeom = std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1319 compIt.second->m_geomVec[j])))
1320 {
1321 if (meshdim == 1)
1322 {
1323 auto expInfo = expansions.find(SegGeom->GetGlobalID());
1324 eInfo = expInfo->second;
1325 }
1326 else // get bkey from Tri or Quad
1327 {
1328 // First, create the element stdExp that the edge belongs to
1330 graph->GetElementsFromEdge(SegGeom);
1331 // elmts -> std::vector<std::pair<GeometrySharedPtr, int> >
1332 // Currently we assume the elements adjacent to the edge
1333 // have the same type. So we directly fetch the first
1334 // element.
1335 SpatialDomains::GeometrySharedPtr geom = elmts->at(0).first;
1336 int edge_id = elmts->at(0).second;
1338 graph->GetExpansionInfo(geom, variable);
1339 LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1340 LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1342
1343 if (geom->GetShapeType() == LibUtilities::eTriangle)
1344 {
1345 elmtStdExp = MemoryManager<
1346 StdRegions::StdTriExp>::AllocateSharedPtr(Ba, Bb);
1347 }
1348 else if (geom->GetShapeType() ==
1350 {
1351 elmtStdExp = MemoryManager<
1352 StdRegions::StdQuadExp>::AllocateSharedPtr(Ba, Bb);
1353 }
1354 else
1355 {
1357 "Fail to cast geom to a known 2D shape.");
1358 }
1359 // Then, get the trace basis key from the element stdExp,
1360 // which may be different from Ba and Bb.
1361 eInfo->m_basisKeyVector.push_back(
1362 elmtStdExp->GetTraceBasisKey(edge_id));
1363 }
1364 }
1365 else if ((TriGeom =
1366 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
1367 compIt.second->m_geomVec[j])))
1368 {
1369 // First, create the element stdExp that the face belongs to
1371 graph->GetElementsFromFace(TriGeom);
1372 // elmts -> std::vector<std::pair<GeometrySharedPtr, int> >
1373 // Currently we assume the elements adjacent to the face have
1374 // the same type. So we directly fetch the first element.
1375 SpatialDomains::GeometrySharedPtr geom = elmts->at(0).first;
1376 int face_id = elmts->at(0).second;
1377 auto expInfo = expansions.find(geom->GetGlobalID());
1378 ASSERTL0(expInfo != expansions.end(),
1379 "Failed to find expansion info");
1381 expInfo->second->m_basisKeyVector[0];
1383 expInfo->second->m_basisKeyVector[1];
1385 expInfo->second->m_basisKeyVector[2];
1387
1388 if (geom->GetShapeType() == LibUtilities::ePrism)
1389 {
1390 elmtStdExp = MemoryManager<
1391 StdRegions::StdPrismExp>::AllocateSharedPtr(Ba, Bb, Bc);
1392 }
1393 else if (geom->GetShapeType() == LibUtilities::eTetrahedron)
1394 {
1395 elmtStdExp =
1397 Ba, Bb, Bc);
1398 }
1399 else if (geom->GetShapeType() == LibUtilities::ePyramid)
1400 {
1401 elmtStdExp =
1403 Ba, Bb, Bc);
1404 }
1405 else // hex cannot have tri surface
1406 {
1408 "Fail to cast geom to a known 3D shape.");
1409 }
1410 // Then, get the trace basis key from the element stdExp,
1411 // which may be different from Ba, Bb and Bc.
1413 elmtStdExp->GetTraceBasisKey(face_id, 0, UseGLLOnTri);
1415 elmtStdExp->GetTraceBasisKey(face_id, 1, UseGLLOnTri);
1416 // swap TriBa and TriBb orientation is transposed
1417 if (geom->GetForient(face_id) >= 9)
1418 {
1419 std::swap(TriBa, TriBb);
1420 }
1421
1422 eInfo->m_basisKeyVector.push_back(TriBa);
1423 eInfo->m_basisKeyVector.push_back(TriBb);
1424 }
1425 else if ((QuadGeom =
1426 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
1427 compIt.second->m_geomVec[j])))
1428 {
1429 // First, create the element stdExp that the face belongs to
1431 graph->GetElementsFromFace(QuadGeom);
1432 // elmts -> std::vector<std::pair<GeometrySharedPtr, int> >
1433 // Currently we assume the elements adjacent to the face have
1434 // the same type. So we directly fetch the first element.
1435 SpatialDomains::GeometrySharedPtr geom = elmts->at(0).first;
1436 int face_id = elmts->at(0).second;
1437 auto expInfo = expansions.find(geom->GetGlobalID());
1438 ASSERTL0(expInfo != expansions.end(),
1439 "Failed to find expansion info");
1441 expInfo->second->m_basisKeyVector[0];
1443 expInfo->second->m_basisKeyVector[1];
1445 expInfo->second->m_basisKeyVector[2];
1447
1448 if (geom->GetShapeType() == LibUtilities::ePrism)
1449 {
1450 elmtStdExp = MemoryManager<
1451 StdRegions::StdPrismExp>::AllocateSharedPtr(Ba, Bb, Bc);
1452 }
1453 else if (geom->GetShapeType() == LibUtilities::eHexahedron)
1454 {
1455 elmtStdExp =
1457 Ba, Bb, Bc);
1458 }
1459 else if (geom->GetShapeType() == LibUtilities::ePyramid)
1460 {
1461 elmtStdExp =
1463 Ba, Bb, Bc);
1464 }
1465 else // Tet cannot have quad surface
1466 {
1468 "Fail to cast geom to a known 3D shape.");
1469 }
1470 // Then, get the trace basis key from the element stdExp,
1471 // which may be different from Ba, Bb and Bc.
1472 LibUtilities::BasisKey QuadBa =
1473 elmtStdExp->GetTraceBasisKey(face_id, 0);
1474 LibUtilities::BasisKey QuadBb =
1475 elmtStdExp->GetTraceBasisKey(face_id, 1);
1476 // swap Ba and Bb if the orientation is transposed
1477 if (geom->GetForient(face_id) >= 9)
1478 {
1479 std::swap(QuadBa, QuadBb);
1480 }
1481
1482 eInfo->m_basisKeyVector.push_back(QuadBa);
1483 eInfo->m_basisKeyVector.push_back(QuadBb);
1484 }
1485 else
1486 {
1487 IsNot0D = false;
1488 }
1489
1490 if (DoOptOnCollection && IsNot0D)
1491 {
1492 int i;
1493 for (i = 0; i < cnt; ++i)
1494 {
1495 if ((eInfo->m_basisKeyVector ==
1496 ExpOrder[i][0]->m_basisKeyVector) &&
1497 (eInfo->m_geomShPtr->GetGeomFactors()->GetGtype() ==
1498 ExpOrder[i][0]
1499 ->m_geomShPtr->GetGeomFactors()
1500 ->GetGtype()))
1501 {
1502 ExpOrder[i].push_back(eInfo);
1503 break;
1504 }
1505 }
1506
1507 if (i == cnt)
1508 {
1509 ExpOrder[cnt++].push_back(eInfo);
1510 }
1511 }
1512 else
1513 {
1514 ExpOrder[0].push_back(eInfo);
1515 }
1516 }
1517 }
1518
1519 // decalare expansions in provided order using geom and basis info
1520 for (auto &ordIt : ExpOrder)
1521 {
1522 for (auto &eit : ordIt.second)
1523 {
1524 // Process each expansion in the region.
1525 if ((PtGeom = std::dynamic_pointer_cast<SpatialDomains::PointGeom>(
1526 eit->m_geomShPtr)))
1527 {
1528 m_expType = e0D;
1529
1531 PtGeom);
1532 }
1533 else if ((SegGeom =
1534 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1535 eit->m_geomShPtr)))
1536 {
1537 m_expType = e1D;
1538
1539 if (SetToOneSpaceDimension)
1540 {
1541 SpatialDomains::SegGeomSharedPtr OneDSegmentGeom =
1542 SegGeom->GenerateOneSpaceDimGeom();
1543
1544 exp =
1546 eit->m_basisKeyVector[0], OneDSegmentGeom);
1547 }
1548 else
1549 {
1550 exp =
1552 eit->m_basisKeyVector[0], SegGeom);
1553 }
1554 }
1555 else if ((TriGeom =
1556 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
1557 eit->m_geomShPtr)))
1558 {
1559 m_expType = e2D;
1560
1561 if (eit->m_basisKeyVector[0].GetBasisType() ==
1563 {
1565
1567 AllocateSharedPtr(eit->m_basisKeyVector[0],
1568 eit->m_basisKeyVector[1], TriNb,
1569 TriGeom);
1570 }
1571 else
1572 {
1573 exp =
1575 eit->m_basisKeyVector[0], eit->m_basisKeyVector[1],
1576 TriGeom);
1577 }
1578 }
1579 else if ((QuadGeom =
1580 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
1581 eit->m_geomShPtr)))
1582 {
1583 m_expType = e2D;
1584
1586 eit->m_basisKeyVector[0], eit->m_basisKeyVector[1],
1587 QuadGeom);
1588 }
1589 else
1590 {
1592 "dynamic cast to a Geom (possibly 3D) failed");
1593 }
1594
1595 exp->SetElmtId(elmtid++);
1596 (*m_exp).push_back(exp);
1597 }
1598 }
1599
1600 // Set up m_coeffs, m_phys and offset arrays.
1601 SetupCoeffPhys(DeclareCoeffPhysArrays);
1602
1603 if (m_expType != e0D)
1604 {
1605 CreateCollections(ImpType);
1606 }
1607}
1608
1609/**
1610 * Each expansion (local element) is processed in turn to
1611 * determine the number of coefficients and physical data
1612 * points it contributes to the domain. Two arrays,
1613 * #m_coeff_offset are #m_phys_offset are also initialised and
1614 * updated to store the data offsets of each element in the
1615 * #m_coeffs and #m_phys arrays, and the element id that each
1616 * consecutive block is associated respectively.
1617 * Finally we initialise #m_coeffs and #m_phys
1618 */
1619void ExpList::SetupCoeffPhys(bool DeclareCoeffPhysArrays, bool SetupOffsets)
1620{
1621 if (SetupOffsets)
1622 {
1623 int i;
1624
1625 // Set up offset information and array sizes
1628
1629 m_ncoeffs = m_npoints = 0;
1630
1631 for (i = 0; i < m_exp->size(); ++i)
1632 {
1635 m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
1636 m_npoints += (*m_exp)[i]->GetTotPoints();
1637 }
1638 }
1639
1640 if (DeclareCoeffPhysArrays)
1641 {
1644 }
1645
1647
1648 for (int i = 0; i < m_exp->size(); ++i)
1649 {
1650 int coeffs_offset = m_coeff_offset[i];
1651
1652 int loccoeffs = (*m_exp)[i]->GetNcoeffs();
1653
1654 for (int j = 0; j < loccoeffs; ++j)
1655 {
1656 m_coeffsToElmt[coeffs_offset + j].first = i;
1657 m_coeffsToElmt[coeffs_offset + j].second = j;
1658 }
1659 }
1660}
1661
1662/**
1663 * Initialise an expansion vector (m_exp) given an expansion map
1664 * expmap which contains a list of basiskeys and geometries pointers.
1665 * This routine is called from a number of ExpList constructors mainly
1666 * handling the domain definitions. Boundary condition expansions are
1667 * handled with specialised operators.
1668 *
1669 * @param expmap The expansion info map contaiining map of basiskeys
1670 * and geometry pointers
1671 *
1672 * By default the routine will try and order the expansions in a
1673 * manner which is optimal for collection type operations. This can be
1674 * disabled by the command line option --no-exp-opt
1675 */
1678{
1681 SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
1686
1687 int id = 0;
1689
1690 bool DoOptOnCollection =
1691 m_session->DefinesCmdLineArgument("no-exp-opt") ? false : true;
1692 map<int, vector<int>> ExpOrder;
1693 if (DoOptOnCollection)
1694 {
1695 auto expIt = expmap.begin();
1696 int cnt = 0;
1697
1698 ExpOrder[cnt++].push_back(expIt->first);
1699 expIt++;
1700
1701 // sort base on basis key and deformed or regular
1702 for (; expIt != expmap.end(); ++expIt)
1703 {
1704 int i;
1705 for (i = 0; i < cnt; ++i)
1706 {
1708 expmap.find(ExpOrder[i][0])->second;
1709
1710 if ((expIt->second->m_basisKeyVector ==
1711 expInfo->m_basisKeyVector) &&
1712 (expIt->second->m_geomShPtr->GetGeomFactors()->GetGtype() ==
1713 expInfo->m_geomShPtr->GetGeomFactors()->GetGtype()))
1714 {
1715 ExpOrder[i].push_back(expIt->first);
1716 break;
1717 }
1718 }
1719
1720 if (i == cnt) // new expansion
1721 {
1722 ExpOrder[cnt++].push_back(expIt->first);
1723 }
1724 }
1725 }
1726 else
1727 {
1728 for (auto &expIt : expmap) // process in order or global id
1729 {
1730 ExpOrder[0].push_back(expIt.first);
1731 }
1732 }
1733
1735
1736 // Process each expansion in the graph
1737 for (auto &it : ExpOrder)
1738 {
1739 for (int c = 0; c < it.second.size(); ++c)
1740 {
1741 auto expIt = expmap.find(it.second[c]);
1742
1743 const SpatialDomains::ExpansionInfoShPtr expInfo = expIt->second;
1744
1745 switch (expInfo->m_basisKeyVector.size())
1746 {
1747 case 1: // Segment Expansions
1748 {
1750 "Cannot mix expansion dimensions in one vector");
1751 m_expType = e1D;
1752
1753 if ((SegmentGeom =
1754 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1755 expInfo->m_geomShPtr)))
1756 {
1757 // Retrieve basis key from expansion
1759 expInfo->m_basisKeyVector[0];
1760
1762 AllocateSharedPtr(bkey, SegmentGeom);
1763 }
1764 else
1765 {
1767 "dynamic cast to a 1D Geom failed");
1768 }
1769 }
1770 break;
1771 case 2:
1772 {
1774 "Cannot mix expansion dimensions in one vector");
1775 m_expType = e2D;
1776
1777 LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1778 LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1779
1780 if ((TriangleGeom = std::dynamic_pointer_cast<
1781 SpatialDomains ::TriGeom>(expInfo->m_geomShPtr)))
1782 {
1783 // This is not elegantly implemented needs re-thinking.
1785 {
1787 Ba.GetNumModes(),
1788 Ba.GetPointsKey());
1789
1793 AllocateSharedPtr(newBa, Bb, TriNb,
1794 TriangleGeom);
1795 }
1796 else
1797 {
1799 AllocateSharedPtr(Ba, Bb, TriangleGeom);
1800 }
1801 }
1802 else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
1804 expInfo->m_geomShPtr)))
1805 {
1807 AllocateSharedPtr(Ba, Bb, QuadrilateralGeom);
1808 }
1809 else
1810 {
1812 "dynamic cast to a 2D Geom failed");
1813 }
1814 }
1815 break;
1816 case 3:
1817 {
1819 "Cannot mix expansion dimensions in one vector");
1820 m_expType = e3D;
1821
1822 LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1823 LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1824 LibUtilities::BasisKey Bc = expInfo->m_basisKeyVector[2];
1825
1826 if ((TetGeom =
1827 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
1828 expInfo->m_geomShPtr)))
1829 {
1832 {
1833 NEKERROR(
1835 "LocalRegions::NodalTetExp is not implemented "
1836 "yet");
1837 }
1838 else
1839 {
1841 AllocateSharedPtr(Ba, Bb, Bc, TetGeom);
1842 }
1843 }
1844 else if ((PrismGeom = std::dynamic_pointer_cast<
1845 SpatialDomains ::PrismGeom>(
1846 expInfo->m_geomShPtr)))
1847 {
1849 AllocateSharedPtr(Ba, Bb, Bc, PrismGeom);
1850 }
1851 else if ((PyrGeom = std::dynamic_pointer_cast<
1853 expInfo->m_geomShPtr)))
1854 {
1855 exp = MemoryManager<
1856 LocalRegions::PyrExp>::AllocateSharedPtr(Ba, Bb, Bc,
1857 PyrGeom);
1858 }
1859 else if ((HexGeom = std::dynamic_pointer_cast<
1861 expInfo->m_geomShPtr)))
1862 {
1863 exp = MemoryManager<
1864 LocalRegions::HexExp>::AllocateSharedPtr(Ba, Bb, Bc,
1865 HexGeom);
1866 }
1867 else
1868 {
1870 "dynamic cast to a Geom failed");
1871 }
1872 }
1873 break;
1874 default:
1876 "Dimension of basis key is greater than 3");
1877 }
1878
1879 // Assign next id
1880 m_elmtToExpId[exp->GetGeom()->GetGlobalID()] = id;
1881 exp->SetElmtId(id++);
1882
1883 // Add the expansion
1884 (*m_exp).push_back(exp);
1885 }
1886 }
1887}
1888
1889/**
1890 *
1891 */
1893{
1894 return m_expType;
1895}
1896
1898{
1899}
1900
1901/**
1902 * Retrieves the block matrix specified by \a bkey, and computes
1903 * \f$ y=Mx \f$.
1904 * @param gkey GlobalMatrixKey specifying the block matrix to
1905 * use in the matrix-vector multiply.
1906 * @param inarray Input vector \f$ x \f$.
1907 * @param outarray Output vector \f$ y \f$.
1908 */
1910 const Array<OneD, const NekDouble> &inarray,
1911 Array<OneD, NekDouble> &outarray)
1912{
1913 // Retrieve the block matrix using the given key.
1914 const DNekScalBlkMatSharedPtr &blockmat = GetBlockMatrix(gkey);
1915 int nrows = blockmat->GetRows();
1916 int ncols = blockmat->GetColumns();
1917
1918 // Create NekVectors from the given data arrays
1919 NekVector<NekDouble> in(ncols, inarray, eWrapper);
1920 NekVector<NekDouble> out(nrows, outarray, eWrapper);
1921
1922 // Perform matrix-vector multiply.
1923 out = (*blockmat) * in;
1924}
1925
1926/**
1927 * multiply the metric jacobi and quadrature weights
1928 */
1930 const Array<OneD, const NekDouble> &inarray,
1931 Array<OneD, NekDouble> &outarray)
1932{
1933 Array<OneD, NekDouble> e_outarray;
1934
1935 for (int i = 0; i < (*m_exp).size(); ++i)
1936 {
1937 (*m_exp)[i]->MultiplyByQuadratureMetric(inarray + m_phys_offset[i],
1938 e_outarray = outarray +
1939 m_phys_offset[i]);
1940 }
1941}
1942
1943/**
1944 * Divided by the metric jacobi and quadrature weights
1945 */
1947 const Array<OneD, const NekDouble> &inarray,
1948 Array<OneD, NekDouble> &outarray)
1949{
1950 Array<OneD, NekDouble> e_outarray;
1951
1952 for (int i = 0; i < (*m_exp).size(); ++i)
1953 {
1954 (*m_exp)[i]->DivideByQuadratureMetric(inarray + m_phys_offset[i],
1955 e_outarray =
1956 outarray + m_phys_offset[i]);
1957 }
1958}
1959
1960/**
1961 * The operation is evaluated locally for every element by the function
1962 * StdRegions#Expansion#IProductWRTDerivBase.
1963 *
1964 * @param dir {0,1} is the direction in which the
1965 * derivative of the basis should be taken
1966 * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1967 * containing the values of the function
1968 * \f$f(\boldsymbol{x})\f$ at the quadrature
1969 * points \f$\boldsymbol{x}_i\f$.
1970 * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
1971 * used to store the result.
1972 */
1974 const int dir, const Array<OneD, const NekDouble> &inarray,
1975 Array<OneD, NekDouble> &outarray)
1976{
1977 int i;
1978
1979 Array<OneD, NekDouble> e_outarray;
1980 for (i = 0; i < (*m_exp).size(); ++i)
1981 {
1982 (*m_exp)[i]->IProductWRTDerivBase(dir, inarray + m_phys_offset[i],
1983 e_outarray =
1984 outarray + m_coeff_offset[i]);
1985 }
1986}
1987
1988/**
1989 * @brief Directional derivative along a given direction
1990 *
1991 */
1993 const Array<OneD, const NekDouble> &direction,
1994 const Array<OneD, const NekDouble> &inarray,
1995 Array<OneD, NekDouble> &outarray)
1996{
1997 int npts_e;
1998 int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1999 int nq = direction.size() / coordim;
2000
2001 Array<OneD, NekDouble> e_outarray;
2002 Array<OneD, NekDouble> e_MFdiv;
2003
2005
2006 for (int i = 0; i < (*m_exp).size(); ++i)
2007 {
2008 npts_e = (*m_exp)[i]->GetTotPoints();
2009 locdir = Array<OneD, NekDouble>(npts_e * coordim);
2010
2011 for (int k = 0; k < coordim; ++k)
2012 {
2013 Vmath::Vcopy(npts_e, &direction[k * nq + m_phys_offset[i]], 1,
2014 &locdir[k * npts_e], 1);
2015 }
2016
2017 (*m_exp)[i]->IProductWRTDirectionalDerivBase(
2018 locdir, inarray + m_phys_offset[i],
2019 e_outarray = outarray + m_coeff_offset[i]);
2020 }
2021}
2022
2023/**
2024 * The operation is evaluated locally for every element by the function
2025 * StdRegions#StdExpansion#IProductWRTDerivBase.
2026 *
2027 * @param inarray An array of arrays of size \f$Q_{\mathrm{tot}}\f$
2028 * containing the values of the function
2029 * \f$f(\boldsymbol{x})\f$ at the quadrature
2030 * points \f$\boldsymbol{x}_i\f$ in dir directions.
2031 * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
2032 * used to store the result.
2033 */
2035 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
2036 Array<OneD, NekDouble> &outarray)
2037{
2038 Array<OneD, NekDouble> tmp0, tmp1, tmp2;
2039 // assume coord dimension defines the size of Deriv Base
2040 int dim = GetCoordim(0);
2041
2042 ASSERTL1(inarray.size() >= dim, "inarray is not of sufficient dimension");
2043
2044 // initialise if required
2046 {
2047 for (int i = 0; i < m_collections.size(); ++i)
2048 {
2050 }
2052 }
2053
2054 LibUtilities::Timer timer;
2055 int input_offset{0};
2056 int output_offset{0};
2057 LIKWID_MARKER_START("IProductWRTDerivBase_coll");
2058 timer.Start();
2059
2060 switch (dim)
2061 {
2062 case 1:
2063 for (int i = 0; i < m_collections.size(); ++i)
2064 {
2065 m_collections[i].ApplyOperator(
2067 inarray[0] + input_offset, tmp0 = outarray + output_offset);
2068 input_offset += m_collections[i].GetInputSize(
2070 output_offset += m_collections[i].GetOutputSize(
2072 }
2073 break;
2074 case 2:
2075 for (int i = 0; i < m_collections.size(); ++i)
2076 {
2077 m_collections[i].ApplyOperator(
2079 inarray[0] + input_offset, tmp0 = inarray[1] + input_offset,
2080 tmp1 = outarray + output_offset);
2081 input_offset += m_collections[i].GetInputSize(
2083 output_offset += m_collections[i].GetOutputSize(
2085 }
2086 break;
2087 case 3:
2088 for (int i = 0; i < m_collections.size(); ++i)
2089 {
2090 m_collections[i].ApplyOperator(
2092 inarray[0] + input_offset, tmp0 = inarray[1] + input_offset,
2093 tmp1 = inarray[2] + input_offset,
2094 tmp2 = outarray + output_offset);
2095 input_offset += m_collections[i].GetInputSize(
2097 output_offset += m_collections[i].GetOutputSize(
2099 }
2100 break;
2101 default:
2102 NEKERROR(ErrorUtil::efatal, "Dimension of inarray not correct");
2103 break;
2104 }
2105
2106 timer.Stop();
2107 LIKWID_MARKER_STOP("IProductWRTDerivBase_coll");
2108
2109 // Elapsed time
2110 timer.AccumulateRegion("Collections:IProductWRTDerivBase", 10);
2111}
2112/**
2113 * Given a function \f$f(\boldsymbol{x})\f$ evaluated at
2114 * the quadrature points, this function calculates the
2115 * derivatives \f$\frac{d}{dx_1}\f$, \f$\frac{d}{dx_2}\f$
2116 * and \f$\frac{d}{dx_3}\f$ of the function
2117 * \f$f(\boldsymbol{x})\f$ at the same quadrature
2118 * points. The local distribution of the quadrature points
2119 * allows an elemental evaluation of the derivative. This
2120 * is done by a call to the function
2121 * StdRegions#StdExpansion#PhysDeriv.
2122 *
2123 * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
2124 * containing the values of the function
2125 * \f$f(\boldsymbol{x})\f$ at the quadrature
2126 * points \f$\boldsymbol{x}_i\f$.
2127 * @param out_d0 The discrete evaluation of the
2128 * derivative\f$\frac{d}{dx_1}\f$ will
2129 * be stored in this array of size
2130 * \f$Q_{\mathrm{tot}}\f$.
2131 * @param out_d1 The discrete evaluation of the
2132 * derivative\f$\frac{d}{dx_2}\f$ will be
2133 * stored in this array of size
2134 * \f$Q_{\mathrm{tot}}\f$. Note that if no
2135 * memory is allocated for \a out_d1, the
2136 * derivative \f$\frac{d}{dx_2}\f$ will not be
2137 * calculated.
2138 * @param out_d2 The discrete evaluation of the
2139 * derivative\f$\frac{d}{dx_3}\f$ will be
2140 * stored in this array of size
2141 * \f$Q_{\mathrm{tot}}\f$. Note that if no
2142 * memory is allocated for \a out_d2, the
2143 * derivative \f$\frac{d}{dx_3}\f$ will not be
2144 * calculated.
2145 */
2147 Array<OneD, NekDouble> &out_d0,
2148 Array<OneD, NekDouble> &out_d1,
2149 Array<OneD, NekDouble> &out_d2)
2150{
2151 Array<OneD, NekDouble> e_out_d0;
2152 Array<OneD, NekDouble> e_out_d1;
2153 Array<OneD, NekDouble> e_out_d2;
2154
2155 // initialise if required
2157 {
2158 for (int i = 0; i < m_collections.size(); ++i)
2159 {
2161 }
2163 }
2164
2165 int offset{0};
2166 LibUtilities::Timer timer;
2167 timer.Start();
2168 for (int i = 0; i < m_collections.size(); ++i)
2169 {
2170 e_out_d0 = out_d0 + offset;
2171 e_out_d1 = out_d1 + offset;
2172 e_out_d2 = out_d2 + offset;
2173 m_collections[i].ApplyOperator(Collections::ePhysDeriv,
2174 inarray + offset, e_out_d0, e_out_d1,
2175 e_out_d2);
2176 offset += m_collections[i].GetInputSize(Collections::ePhysDeriv);
2177 }
2178 timer.Stop();
2179 // Elapsed time
2180 timer.AccumulateRegion("Collections:PhysDeriv", 10);
2181}
2182
2183void ExpList::v_PhysDeriv(const int dir,
2184 const Array<OneD, const NekDouble> &inarray,
2186{
2187 Direction edir = DirCartesianMap[dir];
2188 v_PhysDeriv(edir, inarray, out_d);
2189}
2190
2192 const Array<OneD, const NekDouble> &inarray,
2194{
2195 int i;
2196 if (edir == MultiRegions::eS)
2197 {
2198 Array<OneD, NekDouble> e_out_ds;
2199 for (i = 0; i < (*m_exp).size(); ++i)
2200 {
2201 e_out_ds = out_d + m_phys_offset[i];
2202 (*m_exp)[i]->PhysDeriv_s(inarray + m_phys_offset[i], e_out_ds);
2203 }
2204 }
2205 else if (edir == MultiRegions::eN)
2206 {
2207 Array<OneD, NekDouble> e_out_dn;
2208 for (i = 0; i < (*m_exp).size(); i++)
2209 {
2210 e_out_dn = out_d + m_phys_offset[i];
2211 (*m_exp)[i]->PhysDeriv_n(inarray + m_phys_offset[i], e_out_dn);
2212 }
2213 }
2214 else
2215 {
2216 // initialise if required
2218 {
2219 for (int i = 0; i < m_collections.size(); ++i)
2220 {
2222 }
2224 }
2225
2226 // convert enum into int
2227 int intdir = (int)edir;
2228 Array<OneD, NekDouble> e_out_d;
2229 int offset{0};
2230 for (int i = 0; i < m_collections.size(); ++i)
2231 {
2232 e_out_d = out_d + offset;
2233 m_collections[i].ApplyOperator(Collections::ePhysDeriv, intdir,
2234 inarray + offset, e_out_d);
2235 offset += m_collections[i].GetInputSize(Collections::ePhysDeriv);
2236 }
2237 }
2238}
2239
2240/* Computes the curl of velocity = \nabla \times u
2241 * if m_expType == 2D, Q = [omg_z, (nothing done)]
2242 * if m_expType == 3D, Q = [omg_x, omg_y, omg_z]
2243 */
2246{
2247 int nq = GetTotPoints();
2251
2252 switch (m_expType)
2253 {
2254 case e2D:
2255 {
2256 PhysDeriv(xDir, Vel[yDir], Vx);
2257 PhysDeriv(yDir, Vel[xDir], Uy);
2258
2259 Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[0], 1);
2260 }
2261 break;
2262
2263 case e3D:
2264 {
2269
2270 PhysDeriv(Vel[xDir], Dummy, Uy, Uz);
2271 PhysDeriv(Vel[yDir], Vx, Dummy, Vz);
2272 PhysDeriv(Vel[zDir], Wx, Wy, Dummy);
2273
2274 Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
2275 Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
2276 Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
2277 }
2278 break;
2279 default:
2280 ASSERTL0(0, "Dimension not supported by ExpList::Curl");
2281 break;
2282 }
2283}
2284
2285/* Computes the curl of vorticity = \nabla \times \nabla \times u
2286 * if m_expType == 2D, Q = [dy omg_z, -dx omg_z, 0]
2287 *
2288 * if m_expType == 3D, Q = [dy omg_z - dz omg_y,
2289 * dz omg_x - dx omg_z,
2290 * dx omg_y - dy omg_x]
2291 *
2292 */
2295{
2296 int nq = GetTotPoints();
2300
2301 bool halfMode = false;
2302 if (GetExpType() == e3DH1D)
2303 {
2304 m_session->MatchSolverInfo("ModeType", "HalfMode", halfMode, false);
2305 }
2306
2307 switch (m_expType)
2308 {
2309 case e2D:
2310 {
2311 PhysDeriv(xDir, Vel[yDir], Vx);
2312 PhysDeriv(yDir, Vel[xDir], Uy);
2313
2314 Vmath::Vsub(nq, Vx, 1, Uy, 1, Dummy, 1);
2315
2316 PhysDeriv(Dummy, Q[1], Q[0]);
2317
2318 Vmath::Smul(nq, -1.0, Q[1], 1, Q[1], 1);
2319 }
2320 break;
2321
2322 case e3D:
2323 case e3DH1D:
2324 case e3DH2D:
2325 {
2330
2331 PhysDeriv(Vel[xDir], Dummy, Uy, Uz);
2332 PhysDeriv(Vel[yDir], Vx, Dummy, Vz);
2333 PhysDeriv(Vel[zDir], Wx, Wy, Dummy);
2334
2335 Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
2336 Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
2337 Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
2338
2339 PhysDeriv(Q[0], Dummy, Uy, Uz);
2340 PhysDeriv(Q[1], Vx, Dummy, Vz);
2341 PhysDeriv(Q[2], Wx, Wy, Dummy);
2342
2343 // For halfmode, need to change the sign of z derivatives
2344 if (halfMode)
2345 {
2346 Vmath::Neg(nq, Uz, 1);
2347 Vmath::Neg(nq, Vz, 1);
2348 }
2349
2350 Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
2351 Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
2352 Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
2353 }
2354 break;
2355 default:
2356 ASSERTL0(0, "Dimension not supported");
2357 break;
2358 }
2359}
2360
2362 const Array<OneD, const NekDouble> &direction,
2363 const Array<OneD, const NekDouble> &inarray,
2364 Array<OneD, NekDouble> &outarray)
2365{
2366 int npts_e;
2367 int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
2368 int nq = direction.size() / coordim;
2369
2370 Array<OneD, NekDouble> e_outarray;
2371 Array<OneD, NekDouble> e_MFdiv;
2373
2374 for (int i = 0; i < (*m_exp).size(); ++i)
2375 {
2376 npts_e = (*m_exp)[i]->GetTotPoints();
2377 locdir = Array<OneD, NekDouble>(npts_e * coordim);
2378
2379 for (int k = 0; k < coordim; ++k)
2380 {
2381 Vmath::Vcopy(npts_e, &direction[k * nq + m_phys_offset[i]], 1,
2382 &locdir[k * npts_e], 1);
2383 }
2384
2385 (*m_exp)[i]->PhysDirectionalDeriv(inarray + m_phys_offset[i], locdir,
2386 e_outarray =
2387 outarray + m_phys_offset[i]);
2388 }
2389}
2390
2392 const NekDouble alpha, const NekDouble exponent,
2393 const NekDouble cutoff)
2394{
2395 Array<OneD, NekDouble> e_array;
2396
2397 for (int i = 0; i < (*m_exp).size(); ++i)
2398 {
2399 (*m_exp)[i]->ExponentialFilter(e_array = array + m_phys_offset[i],
2400 alpha, exponent, cutoff);
2401 }
2402}
2403
2404/**
2405 * The coefficients of the function to be acted upon
2406 * should be contained in the \param inarray. The
2407 * resulting coefficients are stored in \param outarray
2408 *
2409 * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
2410 * containing the inner product.
2411 */
2413 Array<OneD, NekDouble> &outarray)
2414{
2416 const DNekScalBlkMatSharedPtr &InvMass = GetBlockMatrix(mkey);
2417
2418 // Inverse mass matrix
2419 NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
2420 if (inarray.data() == outarray.data())
2421 {
2422 NekVector<NekDouble> in(m_ncoeffs, inarray); // copy data
2423 out = (*InvMass) * in;
2424 }
2425 else
2426 {
2428 out = (*InvMass) * in;
2429 }
2430}
2431
2432/**
2433 * Given a function \f$u(\boldsymbol{x})\f$ defined at the
2434 * quadrature points, this function determines the
2435 * transformed elemental coefficients \f$\hat{u}_n^e\f$
2436 * employing a discrete elemental Galerkin projection from
2437 * physical space to coefficient space. For each element,
2438 * the operation is evaluated locally by the function
2439 * StdRegions#StdExpansion#IproductWRTBase followed by a
2440 * call to #MultiRegions#MultiplyByElmtInvMass.
2441 *
2442 * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
2443 * containing the values of the function
2444 * \f$f(\boldsymbol{x})\f$ at the quadrature
2445 * points \f$\boldsymbol{x}_i\f$.
2446 * @param outarray The resulting coefficients
2447 * \f$\hat{u}_n^e\f$ will be stored in this
2448 * array of size \f$N_{\mathrm{eof}}\f$.
2449 */
2451 Array<OneD, NekDouble> &outarray)
2452{
2454
2455 IProductWRTBase(inarray, f);
2456 MultiplyByElmtInvMass(f, outarray);
2457}
2458
2460 const Array<OneD, const NekDouble> &inarray,
2461 Array<OneD, NekDouble> &outarray)
2462{
2463 int i;
2464
2465 Array<OneD, NekDouble> e_outarray;
2466
2467 for (i = 0; i < (*m_exp).size(); ++i)
2468 {
2469 (*m_exp)[i]->FwdTransBndConstrained(inarray + m_phys_offset[i],
2470 e_outarray =
2471 outarray + m_coeff_offset[i]);
2472 }
2473}
2474
2475/**
2476 * This function smooth a field after some calculaitons which have
2477 * been done elementally.
2478 *
2479 * @param field An array containing the field in physical space
2480 *
2481 */
2483{
2484 // Do nothing unless the method is implemented in the appropriate
2485 // class, i.e. ContField1D,ContField2D, etc.
2486
2487 // So far it has been implemented just for ContField2D and
2488 // ContField3DHomogeneous1D
2489
2490 // Block in case users try the smoothing with a modal expansion.
2491 // Maybe a different techique for the smoothing require
2492 // implementation for modal basis.
2493
2494 ASSERTL0((*m_exp)[0]->GetBasisType(0) == LibUtilities::eGLL_Lagrange ||
2495 (*m_exp)[0]->GetBasisType(0) == LibUtilities::eGauss_Lagrange,
2496 "Smoothing is currently not allowed unless you are using "
2497 "a nodal base for efficiency reasons. The implemented "
2498 "smoothing technique requires the mass matrix inversion "
2499 "which is trivial just for GLL_LAGRANGE_SEM and "
2500 "GAUSS_LAGRANGE_SEMexpansions.");
2501}
2502
2503/**
2504 * This function assembles the block diagonal matrix
2505 * \f$\underline{\boldsymbol{M}}^e\f$, which is the
2506 * concatenation of the local matrices
2507 * \f$\boldsymbol{M}^e\f$ of the type \a mtype, that is
2508 *
2509 * \f[
2510 * \underline{\boldsymbol{M}}^e = \left[
2511 * \begin{array}{cccc}
2512 * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
2513 * 0 & \boldsymbol{M}^2 & 0 & 0 \\
2514 * 0 & 0 & \ddots & 0 \\
2515 * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array}\right].\f]
2516 *
2517 * @param mtype the type of matrix to be assembled
2518 * @param scalar an optional parameter
2519 * @param constant an optional parameter
2520 */
2522 const GlobalMatrixKey &gkey)
2523{
2524 int i, cnt1;
2525 int n_exp = 0;
2526 DNekScalMatSharedPtr loc_mat;
2527 DNekScalBlkMatSharedPtr BlkMatrix;
2528 map<int, int> elmt_id;
2530
2532 {
2533 for (i = 0; i < (*m_exp).size(); ++i)
2534 {
2535 if ((*m_exp)[i]->DetShapeType() == ShapeType)
2536 {
2537 elmt_id[n_exp++] = i;
2538 }
2539 }
2540 }
2541 else
2542 {
2543 n_exp = (*m_exp).size();
2544 for (i = 0; i < n_exp; ++i)
2545 {
2546 elmt_id[i] = i;
2547 }
2548 }
2549
2550 Array<OneD, unsigned int> nrows(n_exp);
2551 Array<OneD, unsigned int> ncols(n_exp);
2552
2553 switch (gkey.GetMatrixType())
2554 {
2556 {
2557 // set up an array of integers for block matrix construction
2558 for (i = 0; i < n_exp; ++i)
2559 {
2560 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
2561 ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
2562 }
2563 }
2564 break;
2566 {
2567 // set up an array of integers for block matrix construction
2568 for (i = 0; i < n_exp; ++i)
2569 {
2570 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
2571 ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
2572 }
2573 }
2574 break;
2575 case StdRegions::eMass:
2580 {
2581 // set up an array of integers for block matrix construction
2582 for (i = 0; i < n_exp; ++i)
2583 {
2584 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
2585 ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
2586 }
2587 }
2588 break;
2589
2591 {
2592 // set up an array of integers for block matrix construction
2593 for (i = 0; i < n_exp; ++i)
2594 {
2595 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
2596 ncols[i] =
2597 (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
2598 }
2599 }
2600 break;
2602 {
2603 // set up an array of integers for block matrix construction
2604 for (i = 0; i < n_exp; ++i)
2605 {
2606 nrows[i] =
2607 (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
2608 ncols[i] =
2609 (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
2610 }
2611 }
2612 break;
2613 default:
2614 {
2616 "Global Matrix creation not defined for this "
2617 "type of matrix");
2618 }
2619 }
2620
2621 MatrixStorage blkmatStorage = eDIAGONAL;
2622 BlkMatrix = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nrows, ncols,
2623 blkmatStorage);
2624
2625 int nvarcoeffs = gkey.GetNVarCoeffs();
2626 int eid;
2627 Array<OneD, NekDouble> varcoeffs_wk;
2628
2629 for (i = cnt1 = 0; i < n_exp; ++i)
2630 {
2631 // need to be initialised with zero size for non
2632 // variable coefficient case
2633 StdRegions::VarCoeffMap varcoeffs;
2634
2635 eid = elmt_id[i];
2636 if (nvarcoeffs > 0)
2637 {
2638 varcoeffs = StdRegions::RestrictCoeffMap(
2639 gkey.GetVarCoeffs(), m_phys_offset[i],
2640 (*m_exp)[i]->GetTotPoints());
2641 }
2642
2644 gkey.GetMatrixType(), (*m_exp)[eid]->DetShapeType(), *(*m_exp)[eid],
2645 gkey.GetConstFactors(), varcoeffs);
2646
2647 loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>(
2648 (*m_exp)[elmt_id.find(i)->second])
2649 ->GetLocMatrix(matkey);
2650 BlkMatrix->SetBlock(i, i, loc_mat);
2651 }
2652
2653 return BlkMatrix;
2654}
2655
2657 const GlobalMatrixKey &gkey)
2658{
2659 auto matrixIter = m_blockMat->find(gkey);
2660
2661 if (matrixIter == m_blockMat->end())
2662 {
2663 return ((*m_blockMat)[gkey] = GenBlockMatrix(gkey));
2664 }
2665 else
2666 {
2667 return matrixIter->second;
2668 }
2669}
2670
2671// Routines for continous matrix solution
2672/**
2673 * This operation is equivalent to the evaluation of
2674 * \f$\underline{\boldsymbol{M}}^e\boldsymbol{\hat{u}}_l\f$, that is,
2675 * \f[ \left[
2676 * \begin{array}{cccc}
2677 * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
2678 * 0 & \boldsymbol{M}^2 & 0 & 0 \\
2679 * 0 & 0 & \ddots & 0 \\
2680 * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array} \right]
2681 *\left [ \begin{array}{c}
2682 * \boldsymbol{\hat{u}}^{1} \\
2683 * \boldsymbol{\hat{u}}^{2} \\
2684 * \vdots \\
2685 * \boldsymbol{\hat{u}}^{{{N_{\mathrm{el}}}}} \end{array} \right ]\f]
2686 * where \f$\boldsymbol{M}^e\f$ are the local matrices of type
2687 * specified by the key \a mkey. The decoupling of the local matrices
2688 * allows for a local evaluation of the operation. However, rather than
2689 * a local matrix-vector multiplication, the local operations are
2690 * evaluated as implemented in the function
2691 * StdRegions#StdExpansion#GeneralMatrixOp.
2692 *
2693 * @param mkey This key uniquely defines the type matrix
2694 * required for the operation.
2695 * @param inarray The vector \f$\boldsymbol{\hat{u}}_l\f$ of
2696 * size \f$N_{\mathrm{eof}}\f$.
2697 * @param outarray The resulting vector of size
2698 * \f$N_{\mathrm{eof}}\f$.
2699 */
2701 const Array<OneD, const NekDouble> &inarray,
2702 Array<OneD, NekDouble> &outarray)
2703{
2704 int nvarcoeffs = gkey.GetNVarCoeffs();
2705
2706 if (gkey.GetMatrixType() == StdRegions::eHelmholtz ||
2708 {
2709 // Map operator type based on matrix type
2714
2715 // initialise if required
2716 if (m_collections.size() && m_collectionsDoInit[opType])
2717 {
2718 for (int i = 0; i < m_collections.size(); ++i)
2719 {
2720 m_collections[i].Initialise(opType, gkey.GetConstFactors());
2721 }
2722 m_collectionsDoInit[opType] = false;
2723 }
2724
2725 // Update factors and varoeffs
2726 int cnt{0};
2727 for (int i = 0; i < m_collections.size(); ++i)
2728 {
2729 m_collections[i].UpdateFactors(opType, gkey.GetConstFactors());
2730
2731 // Restrict varcoeffs to collection size and update
2732 StdRegions::VarCoeffMap varcoeffs;
2733 if (nvarcoeffs)
2734 {
2735 varcoeffs = StdRegions::RestrictCoeffMap(
2736 gkey.GetVarCoeffs(), m_phys_offset[cnt],
2737 m_collections[i].GetInputSize(opType, false));
2738 cnt += m_collections[i].GetNumElmt(opType);
2739 }
2740 m_collections[i].UpdateVarcoeffs(opType, varcoeffs);
2741 }
2742
2744 int input_offset{0};
2745 int output_offset{0};
2746 for (int i = 0; i < m_collections.size(); ++i)
2747 {
2748 // the input_offset is equal to the output_offset - this is
2749 // happenning inside the Helmholtz_Helper or LinearADR_Helper class
2750 m_collections[i].ApplyOperator(opType, inarray + input_offset,
2751 tmp = outarray + output_offset);
2752 input_offset += m_collections[i].GetInputSize(opType);
2753 output_offset += m_collections[i].GetOutputSize(opType);
2754 }
2755 }
2756 else
2757 {
2758 Array<OneD, NekDouble> tmp_outarray;
2759 for (int i = 0; i < (*m_exp).size(); ++i)
2760 {
2761 // need to be initialised with zero size for non
2762 // variable coefficient case
2763 StdRegions::VarCoeffMap varcoeffs;
2764
2765 if (nvarcoeffs > 0)
2766 {
2767 varcoeffs = StdRegions::RestrictCoeffMap(
2768 gkey.GetVarCoeffs(), m_phys_offset[i],
2769 (*m_exp)[i]->GetTotPoints());
2770 }
2771
2773 gkey.GetMatrixType(), (*m_exp)[i]->DetShapeType(),
2774 *((*m_exp)[i]), gkey.GetConstFactors(), varcoeffs);
2775
2776 (*m_exp)[i]->GeneralMatrixOp(
2777 inarray + m_coeff_offset[i],
2778 tmp_outarray = outarray + m_coeff_offset[i], mkey);
2779 }
2780 }
2781}
2782
2783/**
2784 * Retrieves local matrices from each expansion in the expansion list
2785 * and combines them together to generate a global matrix system.
2786 * @param mkey Matrix key for the matrix to be generated.
2787 * @param locToGloMap Local to global mapping.
2788 * @returns Shared pointer to the generated global matrix.
2789 */
2791 const GlobalMatrixKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
2792{
2793 int i, j, n, gid1, gid2, cntdim1, cntdim2;
2794 NekDouble sign1, sign2;
2795 DNekScalMatSharedPtr loc_mat;
2796
2797 unsigned int glob_rows = 0;
2798 unsigned int glob_cols = 0;
2799 unsigned int loc_rows = 0;
2800 unsigned int loc_cols = 0;
2801
2802 bool assembleFirstDim = false;
2803 bool assembleSecondDim = false;
2804
2805 switch (mkey.GetMatrixType())
2806 {
2808 {
2809 glob_rows = m_npoints;
2810 glob_cols = locToGloMap->GetNumGlobalCoeffs();
2811
2812 assembleFirstDim = false;
2813 assembleSecondDim = true;
2814 }
2815 break;
2817 {
2818 glob_rows = locToGloMap->GetNumGlobalCoeffs();
2819 glob_cols = m_npoints;
2820
2821 assembleFirstDim = true;
2822 assembleSecondDim = false;
2823 }
2824 break;
2825 case StdRegions::eMass:
2829 {
2830 glob_rows = locToGloMap->GetNumGlobalCoeffs();
2831 glob_cols = locToGloMap->GetNumGlobalCoeffs();
2832
2833 assembleFirstDim = true;
2834 assembleSecondDim = true;
2835 }
2836 break;
2837 default:
2838 {
2840 "Global Matrix creation not defined for this "
2841 "type of matrix");
2842 }
2843 }
2844
2845 COOMatType spcoomat;
2846 CoordType coord;
2847
2848 int nvarcoeffs = mkey.GetNVarCoeffs();
2849 int eid;
2850
2851 // fill global matrix
2852 for (n = cntdim1 = cntdim2 = 0; n < (*m_exp).size(); ++n)
2853 {
2854 // need to be initialised with zero size for non
2855 // variable coefficient case
2856 StdRegions::VarCoeffMap varcoeffs;
2857
2858 eid = n;
2859 if (nvarcoeffs > 0)
2860 {
2861 varcoeffs = StdRegions::RestrictCoeffMap(
2862 mkey.GetVarCoeffs(), m_phys_offset[eid],
2863 (*m_exp)[eid]->GetTotPoints());
2864 }
2865
2867 mkey.GetMatrixType(), (*m_exp)[eid]->DetShapeType(),
2868 *((*m_exp)[eid]), mkey.GetConstFactors(), varcoeffs);
2869
2870 loc_mat =
2871 std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])
2872 ->GetLocMatrix(matkey);
2873
2874 loc_rows = loc_mat->GetRows();
2875 loc_cols = loc_mat->GetColumns();
2876
2877 for (i = 0; i < loc_rows; ++i)
2878 {
2879 if (assembleFirstDim)
2880 {
2881 gid1 = locToGloMap->GetLocalToGlobalMap(cntdim1 + i);
2882 sign1 = locToGloMap->GetLocalToGlobalSign(cntdim1 + i);
2883 }
2884 else
2885 {
2886 gid1 = cntdim1 + i;
2887 sign1 = 1.0;
2888 }
2889
2890 for (j = 0; j < loc_cols; ++j)
2891 {
2892 if (assembleSecondDim)
2893 {
2894 gid2 = locToGloMap->GetLocalToGlobalMap(cntdim2 + j);
2895 sign2 = locToGloMap->GetLocalToGlobalSign(cntdim2 + j);
2896 }
2897 else
2898 {
2899 gid2 = cntdim2 + j;
2900 sign2 = 1.0;
2901 }
2902
2903 // sparse matrix fill
2904 coord = make_pair(gid1, gid2);
2905 if (spcoomat.count(coord) == 0)
2906 {
2907 spcoomat[coord] = sign1 * sign2 * (*loc_mat)(i, j);
2908 }
2909 else
2910 {
2911 spcoomat[coord] += sign1 * sign2 * (*loc_mat)(i, j);
2912 }
2913 }
2914 }
2915 cntdim1 += loc_rows;
2916 cntdim2 += loc_cols;
2917 }
2918
2920 glob_cols, spcoomat);
2921}
2922
2924 const GlobalLinSysKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
2925{
2926 int i, j, n, gid1, gid2, loc_lda, eid;
2927 NekDouble sign1, sign2, value;
2928 DNekScalMatSharedPtr loc_mat;
2929
2930 int totDofs = locToGloMap->GetNumGlobalCoeffs();
2931 int NumDirBCs = locToGloMap->GetNumGlobalDirBndCoeffs();
2932
2933 unsigned int rows = totDofs - NumDirBCs;
2934 unsigned int cols = totDofs - NumDirBCs;
2935 NekDouble zero = 0.0;
2936
2937 DNekMatSharedPtr Gmat;
2938 int bwidth = locToGloMap->GetFullSystemBandWidth();
2939
2940 int nvarcoeffs = mkey.GetNVarCoeffs();
2941 MatrixStorage matStorage;
2942
2943 map<int, RobinBCInfoSharedPtr> RobinBCInfo = GetRobinBCInfo();
2944
2945 switch (mkey.GetMatrixType())
2946 {
2947 // case for all symmetric matices
2950 if ((2 * (bwidth + 1)) < rows)
2951 {
2954 rows, cols, zero, matStorage, bwidth, bwidth);
2955 }
2956 else
2957 {
2958 matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
2960 rows, cols, zero, matStorage);
2961 }
2962
2963 break;
2964 default: // Assume general matrix - currently only set up
2965 // for full invert
2966 {
2967 matStorage = eFULL;
2969 rows, cols, zero, matStorage);
2970 }
2971 }
2972
2973 // fill global symmetric matrix
2974 for (n = 0; n < (*m_exp).size(); ++n)
2975 {
2976 // need to be initialised with zero size for non
2977 // variable coefficient case
2978 StdRegions::VarCoeffMap varcoeffs;
2979
2980 eid = n;
2981 if (nvarcoeffs > 0)
2982 {
2983 varcoeffs = StdRegions::RestrictCoeffMap(
2984 mkey.GetVarCoeffs(), m_phys_offset[eid],
2985 (*m_exp)[eid]->GetTotPoints());
2986 }
2987
2989 mkey.GetMatrixType(), (*m_exp)[eid]->DetShapeType(),
2990 *((*m_exp)[eid]), mkey.GetConstFactors(), varcoeffs);
2991
2992 loc_mat =
2993 std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])
2994 ->GetLocMatrix(matkey);
2995
2996 if (RobinBCInfo.count(n) != 0) // add robin mass matrix
2997 {
2999
3000 // declare local matrix from scaled matrix.
3001 int rows = loc_mat->GetRows();
3002 int cols = loc_mat->GetColumns();
3003 const NekDouble *dat = loc_mat->GetRawPtr();
3004 DNekMatSharedPtr new_mat =
3006 Blas::Dscal(rows * cols, loc_mat->Scale(), new_mat->GetRawPtr(), 1);
3007
3008 // add local matrix contribution
3009 for (rBC = RobinBCInfo.find(n)->second; rBC; rBC = rBC->next)
3010 {
3011 (*m_exp)[n]->AddRobinMassMatrix(
3012 rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
3013 }
3014
3015 NekDouble one = 1.0;
3016 // redeclare loc_mat to point to new_mat plus the scalar.
3017 loc_mat =
3019 }
3020
3021 loc_lda = loc_mat->GetColumns();
3022
3023 for (i = 0; i < loc_lda; ++i)
3024 {
3025 gid1 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + i) -
3026 NumDirBCs;
3027 sign1 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + i);
3028 if (gid1 >= 0)
3029 {
3030 for (j = 0; j < loc_lda; ++j)
3031 {
3032 gid2 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] +
3033 j) -
3034 NumDirBCs;
3035 sign2 = locToGloMap->GetLocalToGlobalSign(
3036 m_coeff_offset[n] + j);
3037 if (gid2 >= 0)
3038 {
3039 // When global matrix is symmetric,
3040 // only add the value for the upper
3041 // triangular part in order to avoid
3042 // entries to be entered twice
3043 if ((matStorage == eFULL) || (gid2 >= gid1))
3044 {
3045 value = Gmat->GetValue(gid1, gid2) +
3046 sign1 * sign2 * (*loc_mat)(i, j);
3047 Gmat->SetValue(gid1, gid2, value);
3048 }
3049 }
3050 }
3051 }
3052 }
3053 }
3054
3055 return Gmat;
3056}
3057
3058/**
3059 * Consider a linear system
3060 * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{f}\f$ to be solved. Dependent
3061 * on the solution method, this function constructs
3062 * - <b>The full linear system</b><BR>
3063 * A call to the function #GenGlobalLinSysFullDirect
3064 * - <b>The statically condensed linear system</b><BR>
3065 * A call to the function #GenGlobalLinSysStaticCond
3066 *
3067 * @param mkey A key which uniquely defines the global
3068 * matrix to be constructed.
3069 * @param locToGloMap Contains the mapping array and required
3070 * information for the transformation from
3071 * local to global degrees of freedom.
3072 * @return (A shared pointer to) the global linear system in
3073 * required format.
3074 */
3076 const GlobalLinSysKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
3077{
3078 GlobalLinSysSharedPtr returnlinsys;
3079 std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
3080
3082
3083 if (vType >= eSIZE_GlobalSysSolnType)
3084 {
3085 NEKERROR(ErrorUtil::efatal, "Matrix solution type not defined");
3086 }
3087 std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
3088
3089 return GetGlobalLinSysFactory().CreateInstance(vSolnType, mkey, vExpList,
3090 locToGloMap);
3091}
3092
3094 const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
3095{
3096 std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
3097 const map<int, RobinBCInfoSharedPtr> vRobinBCInfo = GetRobinBCInfo();
3098
3100
3101 if (vType >= eSIZE_GlobalSysSolnType)
3102 {
3103 NEKERROR(ErrorUtil::efatal, "Matrix solution type not defined");
3104 }
3105 std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
3106
3107 return GetGlobalLinSysFactory().CreateInstance(vSolnType, mkey, vExpList,
3108 locToGloMap);
3109}
3110
3111/**
3112 * Given the elemental coefficients \f$\hat{u}_n^e\f$ of
3113 * an expansion, this function evaluates the spectral/hp
3114 * expansion \f$u^{\delta}(\boldsymbol{x})\f$ at the
3115 * quadrature points \f$\boldsymbol{x}_i\f$. The operation
3116 * is evaluated locally by the elemental function
3117 * StdRegions#StdExpansion#BwdTrans.
3118 *
3119 * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
3120 * containing the local coefficients
3121 * \f$\hat{u}_n^e\f$.
3122 * @param outarray The resulting physical values at the
3123 * quadrature points
3124 * \f$u^{\delta}(\boldsymbol{x}_i)\f$
3125 * will be stored in this array of size
3126 * \f$Q_{\mathrm{tot}}\f$.
3127 */
3128
3130 Array<OneD, NekDouble> &outarray)
3131{
3132 LibUtilities::Timer timer;
3133
3134 if (m_expType == e0D)
3135 {
3136 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
3137 }
3138 else
3139 {
3140 // initialise if required
3142 {
3143 for (int i = 0; i < m_collections.size(); ++i)
3144 {
3146 }
3148 }
3149
3150 LIKWID_MARKER_START("v_BwdTrans");
3151 timer.Start();
3152
3154 int input_offset{0};
3155 int output_offset{0};
3156 for (int i = 0; i < m_collections.size(); ++i)
3157 {
3158 m_collections[i].ApplyOperator(Collections::eBwdTrans,
3159 inarray + input_offset,
3160 tmp = outarray + output_offset);
3161 input_offset +=
3162 m_collections[i].GetInputSize(Collections::eBwdTrans);
3163 output_offset +=
3164 m_collections[i].GetOutputSize(Collections::eBwdTrans);
3165 }
3166
3167 timer.Stop();
3168 LIKWID_MARKER_STOP("v_BwdTrans");
3169 }
3170 // Elapsed time
3171 timer.AccumulateRegion("Collections:BwdTrans", 10);
3172}
3173
3175 const Array<OneD, const NekDouble> &gloCoord)
3176{
3177 return GetExp(GetExpIndex(gloCoord));
3178}
3179
3180/**
3181 * @todo need a smarter search here that first just looks at bounding
3182 * vertices - suggest first seeing if point is within 10% of
3183 * region defined by vertices. The do point search.
3184 */
3186 NekDouble tol, bool returnNearestElmt, int cachedId,
3187 NekDouble maxDistance)
3188{
3189 Array<OneD, NekDouble> Lcoords(gloCoord.size());
3190
3191 return GetExpIndex(gloCoord, Lcoords, tol, returnNearestElmt, cachedId,
3192 maxDistance);
3193}
3194
3196 Array<OneD, NekDouble> &locCoords, NekDouble tol,
3197 bool returnNearestElmt, int cachedId,
3198 NekDouble maxDistance)
3199{
3200 if (GetNumElmts() == 0)
3201 {
3202 return -1;
3203 }
3204
3205 if (m_elmtToExpId.size() == 0)
3206 {
3207 // Loop in reverse order so that in case where using a
3208 // Homogeneous expansion it sets geometry ids to first part of
3209 // m_exp list. Otherwise will set to second (complex) expansion
3210 for (int i = (*m_exp).size() - 1; i >= 0; --i)
3211 {
3212 m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3213 }
3214 }
3215
3216 NekDouble nearpt = 1e6;
3217 NekDouble nearpt_min = 1e6;
3218 int min_id = -1;
3219 Array<OneD, NekDouble> savLocCoords(locCoords.size());
3220
3221 if (cachedId >= 0 && cachedId < (*m_exp).size())
3222 {
3223 nearpt = 1e12;
3224 if ((*m_exp)[cachedId]->GetGeom()->ContainsPoint(gloCoords, locCoords,
3225 tol, nearpt))
3226 {
3227 return cachedId;
3228 }
3229 else if (returnNearestElmt && (nearpt < nearpt_min))
3230 {
3231 // If it does not lie within, keep track of which element
3232 // is nearest.
3233 min_id = cachedId;
3234 nearpt_min = nearpt;
3235 Vmath::Vcopy(locCoords.size(), locCoords, 1, savLocCoords, 1);
3236 }
3237 }
3238
3239 NekDouble x = (gloCoords.size() > 0 ? gloCoords[0] : 0.0);
3240 NekDouble y = (gloCoords.size() > 1 ? gloCoords[1] : 0.0);
3241 NekDouble z = (gloCoords.size() > 2 ? gloCoords[2] : 0.0);
3244 GetExp(0)->GetCoordim(), -1, x, y, z);
3245
3246 // Get the list of elements whose bounding box contains the desired
3247 // point.
3248 std::vector<int> elmts = m_graph->GetElementsContainingPoint(p);
3249
3250 // Check each element in turn to see if point lies within it.
3251 for (int i = 0; i < elmts.size(); ++i)
3252 {
3253 int id = m_elmtToExpId[elmts[i]];
3254 if (id == cachedId)
3255 {
3256 continue;
3257 }
3258 if ((*m_exp)[id]->GetGeom()->ContainsPoint(gloCoords, locCoords, tol,
3259 nearpt))
3260 {
3261 return id;
3262 }
3263 else if (returnNearestElmt && (nearpt < nearpt_min))
3264 {
3265 // If it does not lie within, keep track of which element
3266 // is nearest.
3267 min_id = id;
3268 nearpt_min = nearpt;
3269 Vmath::Vcopy(locCoords.size(), locCoords, 1, savLocCoords, 1);
3270 }
3271 }
3272
3273 // If the calling function is with just the nearest element, return
3274 // that. Otherwise return -1 to indicate no matching elemenet found.
3275 if (returnNearestElmt && nearpt_min <= maxDistance)
3276 {
3277 Vmath::Vcopy(locCoords.size(), savLocCoords, 1, locCoords, 1);
3278 std::string msg = "Failed to find point within a tolerance of " +
3279 boost::lexical_cast<std::string>(tol) +
3280 ", using local point (";
3281 for (size_t j = 0; j < locCoords.size(); ++j)
3282 {
3283 msg += boost::lexical_cast<std::string>(savLocCoords[j]);
3284 if (j < locCoords.size())
3285 {
3286 msg += ", ";
3287 }
3288 }
3289 msg += ") in element: " + std::to_string(min_id) +
3290 " with a distance of " + std::to_string(nearpt_min);
3291 WARNINGL1(false, msg.c_str());
3292 return min_id;
3293 }
3294 else
3295 {
3296 return -1;
3297 }
3298}
3299
3300/**
3301 * Given some coordinates, output the expansion field value at that
3302 * point
3303 */
3305 const Array<OneD, const NekDouble> &phys)
3306{
3307 int dim = GetCoordim(0);
3308 ASSERTL0(dim == coords.size(), "Invalid coordinate dimension.");
3309
3310 // Grab the element index corresponding to coords.
3311 Array<OneD, NekDouble> xi(dim);
3312 int elmtIdx = GetExpIndex(coords, xi);
3313 ASSERTL0(elmtIdx > 0, "Unable to find element containing point.");
3314
3315 // Grab that element's physical storage.
3316 Array<OneD, NekDouble> elmtPhys = phys + m_phys_offset[elmtIdx];
3317
3318 // Evaluate the element at the appropriate point.
3319 return (*m_exp)[elmtIdx]->StdPhysEvaluate(xi, elmtPhys);
3320}
3321
3322/**
3323 * Configures geometric info, such as tangent direction, on each
3324 * expansion.
3325 * @param graph2D Mesh
3326 */
3328{
3329}
3330
3331/**
3332 * @brief Reset geometry information, metrics, matrix managers and
3333 * geometry information.
3334 *
3335 * This routine clears all matrix managers and resets all geometry
3336 * information, which allows the geometry information to be dynamically
3337 * updated as the solver is run.
3338 */
3340{
3341 // Reset matrix managers.
3343 LocalRegions::MatrixKey::opLess>::ClearManager();
3345 LocalRegions::MatrixKey::opLess>::ClearManager();
3346
3347 // Reset block matrix map
3348 m_blockMat->clear();
3349
3350 // Loop over all elements and reset geometry information.
3351 for (int i = 0; i < m_exp->size(); ++i)
3352 {
3353 (*m_exp)[i]->GetGeom()->Reset(m_graph->GetCurvedEdges(),
3354 m_graph->GetCurvedFaces());
3355 }
3356
3357 // Loop over all elements and rebuild geometric factors.
3358 for (int i = 0; i < m_exp->size(); ++i)
3359 {
3360 (*m_exp)[i]->Reset();
3361 }
3362
3363 CreateCollections(Collections::eNoImpType); // @TODO: Might need to pass in
3364 // correct type here
3365}
3366
3368{
3369 // Reset matrix managers.
3371 LocalRegions::MatrixKey::opLess>::ClearManager();
3373 LocalRegions::MatrixKey::opLess>::ClearManager();
3374
3375 // Reset block matrix map
3376 m_blockMat->clear();
3377}
3378
3379/**
3380 * Write Tecplot Files Header
3381 * @param outfile Output file name.
3382 * @param var variables names
3383 */
3384void ExpList::v_WriteTecplotHeader(std::ostream &outfile, std::string var)
3385{
3386 if (GetNumElmts() == 0)
3387 {
3388 return;
3389 }
3390
3391 int coordim = GetExp(0)->GetCoordim();
3392 char vars[3] = {'x', 'y', 'z'};
3393
3394 if (m_expType == e3DH1D)
3395 {
3396 coordim += 1;
3397 }
3398 else if (m_expType == e3DH2D)
3399 {
3400 coordim += 2;
3401 }
3402
3403 outfile << "Variables = x";
3404 for (int i = 1; i < coordim; ++i)
3405 {
3406 outfile << ", " << vars[i];
3407 }
3408
3409 if (var.size() > 0)
3410 {
3411 outfile << ", " << var;
3412 }
3413
3414 outfile << std::endl << std::endl;
3415}
3416
3417/**
3418 * Write Tecplot Files Zone
3419 * @param outfile Output file name.
3420 * @param expansion Expansion that is considered
3421 */
3422void ExpList::v_WriteTecplotZone(std::ostream &outfile, int expansion)
3423{
3424 int i, j;
3425 int coordim = GetCoordim(0);
3426 int nPoints = GetTotPoints();
3427 int nBases = (*m_exp)[0]->GetNumBases();
3428 int numBlocks = 0;
3429
3431
3432 if (expansion == -1)
3433 {
3434 nPoints = GetTotPoints();
3435
3436 coords[0] = Array<OneD, NekDouble>(nPoints);
3437 coords[1] = Array<OneD, NekDouble>(nPoints);
3438 coords[2] = Array<OneD, NekDouble>(nPoints);
3439
3440 GetCoords(coords[0], coords[1], coords[2]);
3441
3442 for (i = 0; i < m_exp->size(); ++i)
3443 {
3444 int numInt = 1;
3445
3446 for (j = 0; j < nBases; ++j)
3447 {
3448 numInt *= (*m_exp)[i]->GetNumPoints(j) - 1;
3449 }
3450
3451 numBlocks += numInt;
3452 }
3453 }
3454 else
3455 {
3456 nPoints = (*m_exp)[expansion]->GetTotPoints();
3457
3458 coords[0] = Array<OneD, NekDouble>(nPoints);
3459 coords[1] = Array<OneD, NekDouble>(nPoints);
3460 coords[2] = Array<OneD, NekDouble>(nPoints);
3461
3462 (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
3463
3464 numBlocks = 1;
3465 for (j = 0; j < nBases; ++j)
3466 {
3467 numBlocks *= (*m_exp)[expansion]->GetNumPoints(j) - 1;
3468 }
3469 }
3470
3471 if (m_expType == e3DH1D)
3472 {
3473 nBases += 1;
3474 coordim += 1;
3475 int nPlanes = GetZIDs().size();
3476 NekDouble tmp = numBlocks * (nPlanes - 1.0) / nPlanes;
3477 numBlocks = (int)tmp;
3478 }
3479 else if (m_expType == e3DH2D)
3480 {
3481 nBases += 2;
3482 coordim += 1;
3483 }
3484
3485 outfile << "Zone, N=" << nPoints << ", E=" << numBlocks << ", F=FEBlock";
3486
3487 switch (nBases)
3488 {
3489 case 2:
3490 outfile << ", ET=QUADRILATERAL" << std::endl;
3491 break;
3492 case 3:
3493 outfile << ", ET=BRICK" << std::endl;
3494 break;
3495 default:
3496 NEKERROR(ErrorUtil::efatal, "Not set up for this type of output");
3497 break;
3498 }
3499
3500 // Write out coordinates
3501 for (j = 0; j < coordim; ++j)
3502 {
3503 for (i = 0; i < nPoints; ++i)
3504 {
3505 outfile << coords[j][i] << " ";
3506 if (i % 1000 == 0 && i)
3507 {
3508 outfile << std::endl;
3509 }
3510 }
3511 outfile << std::endl;
3512 }
3513}
3514
3515void ExpList::v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
3516{
3517 int i, j, k, l;
3518 int nbase = (*m_exp)[0]->GetNumBases();
3519 int cnt = 0;
3520
3521 std::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
3522
3523 if (expansion != -1)
3524 {
3525 exp = std::shared_ptr<LocalRegions::ExpansionVector>(
3527 (*exp)[0] = (*m_exp)[expansion];
3528 }
3529
3530 if (nbase == 2)
3531 {
3532 for (i = 0; i < (*exp).size(); ++i)
3533 {
3534 const int np0 = (*exp)[i]->GetNumPoints(0);
3535 const int np1 = (*exp)[i]->GetNumPoints(1);
3536
3537 for (j = 1; j < np1; ++j)
3538 {
3539 for (k = 1; k < np0; ++k)
3540 {
3541 outfile << cnt + (j - 1) * np0 + k << " ";
3542 outfile << cnt + (j - 1) * np0 + k + 1 << " ";
3543 outfile << cnt + j * np0 + k + 1 << " ";
3544 outfile << cnt + j * np0 + k << endl;
3545 }
3546 }
3547
3548 cnt += np0 * np1;
3549 }
3550 }
3551 else if (nbase == 3)
3552 {
3553 for (i = 0; i < (*exp).size(); ++i)
3554 {
3555 const int np0 = (*exp)[i]->GetNumPoints(0);
3556 const int np1 = (*exp)[i]->GetNumPoints(1);
3557 const int np2 = (*exp)[i]->GetNumPoints(2);
3558 const int np01 = np0 * np1;
3559
3560 for (j = 1; j < np2; ++j)
3561 {
3562 for (k = 1; k < np1; ++k)
3563 {
3564 for (l = 1; l < np0; ++l)
3565 {
3566 outfile << cnt + (j - 1) * np01 + (k - 1) * np0 + l
3567 << " ";
3568 outfile << cnt + (j - 1) * np01 + (k - 1) * np0 + l + 1
3569 << " ";
3570 outfile << cnt + (j - 1) * np01 + k * np0 + l + 1
3571 << " ";
3572 outfile << cnt + (j - 1) * np01 + k * np0 + l << " ";
3573 outfile << cnt + j * np01 + (k - 1) * np0 + l << " ";
3574 outfile << cnt + j * np01 + (k - 1) * np0 + l + 1
3575 << " ";
3576 outfile << cnt + j * np01 + k * np0 + l + 1 << " ";
3577 outfile << cnt + j * np01 + k * np0 + l << endl;
3578 }
3579 }
3580 }
3581 cnt += np0 * np1 * np2;
3582 }
3583 }
3584 else
3585 {
3586 NEKERROR(ErrorUtil::efatal, "Not set up for this dimension");
3587 }
3588}
3589
3590/**
3591 * Write Tecplot Files Field
3592 * @param outfile Output file name.
3593 * @param expansion Expansion that is considered
3594 */
3595void ExpList::v_WriteTecplotField(std::ostream &outfile, int expansion)
3596{
3597 if (expansion == -1)
3598 {
3599 int totpoints = GetTotPoints();
3600 if (m_physState == false)
3601 {
3603 }
3604
3605 for (int i = 0; i < totpoints; ++i)
3606 {
3607 outfile << m_phys[i] << " ";
3608 if (i % 1000 == 0 && i)
3609 {
3610 outfile << std::endl;
3611 }
3612 }
3613 outfile << std::endl;
3614 }
3615 else
3616 {
3617 int nPoints = (*m_exp)[expansion]->GetTotPoints();
3618
3619 for (int i = 0; i < nPoints; ++i)
3620 {
3621 outfile << m_phys[i + m_phys_offset[expansion]] << " ";
3622 }
3623
3624 outfile << std::endl;
3625 }
3626}
3627
3628void ExpList::WriteVtkHeader(std::ostream &outfile)
3629{
3630 outfile << "<?xml version=\"1.0\"?>" << endl;
3631 outfile << R"(<VTKFile type="UnstructuredGrid" version="0.1" )"
3632 << "byte_order=\"LittleEndian\">" << endl;
3633 outfile << " <UnstructuredGrid>" << endl;
3634}
3635
3636void ExpList::WriteVtkFooter(std::ostream &outfile)
3637{
3638 outfile << " </UnstructuredGrid>" << endl;
3639 outfile << "</VTKFile>" << endl;
3640}
3641
3642void ExpList::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion,
3643 [[maybe_unused]] int istrip)
3644{
3645 int i, j, k;
3646 int nbase = (*m_exp)[expansion]->GetNumBases();
3647 int ntot = (*m_exp)[expansion]->GetTotPoints();
3648 int nquad[3];
3649
3650 int ntotminus = 1;
3651 for (i = 0; i < nbase; ++i)
3652 {
3653 nquad[i] = (*m_exp)[expansion]->GetNumPoints(i);
3654 ntotminus *= (nquad[i] - 1);
3655 }
3656
3657 Array<OneD, NekDouble> coords[3];
3658 coords[0] = Array<OneD, NekDouble>(ntot, 0.0);
3659 coords[1] = Array<OneD, NekDouble>(ntot, 0.0);
3660 coords[2] = Array<OneD, NekDouble>(ntot, 0.0);
3661 (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
3662
3663 outfile << " <Piece NumberOfPoints=\"" << ntot << "\" NumberOfCells=\""
3664 << ntotminus << "\">" << endl;
3665 outfile << " <Points>" << endl;
3666 outfile << " <DataArray type=\"Float64\" "
3667 << R"(NumberOfComponents="3" format="ascii">)" << endl;
3668 outfile << " ";
3669 for (i = 0; i < ntot; ++i)
3670 {
3671 for (j = 0; j < 3; ++j)
3672 {
3673 outfile << setprecision(8) << scientific << (float)coords[j][i]
3674 << " ";
3675 }
3676 outfile << endl;
3677 }
3678 outfile << endl;
3679 outfile << " </DataArray>" << endl;
3680 outfile << " </Points>" << endl;
3681 outfile << " <Cells>" << endl;
3682 outfile << " <DataArray type=\"Int32\" "
3683 << R"(Name="connectivity" format="ascii">)" << endl;
3684
3685 int ns = 0; // pow(2,dim) for later usage
3686 string ostr;
3687 switch (m_expType)
3688 {
3689 case e1D:
3690 {
3691 ns = 2;
3692 ostr = "3 ";
3693 for (i = 0; i < nquad[0] - 1; ++i)
3694 {
3695 outfile << i << " " << i + 1 << endl;
3696 }
3697 }
3698 break;
3699 case e2D:
3700 {
3701 ns = 4;
3702 ostr = "9 ";
3703 for (i = 0; i < nquad[0] - 1; ++i)
3704 {
3705 for (j = 0; j < nquad[1] - 1; ++j)
3706 {
3707 outfile << j * nquad[0] + i << " " << j * nquad[0] + i + 1
3708 << " " << (j + 1) * nquad[0] + i + 1 << " "
3709 << (j + 1) * nquad[0] + i << endl;
3710 }
3711 }
3712 }
3713 break;
3714 case e3D:
3715 {
3716 ns = 8;
3717 ostr = "12 ";
3718 for (i = 0; i < nquad[0] - 1; ++i)
3719 {
3720 for (j = 0; j < nquad[1] - 1; ++j)
3721 {
3722 for (k = 0; k < nquad[2] - 1; ++k)
3723 {
3724 outfile
3725 << k * nquad[0] * nquad[1] + j * nquad[0] + i << " "
3726 << k * nquad[0] * nquad[1] + j * nquad[0] + i + 1
3727 << " "
3728 << k * nquad[0] * nquad[1] + (j + 1) * nquad[0] +
3729 i + 1
3730 << " "
3731 << k * nquad[0] * nquad[1] + (j + 1) * nquad[0] + i
3732 << " "
3733 << (k + 1) * nquad[0] * nquad[1] + j * nquad[0] + i
3734 << " "
3735 << (k + 1) * nquad[0] * nquad[1] + j * nquad[0] +
3736 i + 1
3737 << " "
3738 << (k + 1) * nquad[0] * nquad[1] +
3739 (j + 1) * nquad[0] + i + 1
3740 << " "
3741 << (k + 1) * nquad[0] * nquad[1] +
3742 (j + 1) * nquad[0] + i
3743 << " " << endl;
3744 }
3745 }
3746 }
3747 }
3748 break;
3749 default:
3750 break;
3751 }
3752
3753 outfile << endl;
3754 outfile << " </DataArray>" << endl;
3755 outfile << " <DataArray type=\"Int32\" "
3756 << R"(Name="offsets" format="ascii">)" << endl;
3757 for (i = 0; i < ntotminus; ++i)
3758 {
3759 outfile << i * ns + ns << " ";
3760 }
3761 outfile << endl;
3762 outfile << " </DataArray>" << endl;
3763 outfile << " <DataArray type=\"UInt8\" "
3764 << R"(Name="types" format="ascii">)" << endl;
3765 for (i = 0; i < ntotminus; ++i)
3766 {
3767 outfile << ostr;
3768 }
3769 outfile << endl;
3770 outfile << " </DataArray>" << endl;
3771 outfile << " </Cells>" << endl;
3772 outfile << " <PointData>" << endl;
3773}
3774
3775void ExpList::WriteVtkPieceFooter(std::ostream &outfile,
3776 [[maybe_unused]] int expansion)
3777{
3778 outfile << " </PointData>" << endl;
3779 outfile << " </Piece>" << endl;
3780}
3781
3782void ExpList::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
3783 std::string var)
3784{
3785 int i;
3786 int nq = (*m_exp)[expansion]->GetTotPoints();
3787
3788 // printing the fields of that zone
3789 outfile << R"( <DataArray type="Float64" Name=")" << var << "\">"
3790 << endl;
3791 outfile << " ";
3792
3793 const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
3794
3795 for (i = 0; i < nq; ++i)
3796 {
3797 outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i])
3798 << " ";
3799 }
3800 outfile << endl;
3801 outfile << " </DataArray>" << endl;
3802}
3803
3804/**
3805 * Given a spectral/hp approximation
3806 * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
3807 * (which should be contained in #m_phys), this function calculates the
3808 * \f$L_\infty\f$ error of this approximation with respect to an exact
3809 * solution. The local distribution of the quadrature points allows an
3810 * elemental evaluation of this operation through the functions
3811 * StdRegions#StdExpansion#Linf.
3812 *
3813 * The exact solution, also evaluated at the quadrature
3814 * points, should be contained in the variable #m_phys of
3815 * the ExpList object \a Sol.
3816 *
3817 * @param soln A 1D array, containing the discrete
3818 * evaluation of the exact solution at the
3819 * quadrature points in its array #m_phys.
3820 * @return The \f$L_\infty\f$ error of the approximation.
3821 */
3823 const Array<OneD, const NekDouble> &soln)
3824{
3825 NekDouble err = 0.0;
3826
3827 if (soln == NullNekDouble1DArray)
3828 {
3829 err = Vmath::Vmax(m_npoints, inarray, 1);
3830 }
3831 else
3832 {
3833 for (int i = 0; i < m_npoints; ++i)
3834 {
3835 err = max(err, abs(inarray[i] - soln[i]));
3836 }
3837 }
3838
3839 m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
3840
3841 return err;
3842}
3843
3844/**
3845 * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
3846 * evaluated at the quadrature points (which should be contained in
3847 * #m_phys), this function calculates the \f$L_2\f$ error of this
3848 * approximation with respect to an exact solution. The local
3849 * distribution of the quadrature points allows an elemental evaluation
3850 * of this operation through the functions StdRegions#StdExpansion#L2.
3851 *
3852 * The exact solution, also evaluated at the quadrature points, should
3853 * be contained in the variable #m_phys of the ExpList object \a Sol.
3854 *
3855 * @param Sol An ExpList, containing the discrete
3856 * evaluation of the exact solution at the
3857 * quadrature points in its array #m_phys.
3858 * @return The \f$L_2\f$ error of the approximation.
3859 */
3861 const Array<OneD, const NekDouble> &soln)
3862{
3863 NekDouble err = 0.0, errl2;
3864 int i;
3865
3866 if (soln == NullNekDouble1DArray)
3867 {
3868 for (i = 0; i < (*m_exp).size(); ++i)
3869 {
3870 errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
3871 err += errl2 * errl2;
3872 }
3873 }
3874 else
3875 {
3876 for (i = 0; i < (*m_exp).size(); ++i)
3877 {
3878 errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
3879 soln + m_phys_offset[i]);
3880 err += errl2 * errl2;
3881 }
3882 }
3883
3884 m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
3885
3886 return sqrt(err);
3887}
3888
3889/**
3890 * The integration is evaluated locally, that is
3891 * \f[\int
3892 * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
3893 * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
3894 * where the integration over the separate elements is done by the
3895 * function StdRegions#StdExpansion#Integral, which discretely
3896 * evaluates the integral using Gaussian quadrature.
3897 *
3898 * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
3899 * containing the values of the function
3900 * \f$f(\boldsymbol{x})\f$ at the quadrature
3901 * points \f$\boldsymbol{x}_i\f$.
3902 * @return The value of the discretely evaluated integral
3903 * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
3904 */
3906{
3907 NekDouble sum = 0.0;
3908 int i = 0;
3909
3910 for (i = 0; i < (*m_exp).size(); ++i)
3911 {
3912 sum += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
3913 }
3914 m_comm->GetRowComm()->AllReduce(sum, LibUtilities::ReduceSum);
3915
3916 return sum;
3917}
3918
3920 const Array<OneD, Array<OneD, NekDouble>> &inarray)
3921{
3922 NekDouble flux = 0.0;
3923 int i = 0;
3924 int j;
3925
3926 for (i = 0; i < (*m_exp).size(); ++i)
3927 {
3928 Array<OneD, Array<OneD, NekDouble>> tmp(inarray.size());
3929 for (j = 0; j < inarray.size(); ++j)
3930 {
3931 tmp[j] = Array<OneD, NekDouble>(inarray[j] + m_phys_offset[i]);
3932 }
3933 flux += (*m_exp)[i]->VectorFlux(tmp);
3934 }
3935
3936 return flux;
3937}
3938
3940{
3942 "This method is not defined or valid for this class type");
3943 Array<OneD, NekDouble> NoEnergy(1, 0.0);
3944 return NoEnergy;
3945}
3946
3948{
3950 "This method is not defined or valid for this class type");
3952 return trans;
3953}
3954
3956{
3958 "This method is not defined or valid for this class type");
3959 NekDouble len = 0.0;
3960 return len;
3961}
3962
3963void ExpList::v_SetHomoLen([[maybe_unused]] const NekDouble lhom)
3964{
3966 "This method is not defined or valid for this class type");
3967}
3968
3970{
3972 "This method is not defined or valid for this class type");
3973 Array<OneD, unsigned int> NoModes(1);
3974 return NoModes;
3975}
3976
3978{
3980 "This method is not defined or valid for this class type");
3981 Array<OneD, unsigned int> NoModes(1);
3982 return NoModes;
3983}
3984
3986{
3988 "ClearGlobalLinSysManager not implemented for ExpList.");
3989}
3990
3991int ExpList::v_GetPoolCount([[maybe_unused]] std::string poolName)
3992{
3993 NEKERROR(ErrorUtil::efatal, "GetPoolCount not implemented for ExpList.");
3994 return -1;
3995}
3996
3998 [[maybe_unused]] bool clearLocalMatrices)
3999{
4001 "UnsetGlobalLinSys not implemented for ExpList.");
4002}
4003
4006{
4008 "GetGlobalLinSysManager not implemented for ExpList.");
4010}
4011
4012void ExpList::ExtractCoeffsFromFile(const std::string &fileName,
4014 const std::string &varName,
4015 Array<OneD, NekDouble> &coeffs)
4016{
4017 string varString = fileName.substr(0, fileName.find_last_of("."));
4018 int j, k, len = varString.length();
4019 varString = varString.substr(len - 1, len);
4020
4021 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
4022 std::vector<std::vector<NekDouble>> FieldData;
4023
4024 std::string ft = LibUtilities::FieldIO::GetFileType(fileName, comm);
4027 ft, comm, m_session->GetSharedFilesystem());
4028
4029 f->Import(fileName, FieldDef, FieldData);
4030
4031 bool found = false;
4032 for (j = 0; j < FieldDef.size(); ++j)
4033 {
4034 for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
4035 {
4036 if (FieldDef[j]->m_fields[k] == varName)
4037 {
4038 // Copy FieldData into locExpList
4039 ExtractDataToCoeffs(FieldDef[j], FieldData[j],
4040 FieldDef[j]->m_fields[k], coeffs);
4041 found = true;
4042 }
4043 }
4044 }
4045
4046 ASSERTL0(found, "Could not find variable '" + varName +
4047 "' in file boundary condition " + fileName);
4048}
4049
4050/**
4051 * Given a spectral/hp approximation
4052 * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
4053 * (which should be contained in #m_phys), this function calculates the
4054 * \f$H^1_2\f$ error of this approximation with respect to an exact
4055 * solution. The local distribution of the quadrature points allows an
4056 * elemental evaluation of this operation through the functions
4057 * StdRegions#StdExpansion#H1.
4058 *
4059 * The exact solution, also evaluated at the quadrature points, should
4060 * be contained in the variable #m_phys of the ExpList object \a Sol.
4061 *
4062 * @param soln An 1D array, containing the discrete evaluation
4063 * of the exact solution at the quadrature points.
4064 *
4065 * @return The \f$H^1_2\f$ error of the approximation.
4066 */
4068 const Array<OneD, const NekDouble> &soln)
4069{
4070 NekDouble err = 0.0, errh1;
4071 int i;
4072
4073 for (i = 0; i < (*m_exp).size(); ++i)
4074 {
4075 errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
4076 soln + m_phys_offset[i]);
4077 err += errh1 * errh1;
4078 }
4079
4080 m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
4081
4082 return sqrt(err);
4083}
4084
4086 std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
4087 int NumHomoDir, Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis,
4088 std::vector<NekDouble> &HomoLen, bool homoStrips,
4089 std::vector<unsigned int> &HomoSIDs, std::vector<unsigned int> &HomoZIDs,
4090 std::vector<unsigned int> &HomoYIDs)
4091{
4092 int startenum = (int)LibUtilities::eSegment;
4093 int endenum = (int)LibUtilities::eHexahedron;
4094 int s = 0;
4096
4097 ASSERTL1(NumHomoDir == HomoBasis.size(),
4098 "Homogeneous basis is not the same length as NumHomoDir");
4099 ASSERTL1(NumHomoDir == HomoLen.size(),
4100 "Homogeneous length vector is not the same length as NumHomDir");
4101
4102 // count number of shapes
4103 switch ((*m_exp)[0]->GetShapeDimension())
4104 {
4105 case 1:
4106 startenum = (int)LibUtilities::eSegment;
4107 endenum = (int)LibUtilities::eSegment;
4108 break;
4109 case 2:
4110 startenum = (int)LibUtilities::eTriangle;
4111 endenum = (int)LibUtilities::eQuadrilateral;
4112 break;
4113 case 3:
4114 startenum = (int)LibUtilities::eTetrahedron;
4115 endenum = (int)LibUtilities::eHexahedron;
4116 break;
4117 }
4118
4119 for (s = startenum; s <= endenum; ++s)
4120 {
4121 std::vector<unsigned int> elementIDs;
4122 std::vector<LibUtilities::BasisType> basis;
4123 std::vector<unsigned int> numModes;
4124 std::vector<std::string> fields;
4125
4126 bool first = true;
4127 bool UniOrder = true;
4128 int n;
4129
4130 shape = (LibUtilities::ShapeType)s;
4131
4132 for (int i = 0; i < (*m_exp).size(); ++i)
4133 {
4134 if ((*m_exp)[i]->GetGeom()->GetShapeType() == shape)
4135 {
4136 elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
4137 if (first)
4138 {
4139 for (int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
4140 {
4141 basis.push_back(
4142 (*m_exp)[i]->GetBasis(j)->GetBasisType());
4143 numModes.push_back(
4144 (*m_exp)[i]->GetBasis(j)->GetNumModes());
4145 }
4146
4147 // add homogeneous direction details if defined
4148 for (n = 0; n < NumHomoDir; ++n)
4149 {
4150 basis.push_back(HomoBasis[n]->GetBasisType());
4151 numModes.push_back(HomoBasis[n]->GetNumModes());
4152 }
4153
4154 first = false;
4155 }
4156 else
4157 {
4158 ASSERTL0(
4159 (*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],
4160 "Routine is not set up for multiple bases definitions");
4161
4162 for (int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
4163 {
4164 numModes.push_back(
4165 (*m_exp)[i]->GetBasis(j)->GetNumModes());
4166 if (numModes[j] !=
4167 (*m_exp)[i]->GetBasis(j)->GetNumModes())
4168 {
4169 UniOrder = false;
4170 }
4171 }
4172 // add homogeneous direction details if defined
4173 for (n = 0; n < NumHomoDir; ++n)
4174 {
4175 numModes.push_back(HomoBasis[n]->GetNumModes());
4176 }
4177 }
4178 }
4179 }
4180
4181 if (elementIDs.size() > 0)
4182 {
4185 AllocateSharedPtr(shape, elementIDs, basis, UniOrder,
4186 numModes, fields, NumHomoDir, HomoLen,
4187 homoStrips, HomoSIDs, HomoZIDs, HomoYIDs);
4188 fielddef.push_back(fdef);
4189 }
4190 }
4191}
4192
4193//
4194// Virtual functions
4195//
4196std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::
4198{
4199 std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
4200 v_GetFieldDefinitions(returnval);
4201 return returnval;
4202}
4203
4205 std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
4206{
4208}
4209
4210// Append the element data listed in elements
4211// fielddef->m_ElementIDs onto fielddata
4214 std::vector<NekDouble> &fielddata)
4215{
4216 v_AppendFieldData(fielddef, fielddata, m_coeffs);
4217}
4218
4221 std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
4222{
4223 int i;
4224 // Determine mapping from element ids to location in
4225 // expansion list
4226 // Determine mapping from element ids to location in
4227 // expansion list
4228 map<int, int> ElmtID_to_ExpID;
4229 for (i = 0; i < (*m_exp).size(); ++i)
4230 {
4231 ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
4232 }
4233
4234 for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
4235 {
4236 int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
4237 int datalen = (*m_exp)[eid]->GetNcoeffs();
4238 if ((*m_exp)[eid]->IsNodalNonTensorialExp())
4239 {
4240 // need to convert nodal coeff values into orthonormal expansion
4241 Array<OneD, NekDouble> orthocoeffs((*m_exp)[eid]->GetNcoeffs());
4242 (*m_exp)[eid]->NodalToModal(coeffs + m_coeff_offset[eid],
4243 orthocoeffs);
4244 fielddata.insert(fielddata.end(), &orthocoeffs[0],
4245 &orthocoeffs[0] + datalen);
4246 }
4247 else
4248 {
4249 fielddata.insert(fielddata.end(), &coeffs[m_coeff_offset[eid]],
4250 &coeffs[m_coeff_offset[eid]] + datalen);
4251 }
4252 }
4253}
4254
4255/// Extract the data in fielddata into the coeffs
4258 std::vector<NekDouble> &fielddata, std::string &field,
4259 Array<OneD, NekDouble> &coeffs, std::unordered_map<int, int> zIdToPlane)
4260{
4261 v_ExtractDataToCoeffs(fielddef, fielddata, field, coeffs, zIdToPlane);
4262}
4263
4265 const std::shared_ptr<ExpList> &fromExpList,
4266 const Array<OneD, const NekDouble> &fromCoeffs,
4267 Array<OneD, NekDouble> &toCoeffs)
4268{
4269 v_ExtractCoeffsToCoeffs(fromExpList, fromCoeffs, toCoeffs);
4270}
4271
4272/**
4273 * @brief Extract data from raw field data into expansion list.
4274 *
4275 * @param fielddef Field definitions.
4276 * @param fielddata Data for associated field.
4277 * @param field Field variable name.
4278 * @param coeffs Resulting coefficient array.
4279 */
4282 std::vector<NekDouble> &fielddata, std::string &field,
4283 Array<OneD, NekDouble> &coeffs,
4284 [[maybe_unused]] std::unordered_map<int, int> zIdToPlane)
4285{
4286 int i, expId;
4287 int offset = 0;
4288 int modes_offset = 0;
4289 int datalen = fielddata.size() / fielddef->m_fields.size();
4290
4291 // Find data location according to field definition
4292 for (i = 0; i < fielddef->m_fields.size(); ++i)
4293 {
4294 if (fielddef->m_fields[i] == field)
4295 {
4296 break;
4297 }
4298 offset += datalen;
4299 }
4300
4301 ASSERTL0(i != fielddef->m_fields.size(),
4302 "Field (" + field + ") not found in file.");
4303
4304 if (m_elmtToExpId.size() == 0)
4305 {
4306 // Loop in reverse order so that in case where using a
4307 // Homogeneous expansion it sets geometry ids to first part of
4308 // m_exp list. Otherwise will set to second (complex) expansion
4309 for (i = (*m_exp).size() - 1; i >= 0; --i)
4310 {
4311 m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
4312 }
4313 }
4314
4315 for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
4316 {
4317 // Reset modes_offset in the case where all expansions of
4318 // the same order.
4319 if (fielddef->m_uniOrder == true)
4320 {
4321 modes_offset = 0;
4322 }
4323
4325 fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
4326
4327 const int elmtId = fielddef->m_elementIDs[i];
4328 auto eIt = m_elmtToExpId.find(elmtId);
4329
4330 if (eIt == m_elmtToExpId.end())
4331 {
4332 offset += datalen;
4333 modes_offset += (*m_exp)[0]->GetNumBases();
4334 continue;
4335 }
4336
4337 expId = eIt->second;
4338
4339 bool sameBasis = true;
4340 for (int j = 0; j < fielddef->m_basis.size(); ++j)
4341 {
4342 if (fielddef->m_basis[j] != (*m_exp)[expId]->GetBasisType(j))
4343 {
4344 sameBasis = false;
4345 break;
4346 }
4347 }
4348
4349 if (datalen == (*m_exp)[expId]->GetNcoeffs() && sameBasis)
4350 {
4351 Vmath::Vcopy(datalen, &fielddata[offset], 1,
4352 &coeffs[m_coeff_offset[expId]], 1);
4353 }
4354 else
4355 {
4356 (*m_exp)[expId]->ExtractDataToCoeffs(
4357 &fielddata[offset], fielddef->m_numModes, modes_offset,
4358 &coeffs[m_coeff_offset[expId]], fielddef->m_basis);
4359 }
4360
4361 offset += datalen;
4362 modes_offset += (*m_exp)[0]->GetNumBases();
4363 }
4364
4365 return;
4366}
4367
4369 const std::shared_ptr<ExpList> &fromExpList,
4370 const Array<OneD, const NekDouble> &fromCoeffs,
4371 Array<OneD, NekDouble> &toCoeffs)
4372{
4373 int i;
4374 int offset = 0;
4375
4376 map<int, int> GidToEid;
4377
4378 for (i = 0; i < (*m_exp).size(); ++i)
4379 {
4380 GidToEid[fromExpList->GetExp(i)->GetGeom()->GetGlobalID()] = i;
4381 }
4382
4383 for (i = 0; i < (*m_exp).size(); ++i)
4384 {
4385 std::vector<unsigned int> nummodes;
4386 vector<LibUtilities::BasisType> basisTypes;
4387
4388 int eid = GidToEid[(*m_exp)[i]->GetGeom()->GetGlobalID()];
4389 for (int j = 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j)
4390 {
4391 nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j));
4392 basisTypes.push_back(fromExpList->GetExp(eid)->GetBasisType(j));
4393 }
4394
4395 offset = fromExpList->GetCoeff_Offset(eid);
4396 (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes, 0,
4397 &toCoeffs[m_coeff_offset[i]],
4398 basisTypes);
4399 }
4400}
4401
4402/**
4403 * Get the weight value on boundaries
4404 */
4406 Array<OneD, NekDouble> &weightJump)
4407{
4408 size_t nTracePts = weightAver.size();
4409 // average for interior traces
4410 for (int i = 0; i < nTracePts; ++i)
4411 {
4412 weightAver[i] = 0.5;
4413 weightJump[i] = 1.0;
4414 }
4415 FillBwdWithBwdWeight(weightAver, weightJump);
4416}
4417
4419 const Array<OneD, const NekDouble> &CircCentre,
4420 Array<OneD, Array<OneD, NekDouble>> &outarray)
4421{
4422 int npts;
4423
4424 int MFdim = 3;
4425 int nq = outarray[0].size() / MFdim;
4426
4427 // Assume whole array is of same coordinate dimension
4428 int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
4429
4430 Array<OneD, Array<OneD, NekDouble>> MFloc(MFdim * coordim);
4431 // Process each expansion.
4432 for (int i = 0; i < m_exp->size(); ++i)
4433 {
4434 npts = (*m_exp)[i]->GetTotPoints();
4435
4436 for (int j = 0; j < MFdim * coordim; ++j)
4437 {
4438 MFloc[j] = Array<OneD, NekDouble>(npts, 0.0);
4439 }
4440
4441 // MF from LOCALREGIONS
4442 (*m_exp)[i]->GetMetricInfo()->GetMovingFrames(
4443 (*m_exp)[i]->GetPointsKeys(), MMFdir, CircCentre, MFloc);
4444
4445 // Get the physical data offset for this expansion.
4446 for (int j = 0; j < MFdim; ++j)
4447 {
4448 for (int k = 0; k < coordim; ++k)
4449 {
4450 Vmath::Vcopy(npts, &MFloc[j * coordim + k][0], 1,
4451 &outarray[j][k * nq + m_phys_offset[i]], 1);
4452 }
4453 }
4454 }
4455}
4456
4457/**
4458 * @brief Generate vector v such that v[i] = scalar1 if i is in the
4459 * element < ElementID. Otherwise, v[i] = scalar2.
4460 *
4461 */
4462void ExpList::GenerateElementVector(const int ElementID,
4463 const NekDouble scalar1,
4464 const NekDouble scalar2,
4465 Array<OneD, NekDouble> &outarray)
4466{
4467 int npoints_e;
4468 NekDouble coeff;
4469
4470 Array<OneD, NekDouble> outarray_e;
4471
4472 for (int i = 0; i < (*m_exp).size(); ++i)
4473 {
4474 npoints_e = (*m_exp)[i]->GetTotPoints();
4475
4476 if (i <= ElementID)
4477 {
4478 coeff = scalar1;
4479 }
4480 else
4481 {
4482 coeff = scalar2;
4483 }
4484
4485 outarray_e = Array<OneD, NekDouble>(npoints_e, coeff);
4486 Vmath::Vcopy(npoints_e, &outarray_e[0], 1, &outarray[m_phys_offset[i]],
4487 1);
4488 }
4489}
4490
4493{
4495 "This method is not defined or valid for this class type");
4497 return result;
4498}
4499
4500std::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(
4501 [[maybe_unused]] int i)
4502{
4504 "This method is not defined or valid for this class type");
4505 static std::shared_ptr<ExpList> result;
4506 return result;
4507}
4508
4509/**
4510 * Upwind the left and right states given by the Arrays Fwd and Bwd
4511 * using the vector quantity Vec and ouput the upwinded value in the
4512 * array upwind.
4513 *
4514 * @param Vec Velocity field.
4515 * @param Fwd Left state.
4516 * @param Bwd Right state.
4517 * @param Upwind Output vector.
4518 */
4522 Array<OneD, NekDouble> &Upwind)
4523{
4524 switch (m_expType)
4525 {
4526 case e1D:
4527 {
4528 int i, j, k, e_npoints, offset;
4529 Array<OneD, NekDouble> normals;
4530 NekDouble Vn;
4531
4532 // Assume whole array is of same coordimate dimension
4533 int coordim = GetCoordim(0);
4534
4535 ASSERTL1(Vec.size() >= coordim,
4536 "Input vector does not have sufficient dimensions to "
4537 "match coordim");
4538
4539 // Process each expansion
4540 for (i = 0; i < m_exp->size(); ++i)
4541 {
4542 // Get the number of points in the expansion and the normals.
4543 e_npoints = (*m_exp)[i]->GetNumPoints(0);
4544 normals = (*m_exp)[i]->GetPhysNormals();
4545
4546 // Get the physical data offset of the expansion in m_phys.
4547 offset = m_phys_offset[i];
4548
4549 // Compute each data point.
4550 for (j = 0; j < e_npoints; ++j)
4551 {
4552 // Calculate normal velocity.
4553 Vn = 0.0;
4554 for (k = 0; k < coordim; ++k)
4555 {
4556 Vn += Vec[k][offset + j] * normals[k * e_npoints + j];
4557 }
4558
4559 // Upwind based on direction of normal velocity.
4560 if (Vn > 0.0)
4561 {
4562 Upwind[offset + j] = Fwd[offset + j];
4563 }
4564 else
4565 {
4566 Upwind[offset + j] = Bwd[offset + j];
4567 }
4568 }
4569 }
4570 }
4571 break;
4572 default:
4574 "This method is not defined or valid for this class type");
4575 break;
4576 }
4577}
4578
4579/**
4580 * One-dimensional upwind.
4581 * \see ExpList::Upwind(
4582 * const Array<OneD, const Array<OneD, NekDouble> >,
4583 * const Array<OneD, const NekDouble>,
4584 * const Array<OneD, const NekDouble>,
4585 * Array<OneD, NekDouble>, int)
4586 *
4587 * @param Vn Velocity field.
4588 * @param Fwd Left state.
4589 * @param Bwd Right state.
4590 * @param Upwind Output vector.
4591 */
4595 Array<OneD, NekDouble> &Upwind)
4596{
4597 ASSERTL1(Vn.size() >= m_npoints, "Vn is not of sufficient length");
4598 ASSERTL1(Fwd.size() >= m_npoints, "Fwd is not of sufficient length");
4599 ASSERTL1(Bwd.size() >= m_npoints, "Bwd is not of sufficient length");
4600 ASSERTL1(Upwind.size() >= m_npoints, "Upwind is not of sufficient length");
4601
4602 // Process each point in the expansion.
4603 for (int j = 0; j < m_npoints; ++j)
4604 {
4605 // Upwind based on one-dimensional velocity.
4606 if (Vn[j] > 0.0)
4607 {
4608 Upwind[j] = Fwd[j];
4609 }
4610 else
4611 {
4612 Upwind[j] = Bwd[j];
4613 }
4614 }
4615}
4616
4617std::shared_ptr<ExpList> &ExpList::v_GetTrace()
4618{
4620 "This method is not defined or valid for this class type");
4621 static std::shared_ptr<ExpList> returnVal;
4622 return returnVal;
4623}
4624
4625std::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
4626{
4628 "This method is not defined or valid for this class type");
4629 static std::shared_ptr<AssemblyMapDG> result;
4630 return result;
4631}
4632
4633std::shared_ptr<InterfaceMapDG> &ExpList::v_GetInterfaceMap()
4634{
4636 "This method is not defined or valid for this class type");
4637 static std::shared_ptr<InterfaceMapDG> result;
4638 return result;
4639}
4640
4642{
4643 return GetTraceMap()->GetBndCondIDToGlobalTraceID();
4644}
4645
4647{
4649 "This method is not defined or valid for this class type");
4650 static std::vector<bool> result;
4651 return result;
4652}
4653
4654/**
4655 * @brief Helper function to re-align face to a given orientation.
4656 */
4657void AlignFace(const StdRegions::Orientation orient, const int nquad1,
4658 const int nquad2, const Array<OneD, const NekDouble> &in,
4660{
4661 // Copy transpose.
4666 {
4667 for (int i = 0; i < nquad2; ++i)
4668 {
4669 for (int j = 0; j < nquad1; ++j)
4670 {
4671 out[i * nquad1 + j] = in[j * nquad2 + i];
4672 }
4673 }
4674 }
4675 else
4676 {
4677 for (int i = 0; i < nquad2; ++i)
4678 {
4679 for (int j = 0; j < nquad1; ++j)
4680 {
4681 out[i * nquad1 + j] = in[i * nquad1 + j];
4682 }
4683 }
4684 }
4685
4690 {
4691 // Reverse x direction
4692 for (int i = 0; i < nquad2; ++i)
4693 {
4694 for (int j = 0; j < nquad1 / 2; ++j)
4695 {
4696 swap(out[i * nquad1 + j], out[i * nquad1 + nquad1 - j - 1]);
4697 }
4698 }
4699 }
4700
4705 {
4706 // Reverse y direction
4707 for (int j = 0; j < nquad1; ++j)
4708 {
4709 for (int i = 0; i < nquad2 / 2; ++i)
4710 {
4711 swap(out[i * nquad1 + j], out[(nquad2 - i - 1) * nquad1 + j]);
4712 }
4713 }
4714 }
4715}
4716
4717/**
4718 * For each local element, copy the normals stored in the element list
4719 * into the array \a normals. This function should only be called by a
4720 * trace explist, which has setup the left and right adjacent elements.
4721 * @param normals Two dimensional array in which to copy normals
4722 * to. The first dimension is the coordim. The
4723 * second dimension is the same size as trace phys
4724 * space.
4725 */
4727{
4728 int i, j, k, e_npoints, offset;
4730
4731 // Assume whole array is of same coordinate dimension
4732 int coordim = GetCoordim(0);
4733
4734 ASSERTL1(normals.size() >= coordim,
4735 "Output vector does not have sufficient dimensions to "
4736 "match coordim");
4737
4738 switch (m_expType)
4739 {
4740 case e0D:
4741 {
4742 // Process each expansion.
4743 for (i = 0; i < m_exp->size(); ++i)
4744 {
4745 LocalRegions::ExpansionSharedPtr loc_exp = (*m_exp)[i];
4746
4748 loc_exp->GetLeftAdjacentElementExp();
4749
4750 // Get the number of points and normals for this expansion.
4751 e_npoints = 1;
4752 locnormals = loc_elmt->GetTraceNormal(
4753 loc_exp->GetLeftAdjacentElementTrace());
4754
4755 // Get the physical data offset for this expansion.
4756 offset = m_phys_offset[i];
4757
4758 // Process each point in the expansion.
4759 for (j = 0; j < e_npoints; ++j)
4760 {
4761 // Process each spatial dimension and copy the
4762 // values into the output array.
4763 for (k = 0; k < coordim; ++k)
4764 {
4765 normals[k][offset] = locnormals[k][0];
4766 }
4767 }
4768 }
4769 }
4770 break;
4771 case e1D:
4772 {
4773 // Process each (trace) expansion.
4774 for (i = 0; i < m_exp->size(); ++i)
4775 {
4776 LocalRegions::ExpansionSharedPtr traceExp = (*m_exp)[i];
4777 // location of this normal vector in the output array.
4778 int offset = m_phys_offset[i];
4779
4780 // Get number of points from left expansion.
4782 traceExp->GetLeftAdjacentElementExp();
4783 int edgeId = traceExp->GetLeftAdjacentElementTrace();
4784 LibUtilities::PointsKey edgePoints =
4785 exp2D->GetTraceBasisKey(edgeId).GetPointsKey();
4786 LibUtilities::PointsKey tracePoints =
4787 traceExp->GetBasis(0)->GetPointsKey();
4788
4789 // If right adjacent element exists, then we compare
4790 // the left and right side and take the one with
4791 // the highest number of points to compute the
4792 // local normals. However, it's a question whether
4793 // this effort pays off.
4794 bool useRight = false;
4795 if (traceExp->GetRightAdjacentElementTrace() >= 0)
4796 {
4798 traceExp->GetRightAdjacentElementExp();
4799 int RedgeId = traceExp->GetRightAdjacentElementTrace();
4800 LibUtilities::PointsKey RedgePoints =
4801 Rexp2D->GetTraceBasisKey(RedgeId).GetPointsKey();
4802
4803 if (RedgePoints.GetNumPoints() > edgePoints.GetNumPoints())
4804 {
4805 exp2D = Rexp2D;
4806 edgeId = RedgeId;
4807 edgePoints = RedgePoints;
4808 useRight = true;
4809 }
4810 }
4811
4812 const Array<OneD, const Array<OneD, NekDouble>> &locNormals =
4813 exp2D->GetTraceNormal(edgeId);
4814
4815 // For unknown reason, GetTraceNormal(2D) returns normals
4816 // that has been reoriented to trace order.
4817 // So here we don't need to reorient them again.
4818 for (int d = 0; d < coordim; ++d)
4819 {
4820 LibUtilities::Interp1D(edgePoints, locNormals[d].data(),
4821 tracePoints,
4822 normals[d].data() + offset);
4823 // Trace normal direction is always the outward
4824 // direction of the left element.
4825 if (useRight)
4826 {
4827 Vmath::Neg((int)tracePoints.GetNumPoints(),
4828 &normals[d][offset], 1);
4829 }
4830 }
4831 }
4832 }
4833 break;
4834 case e2D:
4835 {
4837
4838 // Process each expansion.
4839 for (i = 0; i < m_exp->size(); ++i)
4840 {
4841 LocalRegions::ExpansionSharedPtr traceExp = (*m_exp)[i];
4842 // location of this normal vector in the output array.
4843 int offset = m_phys_offset[i];
4844
4845 // Get the normals from left expansion.
4846 // NOTE:
4847 // One can choose to compare the left and right side and take
4848 // the one with the highest number of points to compute the
4849 // local normals. Here are 2 reasons why we don't do so:
4850 // 1.
4851 // in general two adjacent elements must share a common
4852 // cuerved edge/face, which can be precisely described even by
4853 // the lower-order side. Even if the two sides are not exactly
4854 // the same, it should not affect the convergence of solution.
4855 // 2.
4856 // In 3D, it's hard to define which is side has higher order.
4857 // The left-side may have higher order in axis 0 but lower in
4858 // axis 1. It's too complicated and not worth the effort.
4860 traceExp->GetLeftAdjacentElementExp();
4861 int faceId = traceExp->GetLeftAdjacentElementTrace();
4862 const Array<OneD, const Array<OneD, NekDouble>> &locNormals =
4863 exp3D->GetTraceNormal(faceId);
4864
4865 StdRegions::Orientation orient = exp3D->GetTraceOrient(faceId);
4866
4867 // swap local basiskey 0 and 1 if orientation is transposed
4868 // (>=9)
4869 int fromid0, fromid1;
4870
4872 {
4873 fromid0 = 0;
4874 fromid1 = 1;
4875 }
4876 else
4877 {
4878 fromid0 = 1;
4879 fromid1 = 0;
4880 }
4881
4882 LibUtilities::BasisKey faceBasis0 =
4883 exp3D->GetTraceBasisKey(faceId, fromid0);
4884 LibUtilities::BasisKey faceBasis1 =
4885 exp3D->GetTraceBasisKey(faceId, fromid1);
4886 LibUtilities::BasisKey traceBasis0 =
4887 traceExp->GetBasis(0)->GetBasisKey();
4888 LibUtilities::BasisKey traceBasis1 =
4889 traceExp->GetBasis(1)->GetBasisKey();
4890
4891 const int faceNq0 = faceBasis0.GetNumPoints();
4892 const int faceNq1 = faceBasis1.GetNumPoints();
4893
4894 // Reorient normals from stdExp definition onto the same
4895 // orientation as the trace expansion.(also match the
4896 // swapped local basiskey)
4897 Array<OneD, int> map;
4898 exp3D->ReOrientTracePhysMap(orient, map, faceNq0, faceNq1);
4899
4900 // Perform reorientation and interpolation.
4901 Array<OneD, NekDouble> traceNormals(faceNq0 * faceNq1);
4902 for (j = 0; j < coordim; ++j)
4903 {
4904 Vmath::Scatr(faceNq0 * faceNq1, locNormals[j], map,
4905 traceNormals);
4906
4908 faceBasis0.GetPointsKey(), faceBasis1.GetPointsKey(),
4909 traceNormals, traceBasis0.GetPointsKey(),
4910 traceBasis1.GetPointsKey(), tmp = normals[j] + offset);
4911 }
4912 }
4913 }
4914 break;
4915 default:
4916 {
4918 "This method is not defined or valid for this class type");
4919 }
4920 }
4921}
4922
4923/**
4924 * Returns the element normal length for each trace expansion. The array
4925 * has the same size as trace phys space. This function should only be
4926 * called by a trace explist, which has setup the left and right adjacent
4927 * elements.
4928 * This function is only used by DiffusionIP to commpute the penalty.
4929 * However, it's a question whether we need to calculate the lengthFwd
4930 * and lengthBwd separately, since in most cases, they are equavalent.
4931 * Same logic applies to v_GetNormals().
4932 * @param lengthsFwd Output array of normal lengths for left side.
4933 * @param lengthsBwd Output array of normal lengths for right side.
4934 */
4936 Array<OneD, NekDouble> &lengthsBwd)
4937{
4938 int e_npoints;
4939
4940 Array<OneD, NekDouble> locLeng;
4943 Array<OneD, int> LRbndnumbs(2);
4945 lengLR[0] = lengthsFwd;
4946 lengLR[1] = lengthsBwd;
4950 int e_npoints0 = -1;
4951 if (m_expType == e1D)
4952 {
4953 for (int i = 0; i < m_exp->size(); ++i)
4954 {
4955 loc_exp = (*m_exp)[i];
4956 int offset = m_phys_offset[i];
4957
4958 e_npoints = (*m_exp)[i]->GetNumPoints(0);
4959 if (e_npoints0 < e_npoints)
4960 {
4961 for (int nlr = 0; nlr < 2; nlr++)
4962 {
4963 lengintp[nlr] = Array<OneD, NekDouble>(e_npoints, 0.0);
4964 }
4965 e_npoints0 = e_npoints;
4966 }
4967
4968 LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4969 LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4970
4971 LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4972 LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4973 for (int nlr = 0; nlr < 2; ++nlr)
4974 {
4975 Vmath::Zero(e_npoints0, lengintp[nlr], 1);
4976 lengAdd[nlr] = lengintp[nlr];
4977 int bndNumber = LRbndnumbs[nlr];
4978 loc_elmt = LRelmts[nlr];
4979 if (bndNumber >= 0)
4980 {
4981 locLeng = loc_elmt->GetElmtBndNormDirElmtLen(bndNumber);
4982
4984 loc_exp->GetBasis(0)->GetPointsKey();
4985 LibUtilities::PointsKey from_key =
4986 loc_elmt->GetTraceBasisKey(bndNumber).GetPointsKey();
4987
4988 // For unknown reason, GetTraceNormal(2D) returns normals
4989 // that has been reoriented to trace order.
4990 // So here we don't need to reorient them again.
4991
4992 // Always do interpolation
4993 LibUtilities::Interp1D(from_key, locLeng, to_key,
4994 lengintp[nlr]);
4995 lengAdd[nlr] = lengintp[nlr];
4996 }
4997
4998 for (int j = 0; j < e_npoints; ++j)
4999 {
5000 lengLR[nlr][offset + j] = lengAdd[nlr][j];
5001 }
5002 }
5003 }
5004 }
5005 else if (m_expType == e2D)
5006 {
5007 for (int i = 0; i < m_exp->size(); ++i)
5008 {
5009 loc_exp = (*m_exp)[i];
5010 int offset = m_phys_offset[i];
5011
5012 LibUtilities::BasisKey traceBasis0 =
5013 loc_exp->GetBasis(0)->GetBasisKey();
5014 LibUtilities::BasisKey traceBasis1 =
5015 loc_exp->GetBasis(1)->GetBasisKey();
5016 const int TraceNq0 = traceBasis0.GetNumPoints();
5017 const int TraceNq1 = traceBasis1.GetNumPoints();
5018 e_npoints = TraceNq0 * TraceNq1;
5019 if (e_npoints0 < e_npoints)
5020 {
5021 for (int nlr = 0; nlr < 2; nlr++)
5022 {
5023 lengintp[nlr] = Array<OneD, NekDouble>(e_npoints, 0.0);
5024 }
5025 e_npoints0 = e_npoints;
5026 }
5027
5028 LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
5029 LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
5030
5031 LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
5032 LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
5033 for (int nlr = 0; nlr < 2; ++nlr)
5034 {
5035 Vmath::Zero(e_npoints0, lengintp[nlr], 1);
5036 int bndNumber = LRbndnumbs[nlr];
5037 loc_elmt = LRelmts[nlr];
5038 if (bndNumber >= 0)
5039 {
5040 locLeng = loc_elmt->GetElmtBndNormDirElmtLen(bndNumber);
5041 // Project normals from 3D element onto the
5042 // same orientation as the trace expansion.
5044 loc_elmt->GetTraceOrient(bndNumber);
5045
5046 int fromid0, fromid1;
5048 {
5049 fromid0 = 0;
5050 fromid1 = 1;
5051 }
5052 else
5053 {
5054 fromid0 = 1;
5055 fromid1 = 0;
5056 }
5057
5058 LibUtilities::BasisKey faceBasis0 =
5059 loc_elmt->GetTraceBasisKey(bndNumber, fromid0);
5060 LibUtilities::BasisKey faceBasis1 =
5061 loc_elmt->GetTraceBasisKey(bndNumber, fromid1);
5062 const int faceNq0 = faceBasis0.GetNumPoints();
5063 const int faceNq1 = faceBasis1.GetNumPoints();
5064 Array<OneD, NekDouble> alignedLeng(faceNq0 * faceNq1);
5065
5066 AlignFace(orient, faceNq0, faceNq1, locLeng, alignedLeng);
5068 faceBasis0.GetPointsKey(), faceBasis1.GetPointsKey(),
5069 alignedLeng, traceBasis0.GetPointsKey(),
5070 traceBasis1.GetPointsKey(), lengintp[nlr]);
5071 }
5072
5073 for (int j = 0; j < e_npoints; ++j)
5074 {
5075 lengLR[nlr][offset + j] = lengintp[nlr][j];
5076 }
5077 }
5078 }
5079 }
5080}
5081
5083 [[maybe_unused]] const Array<OneD, const NekDouble> &Fn,
5084 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5085{
5087 "This method is not defined or valid for this class type");
5088}
5089
5091 [[maybe_unused]] const Array<OneD, const NekDouble> &Fwd,
5092 [[maybe_unused]] const Array<OneD, const NekDouble> &Bwd,
5093 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5094{
5096 "This method is not defined or valid for this class type");
5097}
5098
5100 [[maybe_unused]] Array<OneD, NekDouble> &Bwd)
5101{
5103 "This method is not defined or valid for this class type");
5104}
5105
5107 [[maybe_unused]] const Array<OneD, const NekDouble> &field,
5108 [[maybe_unused]] Array<OneD, NekDouble> &Fwd,
5109 [[maybe_unused]] Array<OneD, NekDouble> &Bwd, [[maybe_unused]] bool FillBnd,
5110 [[maybe_unused]] bool PutFwdInBwdOnBCs, [[maybe_unused]] bool DoExchange)
5111{
5113 "This method is not defined or valid for this class type");
5114}
5115
5117 [[maybe_unused]] const Array<OneD, NekDouble> &Fwd,
5118 [[maybe_unused]] Array<OneD, NekDouble> &Bwd,
5119 [[maybe_unused]] bool PutFwdInBwdOnBCs)
5120{
5122 "This method is not defined or valid for this class type");
5123}
5124
5126 [[maybe_unused]] const Array<OneD, const NekDouble> &Fwd,
5127 [[maybe_unused]] const Array<OneD, const NekDouble> &Bwd,
5128 [[maybe_unused]] Array<OneD, NekDouble> &field)
5129{
5131 "v_AddTraceQuadPhysToField is not defined for this class type");
5132}
5133
5135 [[maybe_unused]] const Array<OneD, const NekDouble> &Fwd,
5136 [[maybe_unused]] const Array<OneD, const NekDouble> &Bwd,
5137 [[maybe_unused]] Array<OneD, NekDouble> &field)
5138{
5140 "v_AddTraceQuadPhysToOffDiag is not defined for this class");
5141}
5142
5144 [[maybe_unused]] const Array<OneD, const NekDouble> &Fwd,
5145 [[maybe_unused]] const Array<OneD, const NekDouble> &Bwd,
5146 [[maybe_unused]] Array<OneD, NekDouble> &locTraceFwd,
5147 [[maybe_unused]] Array<OneD, NekDouble> &locTraceBwd)
5148{
5150 "v_GetLocTraceFromTracePts is not defined for this class");
5151}
5152
5154{
5156 "v_GetBndCondBwdWeight is not defined for this class type");
5157 static Array<OneD, NekDouble> tmp;
5158 return tmp;
5159}
5160
5161void ExpList::v_SetBndCondBwdWeight([[maybe_unused]] const int index,
5162 [[maybe_unused]] const NekDouble value)
5163{
5165 "v_setBndCondBwdWeight is not defined for this class type");
5166}
5167
5168const vector<bool> &ExpList::v_GetLeftAdjacentFaces(void) const
5169{
5171 "This method is not defined or valid for this class type");
5172 static vector<bool> tmp;
5173 return tmp;
5174}
5175
5177 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5178{
5180 "This method is not defined or valid for this class type");
5181}
5182
5184 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5185 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5186{
5188 "This method is not defined or valid for this class type");
5189}
5190
5192 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5193 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5194{
5196 "This method is not defined or valid for this class type");
5197}
5198
5200 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5201 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5202 [[maybe_unused]] const StdRegions::ConstFactorMap &factors,
5203 [[maybe_unused]] const StdRegions::VarCoeffMap &varcoeff,
5204 [[maybe_unused]] const MultiRegions::VarFactorsMap &varfactors,
5205 [[maybe_unused]] const Array<OneD, const NekDouble> &dirForcing,
5206 [[maybe_unused]] const bool PhysSpaceForcing)
5207{
5208 NEKERROR(ErrorUtil::efatal, "HelmSolve not implemented.");
5209 return NullGlobalLinSysKey;
5210}
5211
5213 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5214 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5215 [[maybe_unused]] const StdRegions::ConstFactorMap &factors,
5216 [[maybe_unused]] const StdRegions::VarCoeffMap &varcoeff,
5217 [[maybe_unused]] const MultiRegions::VarFactorsMap &varfactors,
5218 [[maybe_unused]] const Array<OneD, const NekDouble> &dirForcing,
5219 [[maybe_unused]] const bool PhysSpaceForcing)
5220{
5222 "LinearAdvectionDiffusionReactionSolve not implemented.");
5223 return NullGlobalLinSysKey;
5224}
5225
5227 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &velocity,
5228 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5229 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5230 [[maybe_unused]] const NekDouble lambda,
5231 [[maybe_unused]] const Array<OneD, const NekDouble> &dirForcing)
5232{
5234 "This method is not defined or valid for this class type");
5235}
5236
5238 [[maybe_unused]] const int npts,
5239 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5240 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5241 [[maybe_unused]] bool Shuff, [[maybe_unused]] bool UnShuff)
5242{
5244 "This method is not defined or valid for this class type");
5245}
5246
5248 [[maybe_unused]] const int npts,
5249 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5250 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5251 [[maybe_unused]] bool Shuff, [[maybe_unused]] bool UnShuff)
5252{
5254 "This method is not defined or valid for this class type");
5255}
5256
5258 [[maybe_unused]] const int npts,
5259 [[maybe_unused]] const Array<OneD, NekDouble> &inarray1,
5260 [[maybe_unused]] const Array<OneD, NekDouble> &inarray2,
5261 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5262{
5264 "This method is not defined or valid for this class type");
5265}
5266
5268 [[maybe_unused]] const int npts,
5269 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray1,
5270 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray2,
5271 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &outarray)
5272{
5274 "This method is not defined or valid for this class type");
5275}
5276
5278 [[maybe_unused]] Array<OneD, NekDouble> &BndVals,
5279 [[maybe_unused]] const Array<OneD, NekDouble> &TotField,
5280 [[maybe_unused]] int BndID)
5281{
5283 "This method is not defined or valid for this class type");
5284}
5285
5287 [[maybe_unused]] Array<OneD, const NekDouble> &V1,
5288 [[maybe_unused]] Array<OneD, const NekDouble> &V2,
5289 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5290 [[maybe_unused]] int BndID)
5291{
5293 "This method is not defined or valid for this class type");
5294}
5295
5298{
5300 switch (GetCoordim(0))
5301 {
5302 case 1:
5303 {
5304 for (int i = 0; i < GetExpSize(); ++i)
5305 {
5306 (*m_exp)[i]->NormVectorIProductWRTBase(
5307 V[0] + GetPhys_Offset(i),
5308 tmp = outarray + GetCoeff_Offset(i));
5309 }
5310 }
5311 break;
5312 case 2:
5313 {
5314 for (int i = 0; i < GetExpSize(); ++i)
5315 {
5316 (*m_exp)[i]->NormVectorIProductWRTBase(
5317 V[0] + GetPhys_Offset(i), V[1] + GetPhys_Offset(i),
5318 tmp = outarray + GetCoeff_Offset(i));
5319 }
5320 }
5321 break;
5322 case 3:
5323 {
5324 for (int i = 0; i < GetExpSize(); ++i)
5325 {
5326 (*m_exp)[i]->NormVectorIProductWRTBase(
5327 V[0] + GetPhys_Offset(i), V[1] + GetPhys_Offset(i),
5328 V[2] + GetPhys_Offset(i),
5329 tmp = outarray + GetCoeff_Offset(i));
5330 }
5331 }
5332 break;
5333 default:
5334 NEKERROR(ErrorUtil::efatal, "Dimension not supported");
5335 break;
5336 }
5337}
5338
5340 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5341{
5343 "This method is not defined or valid for this class type");
5344}
5345
5346/**
5347 */
5349 [[maybe_unused]] const Array<OneD, NekDouble> coeffs)
5350{
5352 "This method is not defined or valid for this class type");
5353}
5354
5355/**
5356 */
5358 [[maybe_unused]] const int nreg,
5359 [[maybe_unused]] const Array<OneD, NekDouble> coeffs)
5360{
5362 "This method is not defined or valid for this class type");
5363}
5364
5365void ExpList::v_LocalToGlobal([[maybe_unused]] bool useComm)
5366{
5368 "This method is not defined or valid for this class type");
5369}
5370
5372 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5373 [[maybe_unused]] Array<OneD, NekDouble> &outarray,
5374 [[maybe_unused]] bool useComm)
5375{
5377 "This method is not defined or valid for this class type");
5378}
5379
5381{
5383 "This method is not defined or valid for this class type");
5384}
5385
5387 [[maybe_unused]] const Array<OneD, const NekDouble> &inarray,
5388 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
5389{
5391 "This method is not defined or valid for this class type");
5392}
5393
5395 Array<OneD, NekDouble> &outarray)
5396{
5397 v_FwdTransLocalElmt(inarray, outarray);
5398}
5399
5400/**
5401 * The operation is evaluated locally for every element by the function
5402 * StdRegions#StdExpansion#IProductWRTBase.
5403 *
5404 * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
5405 * containing the values of the function
5406 * \f$f(\boldsymbol{x})\f$ at the quadrature
5407 * points \f$\boldsymbol{x}_i\f$.
5408 * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
5409 * used to store the result.
5410 */
5412 Array<OneD, NekDouble> &outarray)
5413{
5414 LibUtilities::Timer timer;
5415 timer.Start();
5416 // initialise if required
5418 {
5419 for (int i = 0; i < m_collections.size(); ++i)
5420 {
5422 }
5424 }
5425
5427 int input_offset{0};
5428 int output_offset{0};
5429 for (int i = 0; i < m_collections.size(); ++i)
5430 {
5432 inarray + input_offset,
5433 tmp = outarray + output_offset);
5434 input_offset +=
5436 output_offset +=
5438 }
5439 timer.Stop();
5440 // Elapsed time
5441 timer.AccumulateRegion("Collections:IProductWRTBase", 10);
5442}
5443
5444/**
5445 * The operation is evaluated locally by the elemental
5446 * function StdRegions#StdExpansion#GetCoords.
5447 *
5448 * @param coord_0 After calculation, the \f$x_1\f$ coordinate
5449 * will be stored in this array.
5450 * @param coord_1 After calculation, the \f$x_2\f$ coordinate
5451 * will be stored in this array.
5452 * @param coord_2 After calculation, the \f$x_3\f$ coordinate
5453 * will be stored in this array.
5454 */
5456 Array<OneD, NekDouble> &coord_1,
5457 Array<OneD, NekDouble> &coord_2)
5458{
5459 if (GetNumElmts() == 0)
5460 {
5461 return;
5462 }
5463
5464 int i;
5465 Array<OneD, NekDouble> e_coord_0;
5466 Array<OneD, NekDouble> e_coord_1;
5467 Array<OneD, NekDouble> e_coord_2;
5468
5469 switch (GetExp(0)->GetCoordim())
5470 {
5471 case 1:
5472 for (i = 0; i < (*m_exp).size(); ++i)
5473 {
5474 e_coord_0 = coord_0 + m_phys_offset[i];
5475 (*m_exp)[i]->GetCoords(e_coord_0);
5476 }
5477 break;
5478 case 2:
5479 ASSERTL0(coord_1.size() != 0, "output coord_1 is not defined");
5480
5481 for (i = 0; i < (*m_exp).size(); ++i)
5482 {
5483 e_coord_0 = coord_0 + m_phys_offset[i];
5484 e_coord_1 = coord_1 + m_phys_offset[i];
5485 (*m_exp)[i]->GetCoords(e_coord_0, e_coord_1);
5486 }
5487 break;
5488 case 3:
5489 ASSERTL0(coord_1.size() != 0, "output coord_1 is not defined");
5490 ASSERTL0(coord_2.size() != 0, "output coord_2 is not defined");
5491
5492 for (i = 0; i < (*m_exp).size(); ++i)
5493 {
5494 e_coord_0 = coord_0 + m_phys_offset[i];
5495 e_coord_1 = coord_1 + m_phys_offset[i];
5496 e_coord_2 = coord_2 + m_phys_offset[i];
5497 (*m_exp)[i]->GetCoords(e_coord_0, e_coord_1, e_coord_2);
5498 }
5499 break;
5500 }
5501}
5502
5506{
5507 (*m_exp)[eid]->GetCoords(xc0, xc1, xc2);
5508}
5509
5510/**
5511 * @brief: Set up a normal along the trace elements between
5512 * two elements at elemental level
5513 *
5514 */
5516{
5517 for (int i = 0; i < m_exp->size(); ++i)
5518 {
5519 for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
5520 {
5521 (*m_exp)[i]->ComputeTraceNormal(j);
5522 }
5523 }
5524}
5525
5526/**
5527 */
5529 [[maybe_unused]] int i, [[maybe_unused]] std::shared_ptr<ExpList> &result,
5530 [[maybe_unused]] const bool DeclareCoeffPhysArrays)
5531{
5533 "This method is not defined or valid for this class type");
5534}
5535
5536/**
5537 */
5539 const Array<OneD, NekDouble> &element,
5540 Array<OneD, NekDouble> &boundary)
5541{
5542 int n, cnt;
5543 Array<OneD, NekDouble> tmp1, tmp2;
5545
5546 Array<OneD, int> ElmtID, EdgeID;
5547 GetBoundaryToElmtMap(ElmtID, EdgeID);
5548
5549 // Initialise result
5550 boundary =
5552
5553 // Skip other boundary regions
5554 for (cnt = n = 0; n < i; ++n)
5555 {
5556 cnt += GetBndCondExpansions()[n]->GetExpSize();
5557 }
5558
5559 int offsetBnd;
5560 int offsetElmt = 0;
5561 for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
5562 {
5563 offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
5564
5565 elmt = GetExp(ElmtID[cnt + n]);
5566 elmt->GetTracePhysVals(
5567 EdgeID[cnt + n], GetBndCondExpansions()[i]->GetExp(n),
5568 tmp1 = element + offsetElmt, tmp2 = boundary + offsetBnd);
5569
5570 offsetElmt += elmt->GetTotPoints();
5571 }
5572}
5573
5574/**
5575 */
5577 const Array<OneD, const NekDouble> &phys,
5578 Array<OneD, NekDouble> &bndElmt)
5579{
5580 int n, cnt, nq;
5581
5582 Array<OneD, int> ElmtID, EdgeID;
5583 GetBoundaryToElmtMap(ElmtID, EdgeID);
5584
5585 // Skip other boundary regions
5586 for (cnt = n = 0; n < i; ++n)
5587 {
5588 cnt += GetBndCondExpansions()[n]->GetExpSize();
5589 }
5590
5591 // Count number of points
5592 int npoints = 0;
5593 for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
5594 {
5595 npoints += GetExp(ElmtID[cnt + n])->GetTotPoints();
5596 }
5597
5598 // Initialise result
5599 bndElmt = Array<OneD, NekDouble>(npoints, 0.0);
5600
5601 // Extract data
5602 int offsetPhys;
5603 int offsetElmt = 0;
5604 for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
5605 {
5606 nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
5607 offsetPhys = GetPhys_Offset(ElmtID[cnt + n]);
5608 Vmath::Vcopy(nq, &phys[offsetPhys], 1, &bndElmt[offsetElmt], 1);
5609 offsetElmt += nq;
5610 }
5611}
5612
5613/**
5614 */
5616 const Array<OneD, const NekDouble> &phys,
5618{
5619 int n, cnt;
5622
5623 Array<OneD, int> ElmtID, EdgeID;
5624 GetBoundaryToElmtMap(ElmtID, EdgeID);
5625
5626 // Initialise result
5627 bnd =
5629
5630 // Skip other boundary regions
5631 for (cnt = n = 0; n < i; ++n)
5632 {
5633 cnt += GetBndCondExpansions()[n]->GetExpSize();
5634 }
5635
5636 int offsetBnd;
5637 int offsetPhys;
5638 for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
5639 {
5640 offsetPhys = GetPhys_Offset(ElmtID[cnt + n]);
5641 offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
5642
5643 elmt = GetExp(ElmtID[cnt + n]);
5644 elmt->GetTracePhysVals(EdgeID[cnt + n],
5646 phys + offsetPhys, tmp1 = bnd + offsetBnd);
5647 }
5648}
5649
5650/**
5651 */
5654{
5655 int j, n, cnt;
5656 int coordim = GetCoordim(0);
5659
5660 Array<OneD, int> ElmtID, EdgeID;
5661 GetBoundaryToElmtMap(ElmtID, EdgeID);
5662
5663 // Initialise result
5664 normals = Array<OneD, Array<OneD, NekDouble>>(coordim);
5665 for (j = 0; j < coordim; ++j)
5666 {
5667 normals[j] = Array<OneD, NekDouble>(
5668 GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
5669 }
5670
5671 // Skip other boundary regions
5672 for (cnt = n = 0; n < i; ++n)
5673 {
5674 cnt += GetBndCondExpansions()[n]->GetExpSize();
5675 }
5676
5677 int offset;
5678 for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
5679 {
5680 offset = GetBndCondExpansions()[i]->GetPhys_Offset(n);
5681
5682 elmt = GetExp(ElmtID[cnt + n]);
5683 const Array<OneD, const Array<OneD, NekDouble>> normalsElmt =
5684 elmt->GetTraceNormal(EdgeID[cnt + n]);
5685
5686 // Interp/Copy to result
5687 for (j = 0; j < coordim; ++j)
5688 {
5689 GetBndCondExpansions()[i]->GetExp(n)->PhysInterp(
5690 elmt, normalsElmt[j], tmp = normals[j] + offset);
5691 }
5692 }
5693}
5694
5695/**
5696 */
5698 [[maybe_unused]] Array<OneD, int> &EdgeID)
5699{
5701 "This method is not defined or valid for this class type");
5702}
5703
5705 [[maybe_unused]] Array<OneD, NekDouble> &weightave,
5706 [[maybe_unused]] Array<OneD, NekDouble> &weightjmp)
5707{
5708 NEKERROR(ErrorUtil::efatal, "v_FillBwdWithBwdWeight not defined");
5709}
5710
5712 [[maybe_unused]] const Array<OneD, const NekDouble> &Fwd,
5713 [[maybe_unused]] Array<OneD, NekDouble> &Bwd)
5714{
5715 NEKERROR(ErrorUtil::efatal, "v_PeriodicBwdCopy not defined");
5716}
5717
5718/**
5719 */
5722{
5724 "This method is not defined or valid for this class type");
5726 return result;
5727}
5728
5729/**
5730 */
5733{
5735 "This method is not defined or valid for this class type");
5737 return result;
5738}
5739
5740/**
5741 */
5743 [[maybe_unused]] const NekDouble time,
5744 [[maybe_unused]] const std::string varName,
5745 [[maybe_unused]] const NekDouble x2_in,
5746 [[maybe_unused]] const NekDouble x3_in)
5747{
5749 "This method is not defined or valid for this class type");
5750}
5751
5752/**
5753 */
5754map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
5755{
5757 "This method is not defined or valid for this class type");
5758 static map<int, RobinBCInfoSharedPtr> result;
5759 return result;
5760}
5761
5762/**
5763 */
5764void ExpList::v_GetPeriodicEntities([[maybe_unused]] PeriodicMap &periodicVerts,
5765 [[maybe_unused]] PeriodicMap &periodicEdges,
5766 [[maybe_unused]] PeriodicMap &periodicFaces)
5767{
5769 "This method is not defined or valid for this class type");
5770}
5771
5774 unsigned int regionId, const std::string &variable)
5775{
5776 auto collectionIter = collection.find(regionId);
5777 ASSERTL1(collectionIter != collection.end(),
5778 "Unable to locate collection " + std::to_string(regionId));
5779
5781 (*collectionIter).second;
5782 auto conditionMapIter = bndCondMap->find(variable);
5783 ASSERTL1(conditionMapIter != bndCondMap->end(),
5784 "Unable to locate condition map.");
5785
5786 const SpatialDomains::BoundaryConditionShPtr boundaryCondition =
5787 (*conditionMapIter).second;
5788
5789 return boundaryCondition;
5790}
5791
5793{
5795 "This method is not defined or valid for this class type");
5796 return NullExpListSharedPtr;
5797}
5798
5799/**
5800 * @brief Construct collections of elements containing a single element
5801 * type and polynomial order from the list of expansions.
5802 */
5804{
5805 // Set up initialisation flags
5807 std::vector<bool>(Collections::SIZE_OperatorType, true);
5808
5809 // Figure out optimisation parameters if provided in
5810 // session file or default given
5812 m_session, (*m_exp)[0]->GetShapeDimension(), ImpType);
5813
5814 // turn on autotuning if explicitly specified in xml file
5815 // or command line option is set but only do optimisation
5816 // for volumetric elements (not boundary condition)
5817 bool autotuning = colOpt.IsUsingAutotuning();
5818 if ((autotuning == false) && (ImpType == Collections::eNoImpType))
5819 {
5820 // turn on autotuning if write-opt-file specified
5821 // if m_graph available
5822 if (m_session->GetUpdateOptFile() && m_graph)
5823 {
5824 // only turn on autotuning for volumetric elements
5825 // where Mesh Dimension is equal to the Shape
5826 // Dimension of element.
5827 if (m_graph->GetMeshDimension() == (*m_exp)[0]->GetShapeDimension())
5828 {
5829 autotuning = true;
5830 }
5831 }
5832 }
5833 bool verbose = (m_session->DefinesCmdLineArgument("verbose")) &&
5834 (m_session->GetComm()->GetRank() == 0);
5835
5836 // clear vectors in case previously called
5837 m_collections.clear();
5838
5839 /*-------------------------------------------------------------------------
5840 Dividing m_exp into sub groups (collections): original exp order is kept.
5841 Use 3 basiskey + deformed flag to determine if two exp belong to the same
5842 collection or not.
5843 -------------------------------------------------------------------------*/
5844 // the maximum size is either explicitly specified, or set to very large
5845 // value which will not affect the actual collection size
5846 int collmax =
5847 (colOpt.GetMaxCollectionSize() > 0 ? colOpt.GetMaxCollectionSize()
5848 : 2 * m_exp->size());
5849
5850 vector<StdRegions::StdExpansionSharedPtr> collExp;
5851 LocalRegions::ExpansionSharedPtr exp = (*m_exp)[0];
5852 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(exp);
5853
5854 // add the first element to the collection - initialization
5855 collExp.push_back(exp);
5856 // collcnt is the number of elements in current collection
5857 int collcnt = 1;
5858 // initialize the basisKeys to NullBasisKey
5859 std::vector<LibUtilities::BasisKey> thisbasisKeys(
5861 std::vector<LibUtilities::BasisKey> prevbasisKeys(
5863 // fetch basiskeys of the first element
5864 for (int d = 0; d < exp->GetNumBases(); d++)
5865 {
5866 prevbasisKeys[d] = exp->GetBasis(d)->GetBasisKey();
5867 }
5868
5869 // initialize the deformed flag based on the first element
5870 bool prevDef =
5871 exp->GetMetricInfo()->GetGtype() == SpatialDomains::eDeformed;
5872 // collsize is the maximum size among all collections
5873 int collsize = 0;
5874 int mincol = (*m_exp).size();
5875 int maxcol = -1;
5876 int meancol = 0;
5877
5878 for (int i = 1; i < (*m_exp).size(); i++)
5879 {
5880 exp = (*m_exp)[i];
5881
5882 // fetch basiskeys of current element
5883 for (int d = 0; d < exp->GetNumBases(); d++)
5884 {
5885 thisbasisKeys[d] = exp->GetBasis(d)->GetBasisKey();
5886 }
5887 // fetch deformed flag of current element
5888 bool Deformed =
5889 (exp->GetMetricInfo()->GetGtype() == SpatialDomains::eDeformed);
5890
5891 // Check if this element is the same type as the previous one or
5892 // if we have reached the maximum collection size
5893 if (thisbasisKeys != prevbasisKeys || prevDef != Deformed ||
5894 collcnt >= collmax)
5895 {
5896 // if no Imp Type provided and No
5897 // setting in xml file. reset
5898 // impTypes using timings
5899 if (autotuning)
5900 {
5901 // if current collection is larger than previous one
5902 // update impTypes; otherwise, use previous impTypes
5903 if (collExp.size() > collsize)
5904 {
5905 impTypes =
5906 colOpt.SetWithTimings(collExp, impTypes, verbose);
5907 collsize = collExp.size();
5908 }
5909 }
5911 colOpt.GetOperatorImpMap(exp);
5912
5913 Collections::Collection tmp(collExp, impTypes);
5914 m_collections.push_back(tmp);
5915 mincol = min(mincol, (int)collExp.size());
5916 maxcol = max(maxcol, (int)collExp.size());
5917 meancol += collExp.size();
5918
5919 // for the new collection calling the optimization routine based on
5920 // its first element
5921 impTypes = colOpt.GetOperatorImpMap((*m_exp)[i]);
5922
5923 // clean-up current element list - temporary collection
5924 collExp.clear();
5925 collcnt = 0;
5926 collsize = 0;
5927 }
5928
5929 // insert exp and increment count
5930 collExp.push_back(exp);
5931 collcnt++;
5932 // update previous info
5933 prevbasisKeys = thisbasisKeys;
5934 prevDef = Deformed;
5935 }
5936
5937 // execute autotuning for the last collection
5938 if (autotuning)
5939 {
5940 if (collExp.size() > collsize)
5941 {
5942 impTypes = colOpt.SetWithTimings(collExp, impTypes, verbose);
5943 collsize = collExp.size();
5944 }
5945 }
5946 Collections::Collection tmp(collExp, impTypes);
5947 m_collections.push_back(tmp);
5948 if (verbose)
5949 {
5950 mincol = min(mincol, (int)collExp.size());
5951 maxcol = max(maxcol, (int)collExp.size());
5952 meancol += collExp.size();
5953 meancol /= m_collections.size();
5954 cout << "Collection group: num. = " << m_collections.size()
5955 << "; mean len = " << meancol << " (min = " << mincol
5956 << ", max = " << maxcol << ")" << endl;
5957 }
5958
5959 // clean-up current element list - temporary collection
5960 collExp.clear();
5961 collcnt = 0;
5962 collsize = 0;
5963
5964 // update optimisation file
5965 if ((m_session->GetUpdateOptFile()) && (ImpType == Collections::eNoImpType))
5966 {
5967 colOpt.UpdateOptFile(m_session->GetSessionName(), m_comm);
5968 // turn off write-opt-file option so only first
5969 // instance is timed
5970 m_session->SetUpdateOptFile(false);
5971 }
5972}
5973
5975{
5977}
5978
5979/**
5980 * Added for access to the pool count by external code (e.g. UnitTests)
5981 * which can't access the static pool across compilation units on
5982 * Windows builds.
5983 */
5984int ExpList::GetPoolCount(std::string poolName)
5985{
5986 return v_GetPoolCount(poolName);
5987}
5988
5989void ExpList::UnsetGlobalLinSys(GlobalLinSysKey key, bool clearLocalMatrices)
5990{
5991 v_UnsetGlobalLinSys(key, clearLocalMatrices);
5992}
5993
5996{
5997 return v_GetGlobalLinSysManager();
5998}
5999
6000void ExpList::v_PhysInterp1DScaled([[maybe_unused]] const NekDouble scale,
6001 const Array<OneD, NekDouble> &inarray,
6002 Array<OneD, NekDouble> &outarray)
6003{
6004 // the scaling factor for the PhysInterp1DScaled is given as NekDouble
6005 // however inside Collections it is treated as a FactorMap
6006 // defining needed FactorMap to pass the scaling factor as an input to
6007 // Collections
6009 // Updating the FactorMap according to the scale input
6011 LibUtilities::Timer timer;
6012
6013 // initialise if required
6014 if (m_collections.size() &&
6016 {
6017 for (int i = 0; i < m_collections.size(); ++i)
6018 {
6021 }
6022 }
6023 // once the collections are initialized, check for the scaling factor
6024 for (int i = 0; i < m_collections.size(); ++i)
6025
6026 {
6028 factors);
6029 }
6030
6031 LIKWID_MARKER_START("v_PhysInterp1DScaled");
6032 timer.Start();
6034 int input_offset{0};
6035 int output_offset{0};
6036 for (int i = 0; i < m_collections.size(); ++i)
6037 {
6039 inarray + input_offset,
6040 tmp = outarray + output_offset);
6041 input_offset +=
6043 output_offset +=
6045 }
6046 timer.Stop();
6047 LIKWID_MARKER_STOP("v_PhysInterp1DScaled");
6048 timer.AccumulateRegion("Collections:PhysInterp1DScaled", 10);
6049}
6051 [[maybe_unused]] const Array<OneD, const NekDouble> &FwdFlux,
6052 [[maybe_unused]] const Array<OneD, const NekDouble> &BwdFlux,
6053 [[maybe_unused]] Array<OneD, NekDouble> &outarray)
6054{
6055 NEKERROR(ErrorUtil::efatal, "AddTraceIntegralToOffDiag not defined");
6056}
6057
6059 const Array<OneD, const Array<OneD, NekDouble>> &inarray, const int nDirctn,
6061{
6062 int nTotElmt = (*m_exp).size();
6063 int nElmtPnt = (*m_exp)[0]->GetTotPoints();
6064 int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
6065
6066 Array<OneD, NekDouble> tmpCoef(nElmtCoef, 0.0);
6067 Array<OneD, NekDouble> tmpPhys(nElmtPnt, 0.0);
6068
6069 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
6070 {
6071 nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
6072 nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
6073
6074 if (tmpPhys.size() != nElmtPnt || tmpCoef.size() != nElmtCoef)
6075 {
6076 tmpPhys = Array<OneD, NekDouble>(nElmtPnt, 0.0);
6077 tmpCoef = Array<OneD, NekDouble>(nElmtCoef, 0.0);
6078 }
6079
6080 for (int ncl = 0; ncl < nElmtPnt; ncl++)
6081 {
6082 tmpPhys[ncl] = inarray[nelmt][ncl];
6083
6084 (*m_exp)[nelmt]->IProductWRTDerivBase(nDirctn, tmpPhys, tmpCoef);
6085
6086 for (int nrw = 0; nrw < nElmtCoef; nrw++)
6087 {
6088 (*mtxPerVar[nelmt])(nrw, ncl) = tmpCoef[nrw];
6089 }
6090 // to maintain all the other columes are zero.
6091 tmpPhys[ncl] = 0.0;
6092 }
6093 }
6094}
6095
6098{
6099 int nTotElmt = (*m_exp).size();
6100
6101 int nspacedim = m_graph->GetSpaceDimension();
6102 Array<OneD, Array<OneD, NekDouble>> projectedpnts(nspacedim);
6103 Array<OneD, Array<OneD, NekDouble>> tmppnts(nspacedim);
6104 Array<OneD, DNekMatSharedPtr> ArrayStdMat(nspacedim);
6105 Array<OneD, Array<OneD, NekDouble>> ArrayStdMat_data(nspacedim);
6106
6107 Array<OneD, NekDouble> clmnArray, clmnStdMatArray;
6108
6110 int nElmtPntPrevious = 0;
6111 int nElmtCoefPrevious = 0;
6112
6113 int nElmtPnt, nElmtCoef;
6114 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
6115 {
6116 nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
6117 nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
6118 LibUtilities::ShapeType ElmtTypeNow = (*m_exp)[nelmt]->DetShapeType();
6119
6120 if (nElmtPntPrevious != nElmtPnt || nElmtCoefPrevious != nElmtCoef ||
6121 (ElmtTypeNow != ElmtTypePrevious))
6122 {
6123 if (nElmtPntPrevious != nElmtPnt)
6124 {
6125 for (int ndir = 0; ndir < nspacedim; ndir++)
6126 {
6127 projectedpnts[ndir] = Array<OneD, NekDouble>(nElmtPnt, 0.0);
6128 tmppnts[ndir] = Array<OneD, NekDouble>(nElmtPnt, 0.0);
6129 }
6130 }
6132 stdExp = (*m_exp)[nelmt]->GetStdExp();
6134 stdExp->DetShapeType(), *stdExp);
6135
6136 ArrayStdMat[0] = stdExp->GetStdMatrix(matkey);
6137 ArrayStdMat_data[0] = ArrayStdMat[0]->GetPtr();
6138
6139 if (nspacedim > 1)
6140 {
6142 StdRegions::eDerivBase1, stdExp->DetShapeType(), *stdExp);
6143
6144 ArrayStdMat[1] = stdExp->GetStdMatrix(matkey);
6145 ArrayStdMat_data[1] = ArrayStdMat[1]->GetPtr();
6146
6147 if (nspacedim > 2)
6148 {
6150 stdExp->DetShapeType(),
6151 *stdExp);
6152
6153 ArrayStdMat[2] = stdExp->GetStdMatrix(matkey);
6154 ArrayStdMat_data[2] = ArrayStdMat[2]->GetPtr();
6155 }
6156 }
6157
6158 ElmtTypePrevious = ElmtTypeNow;
6159 nElmtPntPrevious = nElmtPnt;
6160 nElmtCoefPrevious = nElmtCoef;
6161 }
6162 else
6163 {
6164 for (int ndir = 0; ndir < nspacedim; ndir++)
6165 {
6166 Vmath::Zero(nElmtPnt, projectedpnts[ndir], 1);
6167 }
6168 }
6169
6170 for (int ndir = 0; ndir < nspacedim; ndir++)
6171 {
6172 (*m_exp)[nelmt]->AlignVectorToCollapsedDir(
6173 ndir, inarray[ndir][nelmt], tmppnts);
6174 for (int n = 0; n < nspacedim; n++)
6175 {
6176 Vmath::Vadd(nElmtPnt, tmppnts[n], 1, projectedpnts[n], 1,
6177 projectedpnts[n], 1);
6178 }
6179 }
6180
6181 for (int ndir = 0; ndir < nspacedim; ndir++)
6182 {
6183 // weight with metric
6184 (*m_exp)[nelmt]->MultiplyByQuadratureMetric(projectedpnts[ndir],
6185 projectedpnts[ndir]);
6186 Array<OneD, NekDouble> MatDataArray = mtxPerVar[nelmt]->GetPtr();
6187
6188 for (int np = 0; np < nElmtPnt; np++)
6189 {
6190 NekDouble factor = projectedpnts[ndir][np];
6191 clmnArray = MatDataArray + np * nElmtCoef;
6192 clmnStdMatArray = ArrayStdMat_data[ndir] + np * nElmtCoef;
6193 Vmath::Svtvp(nElmtCoef, factor, clmnStdMatArray, 1, clmnArray,
6194 1, clmnArray, 1);
6195 }
6196 }
6197 }
6198}
6199
6200// TODO: Reduce cost by getting Bwd Matrix directly
6202 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
6204{
6206 int nElmtPntPrevious = 0;
6207 int nElmtCoefPrevious = 0;
6208 int nTotElmt = (*m_exp).size();
6209 int nElmtPnt = (*m_exp)[0]->GetTotPoints();
6210 int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
6211
6212 Array<OneD, NekDouble> tmpPhys;
6213 Array<OneD, NekDouble> clmnArray, clmnStdMatArray;
6214 Array<OneD, NekDouble> stdMat_data;
6215
6216 for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
6217 {
6218 nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
6219 nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
6220 LibUtilities::ShapeType ElmtTypeNow = (*m_exp)[nelmt]->DetShapeType();
6221
6222 if (nElmtPntPrevious != nElmtPnt || nElmtCoefPrevious != nElmtCoef ||
6223 (ElmtTypeNow != ElmtTypePrevious))
6224 {
6226 stdExp = (*m_exp)[nelmt]->GetStdExp();
6228 stdExp->DetShapeType(), *stdExp);
6229
6230 DNekMatSharedPtr BwdMat = stdExp->GetStdMatrix(matkey);
6231 stdMat_data = BwdMat->GetPtr();
6232
6233 if (nElmtPntPrevious != nElmtPnt)
6234 {
6235 tmpPhys = Array<OneD, NekDouble>(nElmtPnt, 0.0);
6236 }
6237
6238 ElmtTypePrevious = ElmtTypeNow;
6239 nElmtPntPrevious = nElmtPnt;
6240 nElmtCoefPrevious = nElmtCoef;
6241 }
6242
6243 (*m_exp)[nelmt]->MultiplyByQuadratureMetric(
6244 inarray[nelmt],
6245 tmpPhys); // weight with metric
6246
6247 Array<OneD, NekDouble> MatDataArray = mtxPerVar[nelmt]->GetPtr();
6248
6249 for (int np = 0; np < nElmtPnt; np++)
6250 {
6251 NekDouble factor = tmpPhys[np];
6252 clmnArray = MatDataArray + np * nElmtCoef;
6253 clmnStdMatArray = stdMat_data + np * nElmtCoef;
6254 Vmath::Smul(nElmtCoef, factor, clmnStdMatArray, 1, clmnArray, 1);
6255 }
6256 }
6257}
6258
6259/**
6260 * @brief inverse process of v_GetFwdBwdTracePhys. Given Trace integration of
6261 * Fwd and Bwd Jacobian, with dimension NtotalTrace*TraceCoef*TracePhys.
6262 * return Elemental Jacobian matrix with dimension
6263 * NtotalElement*ElementCoef*ElementPhys.
6264 *
6265 *
6266 * @param field is a NekDouble array which contains the 2D data
6267 * from which we wish to extract the backward and
6268 * forward orientated trace/edge arrays.
6269 * @param Fwd The resulting forwards space.
6270 * @param Bwd The resulting backwards space.
6271 */
6276{
6278 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
6279 tracelist->GetExp();
6280 int ntotTrace = (*traceExp).size();
6281 int nTracePnt, nTraceCoef;
6282
6283 std::shared_ptr<LocalRegions::ExpansionVector> fieldExp = GetExp();
6284 int nElmtCoef;
6285
6286 const MultiRegions::LocTraceToTraceMapSharedPtr locTraceToTraceMap =
6288 const Array<OneD, const Array<OneD, int>> LRAdjExpid =
6289 locTraceToTraceMap->GetLeftRightAdjacentExpId();
6290 const Array<OneD, const Array<OneD, bool>> LRAdjflag =
6291 locTraceToTraceMap->GetLeftRightAdjacentExpFlag();
6292
6294 locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffMap();
6296 locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffSign();
6297 DNekMatSharedPtr ElmtMat;
6298 Array<OneD, NekDouble> ElmtMat_data;
6299 // int nclAdjExp;
6300 int nrwAdjExp;
6301 int MatIndex, nPnts;
6302 NekDouble sign = 1.0;
6303
6304 int nTracePntsTtl = tracelist->GetTotPoints();
6305 int nlocTracePts = locTraceToTraceMap->GetNLocTracePts();
6306 int nlocTracePtsFwd = locTraceToTraceMap->GetNFwdLocTracePts();
6307 int nlocTracePtsBwd = nlocTracePts - nlocTracePtsFwd;
6308
6309 Array<OneD, int> nlocTracePtsLR(2);
6310 nlocTracePtsLR[0] = nlocTracePtsFwd;
6311 nlocTracePtsLR[1] = nlocTracePtsBwd;
6312
6313 size_t nFwdBwdNonZero = 0;
6314 Array<OneD, int> tmpIndex{2, -1};
6315 for (int i = 0; i < 2; ++i)
6316 {
6317 if (nlocTracePtsLR[i] > 0)
6318 {
6319 tmpIndex[nFwdBwdNonZero] = i;
6320 nFwdBwdNonZero++;
6321 }
6322 }
6323
6324 Array<OneD, int> nlocTracePtsNonZeroIndex{nFwdBwdNonZero};
6325 for (int i = 0; i < nFwdBwdNonZero; ++i)
6326 {
6327 nlocTracePtsNonZeroIndex[i] = tmpIndex[i];
6328 }
6329
6330 Array<OneD, NekDouble> TraceFwdPhy(nTracePntsTtl);
6331 Array<OneD, NekDouble> TraceBwdPhy(nTracePntsTtl);
6332 Array<OneD, Array<OneD, NekDouble>> tmplocTrace(2);
6333 for (int k = 0; k < 2; ++k)
6334 {
6335 tmplocTrace[k] = NullNekDouble1DArray;
6336 }
6337
6338 for (int k = 0; k < nFwdBwdNonZero; ++k)
6339 {
6340 size_t i = nlocTracePtsNonZeroIndex[k];
6341 tmplocTrace[i] = Array<OneD, NekDouble>(nlocTracePtsLR[i]);
6342 }
6343
6344 int nNumbElmt = fieldMat.size();
6345 Array<OneD, Array<OneD, NekDouble>> ElmtMatDataArray(nNumbElmt);
6346 Array<OneD, int> ElmtCoefArray(nNumbElmt);
6347 for (int i = 0; i < nNumbElmt; i++)
6348 {
6349 ElmtMatDataArray[i] = fieldMat[i]->GetPtr();
6350 ElmtCoefArray[i] = GetNcoeffs(i);
6351 }
6352
6353 int nTraceCoefMax = 0;
6354 int nTraceCoefMin = std::numeric_limits<int>::max();
6355 Array<OneD, int> TraceCoefArray(ntotTrace);
6356 Array<OneD, int> TracePntArray(ntotTrace);
6357 Array<OneD, int> TraceOffArray(ntotTrace);
6358 Array<OneD, Array<OneD, NekDouble>> FwdMatData(ntotTrace);
6359 Array<OneD, Array<OneD, NekDouble>> BwdMatData(ntotTrace);
6360 for (int nt = 0; nt < ntotTrace; nt++)
6361 {
6362 nTraceCoef = (*traceExp)[nt]->GetNcoeffs();
6363 nTracePnt = tracelist->GetTotPoints(nt);
6364 int noffset = tracelist->GetPhys_Offset(nt);
6365 TraceCoefArray[nt] = nTraceCoef;
6366 TracePntArray[nt] = nTracePnt;
6367 TraceOffArray[nt] = noffset;
6368 FwdMatData[nt] = FwdMat[nt]->GetPtr();
6369 BwdMatData[nt] = BwdMat[nt]->GetPtr();
6370 if (nTraceCoef > nTraceCoefMax)
6371 {
6372 nTraceCoefMax = nTraceCoef;
6373 }
6374 if (nTraceCoef < nTraceCoefMin)
6375 {
6376 nTraceCoefMin = nTraceCoef;
6377 }
6378 }
6379 WARNINGL1(nTraceCoefMax == nTraceCoefMin,
6380 "nTraceCoefMax!=nTraceCoefMin: Effeciency may be low ");
6381
6382 int traceID, nfieldPnts, ElmtId, noffset;
6383 const Array<OneD, const Array<OneD, int>> LocTracephysToTraceIDMap =
6384 locTraceToTraceMap->GetLocTracephysToTraceIDMap();
6385 const Array<OneD, const int> fieldToLocTraceMap =
6386 locTraceToTraceMap->GetLocTraceToFieldMap();
6387 Array<OneD, Array<OneD, int>> fieldToLocTraceMapLR(2);
6388 noffset = 0;
6389 for (int k = 0; k < nFwdBwdNonZero; ++k)
6390 {
6391 size_t i = nlocTracePtsNonZeroIndex[k];
6392 fieldToLocTraceMapLR[i] = Array<OneD, int>(nlocTracePtsLR[i]);
6393 Vmath::Vcopy(nlocTracePtsLR[i], &fieldToLocTraceMap[0] + noffset, 1,
6394 &fieldToLocTraceMapLR[i][0], 1);
6395 noffset += nlocTracePtsLR[i];
6396 }
6397
6398 Array<OneD, Array<OneD, int>> MatIndexArray(2);
6399 for (int k = 0; k < nFwdBwdNonZero; ++k)
6400 {
6401 size_t nlr = nlocTracePtsNonZeroIndex[k];
6402 MatIndexArray[nlr] = Array<OneD, int>(nlocTracePtsLR[nlr]);
6403 for (int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
6404 {
6405 traceID = LocTracephysToTraceIDMap[nlr][nloc];
6406 nTraceCoef = TraceCoefArray[traceID];
6407 ElmtId = LRAdjExpid[nlr][traceID];
6408 noffset = GetPhys_Offset(ElmtId);
6409 nElmtCoef = ElmtCoefArray[ElmtId];
6410 nfieldPnts = fieldToLocTraceMapLR[nlr][nloc];
6411 nPnts = nfieldPnts - noffset;
6412
6413 MatIndexArray[nlr][nloc] = nPnts * nElmtCoef;
6414 }
6415 }
6416
6417 for (int nc = 0; nc < nTraceCoefMin; nc++)
6418 {
6419 for (int nt = 0; nt < ntotTrace; nt++)
6420 {
6421 nTraceCoef = TraceCoefArray[nt];
6422 nTracePnt = TracePntArray[nt];
6423 noffset = TraceOffArray[nt];
6424 Vmath::Vcopy(nTracePnt, &FwdMatData[nt][nc], nTraceCoef,
6425 &TraceFwdPhy[noffset], 1);
6426 Vmath::Vcopy(nTracePnt, &BwdMatData[nt][nc], nTraceCoef,
6427 &TraceBwdPhy[noffset], 1);
6428 }
6429
6430 for (int k = 0; k < nFwdBwdNonZero; ++k)
6431 {
6432 size_t i = nlocTracePtsNonZeroIndex[k];
6433 Vmath::Zero(nlocTracePtsLR[i], tmplocTrace[i], 1);
6434 }
6435
6436 GetLocTraceFromTracePts(TraceFwdPhy, TraceBwdPhy, tmplocTrace[0],
6437 tmplocTrace[1]);
6438
6439 for (int k = 0; k < nFwdBwdNonZero; ++k)
6440 {
6441 size_t nlr = nlocTracePtsNonZeroIndex[k];
6442 for (int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
6443 {
6444 traceID = LocTracephysToTraceIDMap[nlr][nloc];
6445 nTraceCoef = TraceCoefArray[traceID];
6446 ElmtId = LRAdjExpid[nlr][traceID];
6447 nrwAdjExp = elmtLRMap[nlr][traceID][nc];
6448 sign = elmtLRSign[nlr][traceID][nc];
6449 MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
6450
6451 ElmtMatDataArray[ElmtId][MatIndex] -=
6452 sign * tmplocTrace[nlr][nloc];
6453 }
6454 }
6455 }
6456
6457 for (int nc = nTraceCoefMin; nc < nTraceCoefMax; nc++)
6458 {
6459 for (int nt = 0; nt < ntotTrace; nt++)
6460 {
6461 nTraceCoef = TraceCoefArray[nt];
6462 nTracePnt = TracePntArray[nt];
6463 noffset = TraceOffArray[nt];
6464 if (nc < nTraceCoef)
6465 {
6466 Vmath::Vcopy(nTracePnt, &FwdMatData[nt][nc], nTraceCoef,
6467 &TraceFwdPhy[noffset], 1);
6468 Vmath::Vcopy(nTracePnt, &BwdMatData[nt][nc], nTraceCoef,
6469 &TraceBwdPhy[noffset], 1);
6470 }
6471 else
6472 {
6473 Vmath::Zero(nTracePnt, &TraceFwdPhy[noffset], 1);
6474 Vmath::Zero(nTracePnt, &TraceBwdPhy[noffset], 1);
6475 }
6476 }
6477
6478 for (int k = 0; k < nFwdBwdNonZero; ++k)
6479 {
6480 size_t i = nlocTracePtsNonZeroIndex[k];
6481 Vmath::Zero(nlocTracePtsLR[i], tmplocTrace[i], 1);
6482 }
6483 GetLocTraceFromTracePts(TraceFwdPhy, TraceBwdPhy, tmplocTrace[0],
6484 tmplocTrace[1]);
6485
6486 for (int k = 0; k < nFwdBwdNonZero; ++k)
6487 {
6488 size_t nlr = nlocTracePtsNonZeroIndex[k];
6489 for (int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
6490 {
6491 traceID = LocTracephysToTraceIDMap[nlr][nloc];
6492 nTraceCoef = TraceCoefArray[traceID];
6493 if (nc < nTraceCoef)
6494 {
6495 ElmtId = LRAdjExpid[nlr][traceID];
6496 nrwAdjExp = elmtLRMap[nlr][traceID][nc];
6497 sign = -elmtLRSign[nlr][traceID][nc];
6498 MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
6499
6500 ElmtMatDataArray[ElmtId][MatIndex] +=
6501 sign * tmplocTrace[nlr][nloc];
6502 }
6503 }
6504 }
6505 }
6506}
6507
6509 const int dir, const Array<OneD, const DNekMatSharedPtr> ElmtJacQuad,
6511{
6512 int nelmt;
6513 int nelmtcoef, nelmtpnts, nelmtcoef0, nelmtpnts0;
6514
6515 nelmtcoef = GetNcoeffs(0);
6516 nelmtpnts = GetTotPoints(0);
6517
6518 Array<OneD, NekDouble> innarray(nelmtpnts, 0.0);
6519 Array<OneD, NekDouble> outarray(nelmtcoef, 0.0);
6520
6521 Array<OneD, NekDouble> MatQ_data;
6522 Array<OneD, NekDouble> MatC_data;
6523
6524 DNekMatSharedPtr tmpMatQ, tmpMatC;
6525
6526 nelmtcoef0 = nelmtcoef;
6527 nelmtpnts0 = nelmtpnts;
6528
6529 for (nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
6530 {
6531 nelmtcoef = GetNcoeffs(nelmt);
6532 nelmtpnts = GetTotPoints(nelmt);
6533
6534 tmpMatQ = ElmtJacQuad[nelmt];
6535 tmpMatC = ElmtJacCoef[nelmt];
6536
6537 MatQ_data = tmpMatQ->GetPtr();
6538 MatC_data = tmpMatC->GetPtr();
6539
6540 if (nelmtcoef != nelmtcoef0)
6541 {
6542 outarray = Array<OneD, NekDouble>(nelmtcoef, 0.0);
6543 nelmtcoef0 = nelmtcoef;
6544 }
6545
6546 if (nelmtpnts != nelmtpnts0)
6547 {
6548 innarray = Array<OneD, NekDouble>(nelmtpnts, 0.0);
6549 nelmtpnts0 = nelmtpnts;
6550 }
6551
6552 for (int np = 0; np < nelmtcoef; np++)
6553 {
6554 Vmath::Vcopy(nelmtpnts, &MatQ_data[0] + np, nelmtcoef, &innarray[0],
6555 1);
6556 (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray, innarray);
6557 (*m_exp)[nelmt]->IProductWRTDerivBase(dir, innarray, outarray);
6558
6559 Vmath::Vadd(nelmtcoef, &outarray[0], 1, &MatC_data[0] + np,
6560 nelmtcoef, &MatC_data[0] + np, nelmtcoef);
6561 }
6562 }
6563}
6564
6566 const Array<OneD, const DNekMatSharedPtr> ElmtJacQuad,
6568{
6569 int nelmt;
6570 int nelmtcoef, nelmtpnts, nelmtcoef0, nelmtpnts0;
6571
6572 nelmtcoef = GetNcoeffs(0);
6573 nelmtpnts = GetTotPoints(0);
6574
6575 Array<OneD, NekDouble> innarray(nelmtpnts, 0.0);
6576 Array<OneD, NekDouble> outarray(nelmtcoef, 0.0);
6577
6578 Array<OneD, NekDouble> MatQ_data;
6579 Array<OneD, NekDouble> MatC_data;
6580
6581 DNekMatSharedPtr tmpMatQ, tmpMatC;
6582
6583 nelmtcoef0 = nelmtcoef;
6584 nelmtpnts0 = nelmtpnts;
6585
6586 for (nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
6587 {
6588 nelmtcoef = GetNcoeffs(nelmt);
6589 nelmtpnts = GetTotPoints(nelmt);
6590
6591 tmpMatQ = ElmtJacQuad[nelmt];
6592 tmpMatC = ElmtJacCoef[nelmt];
6593
6594 MatQ_data = tmpMatQ->GetPtr();
6595 MatC_data = tmpMatC->GetPtr();
6596
6597 if (nelmtcoef != nelmtcoef0)
6598 {
6599 outarray = Array<OneD, NekDouble>(nelmtcoef, 0.0);
6600 nelmtcoef0 = nelmtcoef;
6601 }
6602
6603 if (nelmtpnts != nelmtpnts0)
6604 {
6605 innarray = Array<OneD, NekDouble>(nelmtpnts, 0.0);
6606 nelmtpnts0 = nelmtpnts;
6607 }
6608
6609 for (int np = 0; np < nelmtcoef; np++)
6610 {
6611 Vmath::Vcopy(nelmtpnts, &MatQ_data[0] + np, nelmtcoef, &innarray[0],
6612 1);
6613 (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray, innarray);
6614 (*m_exp)[nelmt]->IProductWRTBase(innarray, outarray);
6615
6616 Vmath::Vadd(nelmtcoef, &outarray[0], 1, &MatC_data[0] + np,
6617 nelmtcoef, &MatC_data[0] + np, nelmtcoef);
6618 }
6619 }
6620}
6621
6623 const NekDouble scale, const Array<OneD, NekDouble> &inarray,
6624 Array<OneD, NekDouble> &outarray)
6625{
6626 int cnt, cnt1;
6627
6628 cnt = cnt1 = 0;
6629
6630 switch (m_expType)
6631 {
6632 case e2D:
6633 {
6634 for (int i = 0; i < GetExpSize(); ++i)
6635 {
6636 // get new points key
6637 int pt0 = (*m_exp)[i]->GetNumPoints(0);
6638 int pt1 = (*m_exp)[i]->GetNumPoints(1);
6639 int npt0 = (int)pt0 * scale;
6640 int npt1 = (int)pt1 * scale;
6641
6642 LibUtilities::PointsKey newPointsKey0(
6643 npt0, (*m_exp)[i]->GetPointsType(0));
6644 LibUtilities::PointsKey newPointsKey1(
6645 npt1, (*m_exp)[i]->GetPointsType(1));
6646
6647 // Project points;
6649 newPointsKey0, newPointsKey1, &inarray[cnt],
6650 (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
6651 (*m_exp)[i]->GetBasis(1)->GetPointsKey(), &outarray[cnt1]);
6652
6653 cnt += npt0 * npt1;
6654 cnt1 += pt0 * pt1;
6655 }
6656 }
6657 break;
6658 case e3D:
6659 {
6660 for (int i = 0; i < GetExpSize(); ++i)
6661 {
6662 // get new points key
6663 int pt0 = (*m_exp)[i]->GetNumPoints(0);
6664 int pt1 = (*m_exp)[i]->GetNumPoints(1);
6665 int pt2 = (*m_exp)[i]->GetNumPoints(2);
6666 int npt0 = (int)pt0 * scale;
6667 int npt1 = (int)pt1 * scale;
6668 int npt2 = (int)pt2 * scale;
6669
6670 LibUtilities::PointsKey newPointsKey0(
6671 npt0, (*m_exp)[i]->GetPointsType(0));
6672 LibUtilities::PointsKey newPointsKey1(
6673 npt1, (*m_exp)[i]->GetPointsType(1));
6674 LibUtilities::PointsKey newPointsKey2(
6675 npt2, (*m_exp)[i]->GetPointsType(2));
6676
6677 // Project points;
6679 newPointsKey0, newPointsKey1, newPointsKey2, &inarray[cnt],
6680 (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
6681 (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
6682 (*m_exp)[i]->GetBasis(2)->GetPointsKey(), &outarray[cnt1]);
6683
6684 cnt += npt0 * npt1 * npt2;
6685 cnt1 += pt0 * pt1 * pt2;
6686 }
6687 }
6688 break;
6689 default:
6690 {
6691 NEKERROR(ErrorUtil::efatal, "not setup for this expansion");
6692 }
6693 break;
6694 }
6695}
6696
6698{
6699 NEKERROR(ErrorUtil::efatal, "v_GetLocTraceToTraceMap not coded");
6701}
6702
6703} // namespace Nektar::MultiRegions
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:243
#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 LIKWID_MARKER_START(regionTag)
Definition: Likwid.hpp:46
#define LIKWID_MARKER_STOP(regionTag)
Definition: Likwid.hpp:47
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
COLLECTIONS_EXPORT OperatorImpMap SetWithTimings(std::vector< StdRegions::StdExpansionSharedPtr > pGeom, OperatorImpMap &impTypes, bool verbose=true)
COLLECTIONS_EXPORT void UpdateOptFile(std::string sessName, LibUtilities::CommSharedPtr &comm)
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(StdRegions::StdExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Describes the specification for a Basis.
Definition: Basis.h:45
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:120
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:131
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:137
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:74
static const std::string GetFileType(const std::string &filename, CommSharedPtr comm)
Determine file type of given input file.
Definition: FieldIO.cpp:94
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Defines a specification for a set of points.
Definition: Points.h:50
size_t GetNumPoints() const
Definition: Points.h:85
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:70
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1083
void AddTraceJacToElmtJac(const Array< OneD, const DNekMatSharedPtr > &FwdMat, const Array< OneD, const DNekMatSharedPtr > &BwdMat, Array< OneD, DNekMatSharedPtr > &fieldMat)
inverse process of v_GetFwdBwdTracePhys. Given Trace integration of Fwd and Bwd Jacobian,...
Definition: ExpList.cpp:6272
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
const Array< OneD, const std::shared_ptr< ExpList > > & GetBndCondExpansions()
Definition: ExpList.h:2112
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
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Definition: ExpList.cpp:4368
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2656
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5090
LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > & GetGlobalLinSysManager(void)
Definition: ExpList.cpp:5995
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
Definition: ExpList.cpp:91
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
Definition: ExpList.cpp:5803
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.cpp:5697
void ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: ExpList.cpp:2391
std::shared_ptr< GlobalMatrix > GenGlobalMatrix(const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Generates a global matrix from the given key and map.
Definition: ExpList.cpp:2790
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5411
void AddRightIPTPhysDerivBase(const int dir, const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
Definition: ExpList.cpp:6508
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5176
int GetExpIndex(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0, bool returnNearestElmt=false, int cachedId=-1, NekDouble maxDistance=1e6)
This function returns the index of the local elemental expansion containing the arbitrary point given...
Definition: ExpList.cpp:3185
virtual 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)
Definition: ExpList.cpp:5199
virtual void v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:3515
BlockMatrixMapShPtr m_blockMat
Definition: ExpList.h:1128
virtual Array< OneD, const unsigned int > v_GetYIDs(void)
Definition: ExpList.cpp:3977
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
multiply the metric jacobi and quadrature weights
Definition: ExpList.cpp:1929
virtual void v_LinearAdvectionReactionSolve(const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:5226
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2261
void WriteVtkPieceFooter(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:3775
virtual void v_ExtractPhysToBnd(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
Definition: ExpList.cpp:5615
void GetBwdWeight(Array< OneD, NekDouble > &weightAver, Array< OneD, NekDouble > &weightJump)
Get the weight value for boundary conditions for boundary average and jump calculations.
Definition: ExpList.cpp:4405
virtual void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:6050
virtual std::vector< bool > & v_GetLeftAdjacentTraces(void)
Definition: ExpList.cpp:4646
void InitialiseExpVector(const SpatialDomains::ExpansionInfoMap &expmap)
Define a list of elements using the geometry and basis key information in expmap;.
Definition: ExpList.cpp:1676
void SetupCoeffPhys(bool DeclareCoeffPhysArrays=true, bool SetupOffsets=true)
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
Definition: ExpList.cpp:1619
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:5772
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2361
virtual void v_GetBCValues(Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
Definition: ExpList.cpp:5277
virtual void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Definition: ExpList.cpp:5143
virtual void v_WriteTecplotHeader(std::ostream &outfile, std::string var="")
Definition: ExpList.cpp:3384
virtual std::shared_ptr< ExpList > & v_GetPlane(int n)
Definition: ExpList.cpp:5792
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
virtual void v_AddTraceQuadPhysToOffDiag(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:5134
virtual void v_ClearGlobalLinSysManager(void)
Definition: ExpList.cpp:3985
virtual void v_UnsetGlobalLinSys(GlobalLinSysKey, bool)
Definition: ExpList.cpp:3997
void GeneralGetFieldDefinitions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef, int NumHomoDir=0, Array< OneD, LibUtilities::BasisSharedPtr > &HomoBasis=LibUtilities::NullBasisSharedPtr1DArray, std::vector< NekDouble > &HomoLen=LibUtilities::NullNekDoubleVector, bool homoStrips=false, std::vector< unsigned int > &HomoSIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoZIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoYIDs=LibUtilities::NullUnsignedIntVector)
Definition: ExpList.cpp:4085
virtual void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
Definition: ExpList.cpp:5764
virtual void v_CurlCurl(Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
Definition: ExpList.cpp:2293
virtual void v_LocalToGlobal(bool UseComm)
Definition: ExpList.cpp:5365
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions()
Definition: ExpList.cpp:5721
virtual LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > & v_GetGlobalLinSysManager(void)
Definition: ExpList.cpp:4005
virtual void v_DealiasedDotProd(const int num_dofs, const Array< OneD, Array< OneD, NekDouble > > &inarray1, const Array< OneD, Array< OneD, NekDouble > > &inarray2, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ExpList.cpp:5267
virtual ~ExpList()
The default destructor.
Definition: ExpList.cpp:1897
virtual void v_FwdTransLocalElmt(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2450
int GetPoolCount(std::string)
Definition: ExpList.cpp:5984
virtual void v_ImposeDirichletConditions(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5339
Array< OneD, const unsigned int > GetZIDs(void)
This function returns a vector containing the wave numbers in z-direction associated with the 3D homo...
Definition: ExpList.h:590
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
Definition: ExpList.cpp:5742
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
Definition: ExpList.cpp:5168
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1123
std::shared_ptr< AssemblyMapDG > & GetTraceMap(void)
Definition: ExpList.h:2157
virtual void v_GetMovingFrames(const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ExpList.cpp:4418
size_t GetNumElmts(void)
This function returns the number of elements in the expansion which may be different for a homogeoeno...
Definition: ExpList.h:653
virtual std::shared_ptr< AssemblyMapDG > & v_GetTraceMap()
Definition: ExpList.cpp:4625
Array< OneD, std::pair< int, int > > m_coeffsToElmt
m_coeffs to elemental value map
Definition: ExpList.h:1127
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:5099
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
Definition: ExpList.cpp:3304
virtual GlobalLinSysKey v_LinearAdvectionDiffusionReactionSolve(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)
Definition: ExpList.cpp:5212
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Definition: ExpList.cpp:3905
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs, std::unordered_map< int, int > zIdToPlane)
Extract data from raw field data into expansion list.
Definition: ExpList.cpp:4280
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo(void)
Definition: ExpList.cpp:5754
virtual const Array< OneD, const int > & v_GetTraceBndMap()
Definition: ExpList.cpp:4641
virtual void v_GetBoundaryNormals(int i, Array< OneD, Array< OneD, NekDouble > > &normals)
Definition: ExpList.cpp:5652
virtual Array< OneD, const NekDouble > v_HomogeneousEnergy(void)
Definition: ExpList.cpp:3939
virtual NekDouble v_GetHomoLen(void)
Definition: ExpList.cpp:3955
const LocTraceToTraceMapSharedPtr & GetLocTraceToTraceMap() const
Definition: ExpList.h:2123
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
Definition: ExpList.cpp:3969
void ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Extract the data from fromField using fromExpList the coeffs using the basic ExpList Elemental expans...
Definition: ExpList.cpp:4264
void GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
This function calculates the coordinates of all the elemental quadrature points .
Definition: ExpList.h:1770
void ResetMatrices()
Reset matrices.
Definition: ExpList.cpp:3367
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
virtual void v_SetBndCondBwdWeight(const int index, const NekDouble value)
Definition: ExpList.cpp:5161
void UnsetGlobalLinSys(GlobalLinSysKey, bool)
Definition: ExpList.cpp:5989
virtual void v_SmoothField(Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:2482
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
Definition: ExpList.cpp:4212
std::vector< bool > m_collectionsDoInit
Vector of bools to act as an initialise on first call flag.
Definition: ExpList.h:1121
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1107
virtual void v_HomogeneousFwdTrans(const int npts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:5237
void ExtractCoeffsFromFile(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, Array< OneD, NekDouble > &coeffs)
Definition: ExpList.cpp:4012
virtual void v_DealiasedProd(const int num_dofs, const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5257
void AddRightIPTBaseMatrix(const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
Definition: ExpList.cpp:6565
void GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Definition: ExpList.h:2209
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
Definition: ExpList.cpp:3782
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
Definition: ExpList.cpp:3642
std::shared_ptr< ExpList > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: ExpList.h:956
int GetShapeDimension()
This function returns the dimension of the shape of the element eid.
Definition: ExpList.h:1912
NekDouble H1(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Calculates the error of the global spectral/hp element approximation.
Definition: ExpList.cpp:4067
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1056
std::shared_ptr< GlobalLinSys > GenGlobalLinSys(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
This operation constructs the global linear system of type mkey.
Definition: ExpList.cpp:3075
void WriteVtkHeader(std::ostream &outfile)
Definition: ExpList.cpp:3628
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 GenerateElementVector(const int ElementID, const NekDouble scalar1, const NekDouble scalar2, Array< OneD, NekDouble > &outarray)
Generate vector v such that v[i] = scalar1 if i is in the element < ElementID. Otherwise,...
Definition: ExpList.cpp:4462
std::shared_ptr< DNekMat > GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Definition: ExpList.cpp:2923
virtual void v_NormVectorIProductWRTBase(Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
Definition: ExpList.cpp:5286
virtual void v_ExtractElmtToBndPhys(const int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
Definition: ExpList.cpp:5538
virtual void v_ExtractPhysToBndElmt(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
Definition: ExpList.cpp:5576
virtual std::shared_ptr< InterfaceMapDG > & v_GetInterfaceMap()
Definition: ExpList.cpp:4633
std::map< int, RobinBCInfoSharedPtr > GetRobinBCInfo()
Definition: ExpList.h:885
virtual void v_Upwind(const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.cpp:4519
virtual void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Definition: ExpList.cpp:5704
virtual void v_Curl(Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
Definition: ExpList.cpp:2244
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:6622
void GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the result of the multiplication of a matrix of type specified by mkey with ...
Definition: ExpList.cpp:2700
virtual void v_SetUpPhysNormals()
: Set up a normal along the trace elements between two elements at elemental level
Definition: ExpList.cpp:5515
virtual void v_GlobalToLocal(void)
Definition: ExpList.cpp:5380
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
Definition: ExpList.cpp:3947
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
virtual void v_MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5191
virtual void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:5125
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight()
Definition: ExpList.cpp:5153
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Definition: ExpList.cpp:3860
void WriteVtkFooter(std::ostream &outfile)
Definition: ExpList.cpp:3636
void MultiplyByBlockMatrix(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1909
virtual std::shared_ptr< ExpList > & v_GetTrace()
Definition: ExpList.cpp:4617
void GetElmtNormalLength(Array< OneD, NekDouble > &lengthsFwd, Array< OneD, NekDouble > &lengthsBwd)
Get the length of elements in boundary normal direction.
Definition: ExpList.cpp:4935
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Definition: ExpList.cpp:3919
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
Definition: ExpList.cpp:4197
void IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Directional derivative along a given direction.
Definition: ExpList.cpp:1992
void FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Fill Bwd with boundary conditions.
Definition: ExpList.h:2216
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:3129
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions()
Definition: ExpList.cpp:5732
void GetMatIpwrtDeriveBase(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const int nDirctn, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
Definition: ExpList.cpp:6058
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
Definition: ExpList.cpp:5455
virtual void v_WriteTecplotField(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:3595
virtual void v_FillBndCondFromField(const Array< OneD, NekDouble > coeffs)
Definition: ExpList.cpp:5348
Collections::CollectionVector m_collections
Definition: ExpList.h:1119
virtual int v_GetPoolCount(std::string)
Definition: ExpList.cpp:3991
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:6000
virtual void v_HomogeneousBwdTrans(const int npts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:5247
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:3327
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: ExpList.cpp:2146
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2459
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
Definition: ExpList.h:1135
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:3422
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1058
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:5711
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
Definition: ExpList.cpp:5116
void DivideByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Divided by the metric jacobi and quadrature weights.
Definition: ExpList.cpp:1946
virtual const Array< OneD, const std::shared_ptr< ExpList > > & v_GetBndCondExpansions(void)
Definition: ExpList.cpp:4492
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1125
virtual void v_SetHomoLen(const NekDouble lhom)
Definition: ExpList.cpp:3963
void Upwind(const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.h:2146
virtual void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
Populate normals with the normals of all expansions.
Definition: ExpList.cpp:4726
const DNekScalBlkMatSharedPtr GenBlockMatrix(const GlobalMatrixKey &gkey)
This function assembles the block diagonal matrix of local matrices of the type mtype.
Definition: ExpList.cpp:2521
void PhysDeriv(Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
Definition: ExpList.h:1813
void MultiplyByElmtInvMass(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally mulplies the coefficient space of Sin my the elemental inverse of the mass ...
Definition: ExpList.cpp:2412
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
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Definition: ExpList.cpp:5528
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5394
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion(int i)
Definition: ExpList.cpp:4500
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
virtual const std::shared_ptr< LocTraceToTraceMap > & v_GetLocTraceToTraceMap(void) const
Definition: ExpList.cpp:6697
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1973
void ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs, std::unordered_map< int, int > zIdToPlane=std::unordered_map< int, int >())
Extract the data in fielddata into the coeffs.
Definition: ExpList.cpp:4256
ExpansionType m_expType
Expansion type.
Definition: ExpList.h:1051
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5082
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1099
ExpansionType GetExpType(void)
Returns the type of the expansion.
Definition: ExpList.cpp:1892
NekDouble Linf(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global spectral/hp element approximation.
Definition: ExpList.cpp:3822
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1557
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:1903
void GetDiagMatIpwrtBase(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
Definition: ExpList.cpp:6201
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Describes a matrix with ordering defined by a local to global map.
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
LibUtilities::ShapeType GetShapeType() const
Return the expansion type associated with key.
const StdRegions::VarCoeffMap & GetVarCoeffs() const
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Class representing a prismatic element in reference space.
Definition: StdPrismExp.h:45
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:149
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition: Operator.h:131
@ eLinearAdvectionDiffusionReaction
Definition: Operator.h:66
void PhysGalerkinProject3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
std::vector< BasisKey > BasisKeyVector
Name for a vector of BasisKeys.
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:47
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:184
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:101
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset=0)
Definition: ShapeType.hpp:305
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
Definition: FieldIO.cpp:69
std::shared_ptr< Transposition > TranspositionSharedPtr
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:81
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
void PhysGalerkinProject2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
std::shared_ptr< PointExp > PointExpSharedPtr
Definition: PointExp.h:132
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::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:68
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
void AlignFace(const StdRegions::Orientation orient, const int nquad1, const int nquad2, const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out)
Helper function to re-align face to a given orientation.
Definition: ExpList.cpp:4657
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1499
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:50
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
const char *const GlobalSysSolnTypeMap[]
static LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > NullGlobalLinSysManager
Definition: ExpList.h:1503
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51
std::shared_ptr< GlobalMatrix > GlobalMatrixSharedPtr
Shared pointer to a GlobalMatrix object.
Definition: GlobalMatrix.h:84
GlobalLinSysFactory & GetGlobalLinSysFactory()
static LocTraceToTraceMapSharedPtr NullLocTraceToTraceMapSharedPtr
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::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:50
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
std::map< GlobalMatrixKey, DNekScalBlkMatSharedPtr > BlockMatrixMap
A map between global matrix keys and their associated block matrices.
Definition: ExpList.h:91
static const NekDouble kNekZeroTol
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:45
GeomMMF
Principle direction for MMF.
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:213
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:82
std::shared_ptr< std::vector< std::pair< GeometrySharedPtr, int > > > GeometryLinkSharedPtr
Definition: MeshGraph.h:166
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:220
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:58
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:75
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition: MeshGraph.h:140
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:219
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:51
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:56
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:61
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:141
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:136
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
VarCoeffMap RestrictCoeffMap(const VarCoeffMap &m, size_t offset, size_t cnt)
Definition: StdRegions.hpp:378
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
std::vector< double > z(NPUPPER)
std::vector< double > d(NPUPPER *NPUPPER)
StdRegions::ConstFactorMap factors
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:68
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::pair< IndexType, IndexType > CoordType
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC
std::map< CoordType, NekDouble > COOMatType
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition: Vmath.hpp:396
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 Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.hpp:539
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
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
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 > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:69
std::shared_ptr< RobinBCInfo > next