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