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