Nektar++
LocTraceToTraceMap.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: LocTraceToTraceMap.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: Local trace to general trace mapping information
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <LibUtilities/Foundations/ManagerAccess.h>
36#include <LibUtilities/LibUtilitiesDeclspec.h>
37#include <LibUtilities/LinearAlgebra/Blas.hpp>
38#include <LocalRegions/Expansion1D.h>
39#include <LocalRegions/Expansion2D.h>
40#include <LocalRegions/Expansion3D.h>
41#include <MultiRegions/AssemblyMap/AssemblyMap.h>
42#include <MultiRegions/AssemblyMap/LocTraceToTraceMap.h>
43#include <MultiRegions/ExpList.h>
44
45using namespace std;
46
47namespace Nektar::MultiRegions
48{
49
50/**
51 * @brief Set up trace to trace mapping components.
52 *
53 * @param locExp Expansion list of full dimension problem.
54 * @param trace Expansion list of one dimension lower trace.
55 * @param elmtToTrace Mapping from elemental facets to trace.
56 * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
57 *
58 * @todo Add 1D support
59 */
60LocTraceToTraceMap::LocTraceToTraceMap(
61 const ExpList &locExp, const ExpListSharedPtr &trace,
62 const Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
63 &elmtToTrace,
64 const vector<bool> &LeftAdjacents)
65{
66 const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
67
68 // Assume that all the elements have same dimension
69 m_expdim = locExpVector[0]->GetShapeDimension();
70
71 // set up interpolation details for all dimension elements.
72 Setup(locExp, trace, elmtToTrace, LeftAdjacents);
73}
74
75LocTraceToTraceMap::~LocTraceToTraceMap()
76{
77}
78
79/**
80 * @brief Set up member variables for a two-dimensional problem.
81 *
82 * @param locExp Expansion list of elements
83 * @param trace Expansion list of the trace.
84 * @param elmtToTrace Mapping from elemental trace to unique trace.
85 * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
86 */
87void LocTraceToTraceMap::Setup(
88 const ExpList &locExp, const ExpListSharedPtr &trace,
89 const Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
90 &elmtToTrace,
91 const vector<bool> &LeftAdjacents)
92{
93 m_locInterpTraceToTraceMap = Array<OneD, Array<OneD, int>>(2);
94 m_locTraceToElmtTraceMap = Array<OneD, Array<OneD, int>>(2);
95 m_interpTrace = Array<OneD, Array<OneD, InterpLocTraceToTrace>>(2);
96 m_interpTraceI0 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
97 m_interpEndPtI0 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
98 m_interpFromTraceI0 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
99 m_interpPoints = Array<OneD, Array<OneD, TraceInterpPoints>>(2);
100 m_interpNtraces = Array<OneD, Array<OneD, int>>(2);
101
102 if (m_expdim == 3)
103 {
104 m_interpTraceI1 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
105 m_interpFromTraceI1 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
106 m_interpEndPtI1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
107 }
108
109 m_traceCoeffsToElmtMap = Array<OneD, Array<OneD, int>>(2);
110 m_traceCoeffsToElmtTrace = Array<OneD, Array<OneD, int>>(2);
111 m_traceCoeffsToElmtSign = Array<OneD, Array<OneD, int>>(2);
112
113 LocalRegions::ExpansionSharedPtr elmt;
114 const std::shared_ptr<LocalRegions::ExpansionVector> exp = locExp.GetExp();
115
116 int cnt, n, e, phys_offset;
117
118 int nexp = exp->size();
119 m_nTracePts = trace->GetTotPoints();
120
121 // Count number of traces and points required for maps
122 int nFwdPts = 0;
123 int nBwdPts = 0;
124 int nFwdCoeffs = 0;
125 int nBwdCoeffs = 0;
126 m_nFwdLocTracePts = 0;
127 m_nLocTracePts = 0;
128
129 for (cnt = n = 0; n < nexp; ++n)
130 {
131 elmt = (*exp)[n];
132
133 for (int i = 0; i < elmt->GetNtraces(); ++i, ++cnt)
134 {
135 int nLocPts = elmt->GetTraceNumPoints(i);
136 m_nLocTracePts += nLocPts;
137
138 if (LeftAdjacents[cnt])
139 {
140 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
141 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
142 m_nFwdLocTracePts += nLocPts;
143 }
144 else
145 {
146 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
147 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
148 }
149 }
150 }
151
152 m_locTraceToFieldMap = Array<OneD, int>(m_nLocTracePts);
153
154 m_locTraceToElmtTraceMap[0] = Array<OneD, int>(m_nFwdLocTracePts);
155 m_locTraceToElmtTraceMap[1] =
156 Array<OneD, int>(m_nLocTracePts - m_nFwdLocTracePts);
157
158 m_locInterpTraceToTraceMap[0] = Array<OneD, int>(nFwdPts);
159 m_locInterpTraceToTraceMap[1] = Array<OneD, int>(nBwdPts);
160
161 m_nTraceCoeffs[0] = nFwdCoeffs;
162 m_nTraceCoeffs[1] = nBwdCoeffs;
163
164 m_traceCoeffsToElmtMap[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
165 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
166 m_traceCoeffsToElmtTrace[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
167 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
168 m_traceCoeffsToElmtSign[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
169 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
170
171 // Gather information about trace interpolations
172 map<TraceInterpPoints, vector<pair<int, int>>, cmpop> TraceInterpMap;
173
174 vector<vector<int>> TraceOrder;
175 TraceOrder.resize(nexp);
176 vector<vector<int>> ElmtPhysTraceOffset;
177 ElmtPhysTraceOffset.resize(nexp);
178 int ntrace;
179 int fwdcnt = 0;
180 int bwdcnt = 0;
181 int neoffset = 0;
182 // Generate a map of similar traces with the same
183 // interpolation requirements
184
185 for (cnt = n = 0; n < nexp; ++n)
186 {
187 elmt = (*exp)[n];
188 ntrace = elmt->GetNtraces();
189 TraceOrder[n].resize(ntrace);
190 ElmtPhysTraceOffset[n].resize(ntrace);
191
192 int coeffoffset = locExp.GetCoeff_Offset(n);
193 for (e = 0; e < ntrace; ++e, ++cnt)
194 {
195 LocalRegions::ExpansionSharedPtr elmttrace = elmtToTrace[n][e];
196 StdRegions::Orientation orient = elmt->GetTraceOrient(e);
197
198 LibUtilities::PointsKey fromPointsKey0, fromPointsKey1;
199 LibUtilities::PointsKey toPointsKey0, toPointsKey1;
200 Array<OneD, int> P(2, -1);
201
202 switch (m_expdim)
203 {
204 case 1:
205 {
206 fromPointsKey0 = elmt->GetBasis(0)->GetPointsKey();
207 fromPointsKey1 =
208 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
209 // dummy info since no interpolation is required in this
210 // case.
211 toPointsKey0 =
212 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
213 toPointsKey1 =
214 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
215 }
216 break;
217 case 2:
218 {
219 int dir0 = elmt->GetGeom()->GetDir(e, 0);
220
221 fromPointsKey0 = elmt->GetBasis(dir0)->GetPointsKey();
222 fromPointsKey1 =
223 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
224
225 toPointsKey0 = elmttrace->GetBasis(0)->GetPointsKey();
226 toPointsKey1 =
227 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
228
229 P[0] = elmttrace->GetBasisNumModes(0);
230 }
231 break;
232 case 3:
233 {
234 int dir0 = elmt->GetGeom()->GetDir(e, 0);
235 int dir1 = elmt->GetGeom()->GetDir(e, 1);
236
237 fromPointsKey0 = elmt->GetBasis(dir0)->GetPointsKey();
238 fromPointsKey1 = elmt->GetBasis(dir1)->GetPointsKey();
239
240 if (orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
241 {
242 toPointsKey0 = elmttrace->GetBasis(0)->GetPointsKey();
243 toPointsKey1 = elmttrace->GetBasis(1)->GetPointsKey();
244 P[0] = elmttrace->GetBasisNumModes(0);
245 P[1] = elmttrace->GetBasisNumModes(1);
246 }
247 else // transpose points key evaluation
248 {
249 toPointsKey0 = elmttrace->GetBasis(1)->GetPointsKey();
250 toPointsKey1 = elmttrace->GetBasis(0)->GetPointsKey();
251 P[0] = elmttrace->GetBasisNumModes(1);
252 P[1] = elmttrace->GetBasisNumModes(0);
253 }
254 }
255 break;
256 }
257
258 TraceInterpPoints fpoint(fromPointsKey0, fromPointsKey1,
259 toPointsKey0, toPointsKey1);
260
261 pair<int, int> epf(n, e);
262 TraceInterpMap[fpoint].push_back(epf);
263 TraceOrder[n][e] = cnt;
264
265 ElmtPhysTraceOffset[n][e] = neoffset;
266 neoffset += elmt->GetTraceNumPoints(e);
267
268 // Setup for coefficient mapping from trace normal flux
269 // to elements
270 Array<OneD, unsigned int> map;
271 Array<OneD, int> sign;
272 // Test shows we should swap P0 and P1 before calling
273 // GetTraceToElementMap if orientation is transposed.
274 // This is contrary to ReOrientTracePhysMap
275 elmt->GetTraceToElementMap(e, map, sign, orient, P[0], P[1]);
276
277 int order_t = elmttrace->GetNcoeffs();
278 int t_offset = trace->GetCoeff_Offset(elmttrace->GetElmtId());
279
280 double fac = 1.0;
281
282 if (elmt->GetTraceExp(e)->GetRightAdjacentElementExp())
283 {
284 if (elmttrace->GetRightAdjacentElementExp()
285 ->GetGeom()
286 ->GetGlobalID() == elmt->GetGeom()->GetGlobalID())
287 {
288 fac = -1.0;
289 }
290 }
291
292 if (LeftAdjacents[cnt])
293 {
294 for (int i = 0; i < order_t; ++i)
295 {
296 m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
297 m_traceCoeffsToElmtTrace[0][fwdcnt] = t_offset + i;
298 m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
299 }
300 }
301 else
302 {
303 for (int i = 0; i < order_t; ++i)
304 {
305 m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
306 m_traceCoeffsToElmtTrace[1][bwdcnt] = t_offset + i;
307 m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
308 }
309 }
310 }
311 }
312
313 int nInterpType = TraceInterpMap.size();
314
315 // need to decide on 1D case here !!!!!
316 for (int i = 0; i < 2; ++i)
317 {
318 m_interpTrace[i] = Array<OneD, InterpLocTraceToTrace>(nInterpType);
319 m_interpTraceI0[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
320 m_interpEndPtI0[i] = Array<OneD, Array<OneD, NekDouble>>(nInterpType);
321 m_interpFromTraceI0[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
322 m_interpPoints[i] = Array<OneD, TraceInterpPoints>(nInterpType);
323 m_interpNtraces[i] = Array<OneD, int>(nInterpType, 0);
324 }
325
326 if (m_expdim > 2)
327 {
328 for (int i = 0; i < 2; ++i)
329 {
330 m_interpTraceI1[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
331 m_interpFromTraceI1[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
332 m_interpEndPtI1[i] =
333 Array<OneD, Array<OneD, NekDouble>>(nInterpType);
334 }
335 }
336
337 int ntracepts, ntracepts1;
338 int cnt1 = 0;
339 int cnt2 = 0;
340 int cntFwd = 0;
341 int cntBwd = 0;
342 int cntFwd1 = 0;
343 int cntBwd1 = 0;
344 int set;
345 Array<OneD, int> traceids;
346 Array<OneD, int> locTraceToTraceMap;
347 cnt = 0;
348
349 for (auto it = TraceInterpMap.begin(); it != TraceInterpMap.end();
350 ++it, ++cnt1)
351 {
352 LibUtilities::PointsKey fromPointsKey0 = std::get<0>(it->first);
353 LibUtilities::PointsKey fromPointsKey1 = std::get<1>(it->first);
354 LibUtilities::PointsKey toPointsKey0 = std::get<2>(it->first);
355 LibUtilities::PointsKey toPointsKey1 = std::get<3>(it->first);
356
357 bool fwdSet = false;
358 bool bwdSet = false;
359
360 for (int f = 0; f < it->second.size(); ++f, ++cnt2)
361 {
362 n = it->second[f].first;
363 e = it->second[f].second;
364
365 StdRegions::StdExpansionSharedPtr elmttrace = elmtToTrace[n][e];
366
367 elmt = (*exp)[n];
368 phys_offset = locExp.GetPhys_Offset(n);
369
370 // Mapping of new edge order to one that loops over elmts
371 // then set up mapping of faces in standard cartesian order
372 elmt->GetTracePhysMap(e, traceids);
373
374 ntracepts = elmt->GetTraceNumPoints(e);
375 ntracepts1 = elmttrace->GetTotPoints();
376
377 StdRegions::Orientation orient = elmt->GetTraceOrient(e);
378
379 // toPoints have already been swapped. But here we need original
380 // elmttrace points (w.r.t local axes). So swap back if orient >= 9
381 if (orient >= 9)
382 {
383 elmt->ReOrientTracePhysMap(orient, locTraceToTraceMap,
384 toPointsKey1.GetNumPoints(),
385 toPointsKey0.GetNumPoints());
386 }
387 else
388 {
389 elmt->ReOrientTracePhysMap(orient, locTraceToTraceMap,
390 toPointsKey0.GetNumPoints(),
391 toPointsKey1.GetNumPoints());
392 }
393
394 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
395
396 if (LeftAdjacents[TraceOrder[n][e]])
397 {
398 for (int i = 0; i < ntracepts; ++i)
399 {
400 m_locTraceToFieldMap[cntFwd + i] =
401 phys_offset + traceids[i];
402 }
403
404 for (int i = 0; i < ntracepts; ++i)
405 {
406 m_locTraceToElmtTraceMap[0][cntFwd + i] =
407 ElmtPhysTraceOffset[n][e] + i;
408 }
409
410 for (int i = 0; i < ntracepts1; ++i)
411 {
412 m_locInterpTraceToTraceMap[0][cntFwd1 + i] =
413 offset + locTraceToTraceMap[i];
414 }
415
416 cntFwd += ntracepts;
417 cntFwd1 += ntracepts1;
418 set = 0;
419 }
420 else
421 {
422 for (int i = 0; i < ntracepts; ++i)
423 {
424 m_locTraceToFieldMap[m_nFwdLocTracePts + cntBwd + i] =
425 phys_offset + traceids[i];
426 }
427
428 for (int i = 0; i < ntracepts; ++i)
429 {
430 m_locTraceToElmtTraceMap[1][cntBwd + i] =
431 ElmtPhysTraceOffset[n][e] + i;
432 }
433
434 for (int i = 0; i < ntracepts1; ++i)
435 {
436 m_locInterpTraceToTraceMap[1][cntBwd1 + i] =
437 offset + locTraceToTraceMap[i];
438 }
439
440 cntBwd += ntracepts;
441 cntBwd1 += ntracepts1;
442 set = 1;
443 }
444
445 m_interpNtraces[set][cnt1] += 1;
446
447 if ((fwdSet == false && set == 0) || (bwdSet == false && set == 1))
448 {
449 m_interpPoints[set][cnt1] = it->first;
450
451 switch (m_expdim)
452 {
453 case 1:
454 {
455 // Always no interplation in this case
456 m_interpTrace[set][cnt1] = eNoInterp;
457 }
458 break;
459 case 2:
460 {
461 if (fromPointsKey0 == toPointsKey0)
462 {
463 m_interpTrace[set][cnt1] = eNoInterp;
464 }
465 else
466 {
467 m_interpTrace[set][cnt1] = eInterpDir0;
468 m_interpTraceI0[set][cnt1] =
469 LibUtilities::PointsManager()[fromPointsKey0]
470 ->GetI(toPointsKey0);
471 m_interpFromTraceI0[set][cnt1] =
472 LibUtilities::PointsManager()[toPointsKey0]
473 ->GetI(fromPointsKey0);
474 // Check to see if we can
475 // just interpolate endpoint
476 if ((fromPointsKey0.GetPointsType() ==
477 LibUtilities::eGaussRadauMAlpha1Beta0) &&
478 (toPointsKey0.GetPointsType() ==
479 LibUtilities::eGaussLobattoLegendre))
480 {
481 if (fromPointsKey0.GetNumPoints() + 1 ==
482 toPointsKey0.GetNumPoints())
483 {
484 m_interpTrace[set][cnt1] = eInterpEndPtDir0;
485
486 int fnp0 = fromPointsKey0.GetNumPoints();
487 int tnp0 = toPointsKey0.GetNumPoints();
488
489 m_interpEndPtI0[set][cnt1] =
490 Array<OneD, NekDouble>(fnp0);
491
492 Vmath::Vcopy(fnp0,
493 m_interpTraceI0[set][cnt1]
494 ->GetPtr()
495 .data() +
496 tnp0 - 1,
497 tnp0,
498 &m_interpEndPtI0[set][cnt1][0],
499 1);
500 }
501 }
502 }
503 }
504 break;
505 case 3:
506 {
507 if (fromPointsKey0 == toPointsKey0)
508 {
509 if (fromPointsKey1 == toPointsKey1)
510 {
511 m_interpTrace[set][cnt1] = eNoInterp;
512 }
513 else
514 {
515 m_interpTrace[set][cnt1] = eInterpDir1;
516 m_interpTraceI1[set][cnt1] =
517 LibUtilities::PointsManager()
518 [fromPointsKey1]
519 ->GetI(toPointsKey1);
520 m_interpFromTraceI1[set][cnt1] =
521 LibUtilities::PointsManager()[toPointsKey1]
522 ->GetI(fromPointsKey1);
523
524 // Check to see if we can just
525 // interpolate endpoint
526 if ((fromPointsKey1.GetPointsType() ==
527 LibUtilities::eGaussRadauMAlpha1Beta0) &&
528 (toPointsKey1.GetPointsType() ==
529 LibUtilities::eGaussLobattoLegendre))
530 {
531 if (fromPointsKey1.GetNumPoints() + 1 ==
532 toPointsKey1.GetNumPoints())
533 {
534 m_interpTrace[set][cnt1] =
535 eInterpEndPtDir1;
536 int fnp1 =
537 fromPointsKey1.GetNumPoints();
538 int tnp1 = toPointsKey1.GetNumPoints();
539 m_interpEndPtI1[set][cnt1] =
540 Array<OneD, NekDouble>(fnp1);
541 Vmath::Vcopy(
542 fnp1,
543 m_interpTraceI1[set][cnt1]
544 ->GetPtr()
545 .data() +
546 tnp1 - 1,
547 tnp1,
548 &m_interpEndPtI1[set][cnt1][0], 1);
549 }
550 }
551 }
552 }
553 else
554 {
555 if (fromPointsKey1 == toPointsKey1)
556 {
557 m_interpTrace[set][cnt1] = eInterpDir0;
558 m_interpTraceI0[set][cnt1] =
559 LibUtilities::PointsManager()
560 [fromPointsKey0]
561 ->GetI(toPointsKey0);
562 m_interpFromTraceI0[set][cnt1] =
563 LibUtilities::PointsManager()[toPointsKey0]
564 ->GetI(fromPointsKey0);
565
566 // Check to see if we can just
567 // interpolate endpoint
568 if ((fromPointsKey0.GetPointsType() ==
569 LibUtilities::eGaussRadauMAlpha1Beta0) &&
570 (toPointsKey0.GetPointsType() ==
571 LibUtilities::eGaussLobattoLegendre))
572 {
573 if (fromPointsKey0.GetNumPoints() + 1 ==
574 toPointsKey0.GetNumPoints())
575 {
576 m_interpTrace[set][cnt1] =
577 eInterpEndPtDir0;
578 int fnp0 =
579 fromPointsKey0.GetNumPoints();
580 int tnp0 = toPointsKey0.GetNumPoints();
581 m_interpEndPtI0[set][cnt1] =
582 Array<OneD, NekDouble>(fnp0);
583 Vmath::Vcopy(
584 fnp0,
585 m_interpTraceI0[set][cnt1]
586 ->GetPtr()
587 .data() +
588 tnp0 - 1,
589 tnp0,
590 &m_interpEndPtI0[set][cnt1][0], 1);
591 }
592 }
593 }
594 else
595 {
596 m_interpTrace[set][cnt1] = eInterpBothDirs;
597 m_interpTraceI0[set][cnt1] =
598 LibUtilities::PointsManager()
599 [fromPointsKey0]
600 ->GetI(toPointsKey0);
601 m_interpFromTraceI0[set][cnt1] =
602 LibUtilities::PointsManager()[toPointsKey0]
603 ->GetI(fromPointsKey0);
604 m_interpTraceI1[set][cnt1] =
605 LibUtilities::PointsManager()
606 [fromPointsKey1]
607 ->GetI(toPointsKey1);
608 m_interpFromTraceI1[set][cnt1] =
609 LibUtilities::PointsManager()[toPointsKey1]
610 ->GetI(fromPointsKey1);
611
612 // check to see if we can just
613 // interpolate endpoint
614 if ((fromPointsKey0.GetPointsType() ==
615 LibUtilities::eGaussRadauMAlpha1Beta0) &&
616 (toPointsKey0.GetPointsType() ==
617 LibUtilities::eGaussLobattoLegendre))
618 {
619 if (fromPointsKey0.GetNumPoints() + 1 ==
620 toPointsKey0.GetNumPoints())
621 {
622 m_interpTrace[set][cnt1] =
623 eInterpEndPtDir0InterpDir1;
624 int fnp0 =
625 fromPointsKey0.GetNumPoints();
626 int tnp0 = toPointsKey0.GetNumPoints();
627 m_interpEndPtI0[set][cnt1] =
628 Array<OneD, NekDouble>(fnp0);
629 Vmath::Vcopy(
630 fnp0,
631 m_interpTraceI0[set][cnt1]
632 ->GetPtr()
633 .data() +
634 tnp0 - 1,
635 tnp0,
636 &m_interpEndPtI0[set][cnt1][0], 1);
637 }
638 }
639 }
640 }
641 }
642 }
643
644 if (set == 0)
645 {
646 fwdSet = true;
647 }
648 else
649 {
650 bwdSet = true;
651 }
652 }
653 }
654 }
655
656 TraceLocToElmtLocCoeffMap(locExp, trace);
657 FindElmtNeighbors(locExp, trace);
658}
659
660void LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap(
661 const ExpListSharedPtr &tracelist, const int ndim)
662{
663 switch (ndim)
664 {
665 case 2:
666 CalcLocTracePhysToTraceIDMap_2D(tracelist);
667 break;
668 case 3:
669 CalcLocTracePhysToTraceIDMap_3D(tracelist);
670 break;
671 default:
672 NEKERROR(ErrorUtil::efatal,
673 "CalcLocTracePhysToTraceIDMap not coded");
674 }
675}
676
677void LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_2D(
678 const ExpListSharedPtr &tracelist)
679{
680 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
681 tracelist->GetExp();
682 int ntotTrace = (*traceExp).size();
683 int ntPnts, noffset;
684
685 m_LocTracephysToTraceIDMap = Array<OneD, Array<OneD, int>>(2);
686 m_LocTracephysToTraceIDMap[0] = Array<OneD, int>(m_nFwdLocTracePts, -1);
687 m_LocTracephysToTraceIDMap[1] =
688 Array<OneD, int>(m_nLocTracePts - m_nFwdLocTracePts, -1);
689
690 Array<OneD, NekDouble> tracePnts(m_nTracePts, 0.0);
691 for (int nt = 0; nt < ntotTrace; nt++)
692 {
693 ntPnts = tracelist->GetTotPoints(nt);
694 noffset = tracelist->GetPhys_Offset(nt);
695 for (int i = 0; i < ntPnts; i++)
696 {
697 tracePnts[noffset + i] = NekDouble(nt);
698 }
699 }
700
701 Array<OneD, Array<OneD, NekDouble>> loctracePntsLR(2);
702 loctracePntsLR[0] = Array<OneD, NekDouble>(m_nFwdLocTracePts, 0.0);
703 loctracePntsLR[1] =
704 Array<OneD, NekDouble>(m_nLocTracePts - m_nFwdLocTracePts, 0.0);
705
706 for (int dir = 0; dir < 2; dir++)
707 {
708 int cnt = 0;
709 int cnt1 = 0;
710
711 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
712 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(),
713 tracePnts.data(), m_locInterpTraceToTraceMap[dir].data(),
714 tmp.data());
715
716 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
717 {
718 if (m_interpNtraces[dir][i])
719 {
720 LibUtilities::PointsKey fromPointsKey0 =
721 std::get<0>(m_interpPoints[dir][i]);
722 LibUtilities::PointsKey toPointsKey0 =
723 std::get<2>(m_interpPoints[dir][i]);
724
725 int fnp = fromPointsKey0.GetNumPoints();
726 int tnp = toPointsKey0.GetNumPoints();
727 int nedges = m_interpNtraces[dir][i];
728
729 for (int ne = 0; ne < nedges; ne++)
730 {
731 Vmath::Fill(fnp, tmp[cnt1], &loctracePntsLR[dir][cnt], 1);
732 cnt += fnp;
733 cnt1 += tnp;
734 }
735 }
736 }
737 }
738
739 NekDouble error = 0.0;
740 for (int nlr = 0; nlr < 2; nlr++)
741 {
742 for (int i = 0; i < loctracePntsLR[nlr].size(); i++)
743 {
744 m_LocTracephysToTraceIDMap[nlr][i] =
745 std::round(loctracePntsLR[nlr][i]);
746 error += abs(loctracePntsLR[nlr][i] -
747 NekDouble(m_LocTracephysToTraceIDMap[nlr][i]));
748 }
749 }
750 error = error / NekDouble(m_nLocTracePts);
751 ASSERTL0(error < NekConstants::kNekZeroTol,
752 "m_LocTracephysToTraceIDMap may not be integer !!");
753}
754
755void LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_3D(
756 const ExpListSharedPtr &tracelist)
757{
758 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
759 tracelist->GetExp();
760 int ntotTrace = (*traceExp).size();
761 int ntPnts, noffset;
762
763 m_LocTracephysToTraceIDMap = Array<OneD, Array<OneD, int>>(2);
764 m_LocTracephysToTraceIDMap[0] = Array<OneD, int>(m_nFwdLocTracePts, -1);
765 m_LocTracephysToTraceIDMap[1] =
766 Array<OneD, int>(m_nLocTracePts - m_nFwdLocTracePts, -1);
767
768 Array<OneD, NekDouble> tracePnts(m_nTracePts, 0.0);
769 for (int nt = 0; nt < ntotTrace; nt++)
770 {
771 ntPnts = tracelist->GetTotPoints(nt);
772 noffset = tracelist->GetPhys_Offset(nt);
773 for (int i = 0; i < ntPnts; i++)
774 {
775 tracePnts[noffset + i] = NekDouble(nt);
776 }
777 }
778
779 Array<OneD, Array<OneD, NekDouble>> loctracePntsLR(2);
780 loctracePntsLR[0] = Array<OneD, NekDouble>(m_nFwdLocTracePts, 0.0);
781 loctracePntsLR[1] =
782 Array<OneD, NekDouble>(m_nLocTracePts - m_nFwdLocTracePts, 0.0);
783
784 for (int dir = 0; dir < 2; dir++)
785 {
786 int cnt = 0;
787 int cnt1 = 0;
788
789 // tmp space assuming forward map is of size of trace
790 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
791 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(),
792 tracePnts.data(), m_locInterpTraceToTraceMap[dir].data(),
793 tmp.data());
794
795 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
796 {
797 if (m_interpNtraces[dir][i])
798 {
799 LibUtilities::PointsKey fromPointsKey0 =
800 std::get<0>(m_interpPoints[dir][i]);
801 LibUtilities::PointsKey fromPointsKey1 =
802 std::get<1>(m_interpPoints[dir][i]);
803 LibUtilities::PointsKey toPointsKey0 =
804 std::get<2>(m_interpPoints[dir][i]);
805 LibUtilities::PointsKey toPointsKey1 =
806 std::get<3>(m_interpPoints[dir][i]);
807
808 int fnp0 = fromPointsKey0.GetNumPoints();
809 int fnp1 = fromPointsKey1.GetNumPoints();
810 int tnp0 = toPointsKey0.GetNumPoints();
811 int tnp1 = toPointsKey1.GetNumPoints();
812
813 int nfttl = fnp0 * fnp1;
814
815 for (int ne = 0; ne < m_interpNtraces[dir][i]; ne++)
816 {
817 Vmath::Fill(nfttl, tmp[cnt1], &loctracePntsLR[dir][cnt], 1);
818 cnt += nfttl;
819 cnt1 += tnp0 * tnp1;
820 }
821 }
822 }
823 }
824
825 NekDouble error = 0.0;
826 for (int nlr = 0; nlr < 2; nlr++)
827 {
828 for (int i = 0; i < loctracePntsLR[nlr].size(); i++)
829 {
830 m_LocTracephysToTraceIDMap[nlr][i] =
831 std::round(loctracePntsLR[nlr][i]);
832 error += abs(loctracePntsLR[nlr][i] -
833 NekDouble(m_LocTracephysToTraceIDMap[nlr][i]));
834 }
835 }
836 error = error / NekDouble(m_nLocTracePts);
837 ASSERTL0(error < NekConstants::kNekZeroTol,
838 "m_LocTracephysToTraceIDMap may not be integer !!");
839}
840
841/**
842 * @brief Set up maps between coefficients on trace and in cells.
843 *
844 * @param locExp Expansion list in elements
845 * @param trace Expansion list on traces.
846 */
847void LocTraceToTraceMap::TraceLocToElmtLocCoeffMap(
848 const ExpList &locExp, const ExpListSharedPtr &trace)
849{
850 const std::shared_ptr<LocalRegions::ExpansionVector> exptrac =
851 trace->GetExp();
852 size_t ntrace = exptrac->size();
853
854 Array<OneD, Array<OneD, int>> LRAdjExpid{2};
855 Array<OneD, Array<OneD, bool>> LRAdjflag{2};
856
857 TensorOfArray3D<int> elmtLRMap{2};
858 TensorOfArray3D<int> elmtLRSign{2};
859
860 for (int lr = 0; lr < 2; ++lr)
861 {
862 LRAdjExpid[lr] = Array<OneD, int>{ntrace, 0};
863 LRAdjflag[lr] = Array<OneD, bool>{ntrace, false};
864 elmtLRMap[lr] = Array<OneD, Array<OneD, int>>{ntrace};
865 elmtLRSign[lr] = Array<OneD, Array<OneD, int>>{ntrace};
866 for (int i = 0; i < ntrace; ++i)
867 {
868 size_t ncoeff = trace->GetNcoeffs(i);
869 elmtLRMap[lr][i] = Array<OneD, int>{ncoeff, 0};
870 elmtLRSign[lr][i] = Array<OneD, int>{ncoeff, 0};
871 }
872 }
873
874 const Array<OneD, const pair<int, int>> field_coeffToElmt =
875 locExp.GetCoeffsToElmt();
876 const Array<OneD, const pair<int, int>> trace_coeffToElmt =
877 trace->GetCoeffsToElmt();
878
879 for (int lr = 0; lr < 2; ++lr)
880 {
881 int ntotcoeffs = m_nTraceCoeffs[lr];
882 for (int i = 0; i < ntotcoeffs; ++i)
883 {
884 int ncoeffField = m_traceCoeffsToElmtMap[lr][i];
885 int ncoeffTrace = m_traceCoeffsToElmtTrace[lr][i];
886 int sign = m_traceCoeffsToElmtSign[lr][i];
887
888 int ntraceelmt = trace_coeffToElmt[ncoeffTrace].first;
889 int ntracelocN = trace_coeffToElmt[ncoeffTrace].second;
890
891 int nfieldelmt = field_coeffToElmt[ncoeffField].first;
892 int nfieldlocN = field_coeffToElmt[ncoeffField].second;
893
894 LRAdjflag[lr][ntraceelmt] = true;
895 LRAdjExpid[lr][ntraceelmt] = nfieldelmt;
896
897 elmtLRMap[lr][ntraceelmt][ntracelocN] = nfieldlocN;
898 elmtLRSign[lr][ntraceelmt][ntracelocN] = sign;
899 }
900 }
901 m_leftRightAdjacentExpId = LRAdjExpid;
902 m_leftRightAdjacentExpFlag = LRAdjflag;
903 m_traceCoeffToLeftRightExpCoeffMap = elmtLRMap;
904 m_traceCoeffToLeftRightExpCoeffSign = elmtLRSign;
905}
906
907void LocTraceToTraceMap::FindElmtNeighbors(const ExpList &locExp,
908 const ExpListSharedPtr &trace)
909{
910 const std::shared_ptr<LocalRegions::ExpansionVector> exptrac =
911 trace->GetExp();
912 int ntrace = exptrac->size();
913
914 const std::shared_ptr<LocalRegions::ExpansionVector> exp = locExp.GetExp();
915 int nexp = exp->size();
916
917 Array<OneD, Array<OneD, int>> LRAdjExpid(2);
918 Array<OneD, Array<OneD, bool>> LRAdjflag(2);
919 LRAdjExpid = m_leftRightAdjacentExpId;
920 LRAdjflag = m_leftRightAdjacentExpFlag;
921
922 std::set<std::pair<int, int>> neighborSet;
923 int ntmp0, ntmp1;
924 for (int nt = 0; nt < ntrace; nt++)
925 {
926 if (LRAdjflag[0][nt] && LRAdjflag[1][nt])
927 {
928 ntmp0 = LRAdjExpid[0][nt];
929 ntmp1 = LRAdjExpid[1][nt];
930
931 ASSERTL0(ntmp0 != ntmp1,
932 " ntmp0==ntmp1, trace inside a element?? ");
933
934 std::set<std::pair<int, int>>::iterator it = neighborSet.begin();
935 neighborSet.insert(it, std::make_pair(ntmp0, ntmp1));
936 neighborSet.insert(it, std::make_pair(ntmp1, ntmp0));
937 }
938 }
939
940 Array<OneD, int> ElemIndex(nexp, 0);
941 for (std::set<std::pair<int, int>>::iterator it = neighborSet.begin();
942 it != neighborSet.end(); ++it)
943 {
944 int ncurrent = it->first;
945 ElemIndex[ncurrent]++;
946 }
947
948 Array<OneD, Array<OneD, int>> ElemNeighbsId(nexp);
949 Array<OneD, Array<OneD, int>> tmpId(nexp);
950 Array<OneD, int> ElemNeighbsNumb(nexp, -1);
951 Vmath::Vcopy(nexp, ElemIndex, 1, ElemNeighbsNumb, 1);
952 for (int ne = 0; ne < nexp; ne++)
953 {
954 int neighb = ElemNeighbsNumb[ne];
955 ElemNeighbsId[ne] = Array<OneD, int>(neighb, -1);
956 tmpId[ne] = Array<OneD, int>(neighb, -1);
957 }
958
959 for (int ne = 0; ne < nexp; ne++)
960 {
961 ElemIndex[ne] = 0;
962 }
963 for (std::set<std::pair<int, int>>::iterator it = neighborSet.begin();
964 it != neighborSet.end(); ++it)
965 {
966 int ncurrent = it->first;
967 int neighbor = it->second;
968 ElemNeighbsId[ncurrent][ElemIndex[ncurrent]] = neighbor;
969 ElemIndex[ncurrent]++;
970 }
971
972 // pickout repeated indexes
973 for (int ne = 0; ne < nexp; ne++)
974 {
975 ElemIndex[ne] = 0;
976 for (int nb = 0; nb < ElemNeighbsNumb[ne]; nb++)
977 {
978 int neighbId = ElemNeighbsId[ne][nb];
979 bool found = false;
980 for (int nc = 0; nc < ElemIndex[ne]; nc++)
981 {
982 if (ElemNeighbsId[ne][nb] == tmpId[ne][nc])
983 {
984 found = true;
985 }
986 }
987 if (!found)
988 {
989 tmpId[ne][ElemIndex[ne]] = neighbId;
990 ElemIndex[ne]++;
991 }
992 }
993 }
994 ElemNeighbsNumb = ElemIndex;
995 for (int ne = 0; ne < nexp; ne++)
996 {
997 int neighb = ElemNeighbsNumb[ne];
998 if (neighb > 0)
999 {
1000 ElemNeighbsId[ne] = Array<OneD, int>(neighb, -1);
1001 Vmath::Vcopy(neighb, tmpId[ne], 1, ElemNeighbsId[ne], 1);
1002 }
1003 }
1004
1005 // check errors
1006 for (int ne = 0; ne < nexp; ne++)
1007 {
1008 for (int nb = 0; nb < ElemNeighbsNumb[ne]; nb++)
1009 {
1010 ASSERTL0((ElemNeighbsId[ne][nb] >= 0) &&
1011 (ElemNeighbsId[ne][nb] <= nexp),
1012 "Element id <0 or >number of total elements")
1013 }
1014 }
1015
1016 m_ElemNeighbsNumb = ElemNeighbsNumb;
1017 m_ElemNeighbsId = ElemNeighbsId;
1018}
1019
1020/**
1021 * @brief Gather the local elemental traces in physical space from
1022 * field using #m_locTraceToFieldMap. Note traces are blocked together
1023 * in similar trace point ordering
1024 *
1025 * @param field Solution field in physical space
1026 * @param faces Resulting local traces.
1027 */
1028void LocTraceToTraceMap::LocTracesFromField(
1029 const Array<OneD, const NekDouble> &field, Array<OneD, NekDouble> faces)
1030{
1031 // The static cast is necessary because m_locTraceToFieldMap should be
1032 // Array<OneD, size_t> ... or at least the same type as
1033 // m_locTraceToFieldMap.size() ...
1034 Vmath::Gathr(static_cast<int>(m_locTraceToFieldMap.size()), field,
1035 m_locTraceToFieldMap, faces);
1036}
1037
1038/**
1039 * @brief Reverse process of LocTracesFromField()
1040 * Add the local traces in physical space to field using
1041 * #m_locTraceToFieldMap.
1042 *
1043 * @param field Solution field in physical space
1044 * @param faces local traces.
1045 */
1046void LocTraceToTraceMap::AddLocTracesToField(
1047 const Array<OneD, const NekDouble> &faces, Array<OneD, NekDouble> &field)
1048{
1049 size_t nfield = field.size();
1050 Array<OneD, NekDouble> tmp{nfield, 0.0};
1051 Vmath::Assmb(m_locTraceToFieldMap.size(), faces.data(),
1052 m_locTraceToFieldMap.data(), tmp.data());
1053 Vmath::Vadd(nfield, tmp, 1, field, 1, field, 1);
1054}
1055
1056/**
1057 * @brief Gather the forwards-oriented local traces in physical space from field
1058 * using #m_locTraceToFieldMap.
1059 *
1060 * @param field Solution field in physical space
1061 * @param faces Resulting local forwards-oriented traces.
1062 */
1063void LocTraceToTraceMap::FwdLocTracesFromField(
1064 const Array<OneD, const NekDouble> &field, Array<OneD, NekDouble> faces)
1065{
1066 Vmath::Gathr(m_nFwdLocTracePts, field, m_locTraceToFieldMap, faces);
1067}
1068
1069/**
1070 * @brief Reshuffle local elemental traces in physical space so that
1071 * similar faces points are blocked together so they can then be
1072 * interpolated with InterpLocTraceToTrace method.
1073 *
1074 * @param loctrace local traces in physical space
1075 * @param reshuffle traces ordered in reshuffled format of similar patterns.
1076 */
1077void LocTraceToTraceMap::ReshuffleLocTracesForInterp(
1078 const int dir, const Array<OneD, const NekDouble> &loctraces,
1079 Array<OneD, NekDouble> reshuffle)
1080{
1081 ASSERTL1(dir < 2, "option dir out of range, "
1082 " dir=0 is fwd, dir=1 is bwd");
1083
1084 Vmath::Gathr(static_cast<int>(m_locTraceToElmtTraceMap[dir].size()),
1085 loctraces, m_locTraceToElmtTraceMap[dir], reshuffle);
1086}
1087
1088/**
1089 * @brief Unshuffle local elemental traces in physical space from
1090 * similar faces points are blocked together to the local elemental trace format
1091 *
1092 * @param loctrace local traces in physical space
1093 * @param reshuffle traces ordered in reshuffled format of similar patterns.
1094 */
1095void LocTraceToTraceMap::UnshuffleLocTraces(
1096 const int dir, const Array<OneD, const NekDouble> &loctraces,
1097 Array<OneD, NekDouble> unshuffle)
1098{
1099 ASSERTL1(dir < 2, "option dir out of range, "
1100 " dir=0 is fwd, dir=1 is bwd");
1101
1102 if (m_locTraceToElmtTraceMap[dir].size()) // single elemt check
1103 {
1104 Vmath::Scatr(m_locTraceToElmtTraceMap[dir].size(), loctraces,
1105 m_locTraceToElmtTraceMap[dir], unshuffle);
1106 }
1107}
1108
1109void LocTraceToTraceMap::InterpLocTracesToTrace(
1110 const int dir, const Array<OneD, const NekDouble> &loctraces,
1111 Array<OneD, NekDouble> &traces)
1112{
1113 switch (m_expdim)
1114 {
1115 case 1: // Essentially do copy
1116 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(),
1117 loctraces.data(),
1118 m_locInterpTraceToTraceMap[dir].data(), traces.data());
1119 break;
1120 case 2:
1121 InterpLocEdgesToTrace(dir, loctraces, traces);
1122 break;
1123 case 3:
1124 InterpLocFacesToTrace(dir, loctraces, traces);
1125 break;
1126 default:
1127 NEKERROR(ErrorUtil::efatal, "Not set up");
1128 break;
1129 }
1130}
1131
1132/**
1133 * @brief Interpolate local trace edges to global trace edge point distributions
1134 * where required.
1135 *
1136 * @param dir Selects forwards (0) or backwards (1) direction.
1137 * @param locfaces Local trace edge storage.
1138 * @param faces Global trace edge storage
1139 */
1140void LocTraceToTraceMap::InterpLocEdgesToTrace(
1141 const int dir, const Array<OneD, const NekDouble> &locedges,
1142 Array<OneD, NekDouble> &edges)
1143{
1144 ASSERTL1(dir < 2, "option dir out of range, "
1145 " dir=0 is fwd, dir=1 is bwd");
1146
1147 int cnt = 0;
1148 int cnt1 = 0;
1149
1150 // tmp space assuming forward map is of size of trace
1151 Array<OneD, NekDouble> tmp(m_nTracePts);
1152
1153 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1154 {
1155 // Check if there are edges to interpolate
1156 if (m_interpNtraces[dir][i])
1157 {
1158 // Get to/from points
1159 LibUtilities::PointsKey fromPointsKey0 =
1160 std::get<0>(m_interpPoints[dir][i]);
1161 LibUtilities::PointsKey toPointsKey0 =
1162 std::get<2>(m_interpPoints[dir][i]);
1163
1164 int fnp = fromPointsKey0.GetNumPoints();
1165 int tnp = toPointsKey0.GetNumPoints();
1166 int nedges = m_interpNtraces[dir][i];
1167
1168 // Do interpolation here if required
1169 switch (m_interpTrace[dir][i])
1170 {
1171 case eNoInterp: // Just copy
1172 {
1173 Vmath::Vcopy(nedges * fnp, locedges.data() + cnt, 1,
1174 tmp.data() + cnt1, 1);
1175 }
1176 break;
1177 case eInterpDir0:
1178 {
1179 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1180 Blas::Dgemm('N', 'N', tnp, nedges, fnp, 1.0,
1181 I0->GetPtr().data(), tnp, locedges.data() + cnt,
1182 fnp, 0.0, tmp.data() + cnt1, tnp);
1183 }
1184 break;
1185 case eInterpEndPtDir0:
1186 {
1187 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1188
1189 for (int k = 0; k < nedges; ++k)
1190 {
1191 Vmath::Vcopy(fnp, &locedges[cnt + k * fnp], 1,
1192 &tmp[cnt1 + k * tnp], 1);
1193
1194 tmp[cnt1 + k * tnp + tnp - 1] = Blas::Ddot(
1195 fnp, locedges.data() + cnt + k * fnp, 1, &I0[0], 1);
1196 }
1197 }
1198 break;
1199 default:
1200 NEKERROR(ErrorUtil::efatal,
1201 "Invalid interpolation type for 2D elements");
1202 break;
1203 }
1204
1205 cnt += nedges * fnp;
1206 cnt1 += nedges * tnp;
1207 }
1208 }
1209
1210 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(), tmp.data(),
1211 m_locInterpTraceToTraceMap[dir].data(), edges.data());
1212}
1213
1214/**
1215 * @brief Interpolate local faces to trace face point distributions where
1216 * required.
1217 *
1218 * @param dir Selects forwards (0) or backwards (1) direction.
1219 * @param locfaces Local trace face storage.
1220 * @param faces Global trace face storage
1221 */
1222void LocTraceToTraceMap::InterpLocFacesToTrace(
1223 const int dir, const Array<OneD, const NekDouble> &locfaces,
1224 Array<OneD, NekDouble> faces)
1225{
1226 ASSERTL1(dir < 2, "option dir out of range, "
1227 " dir=0 is fwd, dir=1 is bwd");
1228
1229 int cnt1 = 0;
1230 int cnt = 0;
1231
1232 // tmp space assuming forward map is of size of trace
1233 Array<OneD, NekDouble> tmp(m_nTracePts);
1234
1235 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1236 {
1237 // Check if there are faces to interpolate
1238 if (m_interpNtraces[dir][i])
1239 {
1240 // Get to/from points
1241 LibUtilities::PointsKey fromPointsKey0 =
1242 std::get<0>(m_interpPoints[dir][i]);
1243 LibUtilities::PointsKey fromPointsKey1 =
1244 std::get<1>(m_interpPoints[dir][i]);
1245 LibUtilities::PointsKey toPointsKey0 =
1246 std::get<2>(m_interpPoints[dir][i]);
1247 LibUtilities::PointsKey toPointsKey1 =
1248 std::get<3>(m_interpPoints[dir][i]);
1249
1250 int fnp0 = fromPointsKey0.GetNumPoints();
1251 int fnp1 = fromPointsKey1.GetNumPoints();
1252 int tnp0 = toPointsKey0.GetNumPoints();
1253 int tnp1 = toPointsKey1.GetNumPoints();
1254 int nfaces = m_interpNtraces[dir][i];
1255 int nfromfacepts = nfaces * fnp0 * fnp1;
1256
1257 // Do interpolation here if required
1258 switch (m_interpTrace[dir][i])
1259 {
1260 case eNoInterp: // Just copy
1261 {
1262 Vmath::Vcopy(nfromfacepts, locfaces.data() + cnt, 1,
1263 tmp.data() + cnt1, 1);
1264 }
1265 break;
1266 case eInterpDir0:
1267 {
1268 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1269 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1270 I0->GetPtr().data(), tnp0,
1271 locfaces.data() + cnt, fnp0, 0.0,
1272 tmp.data() + cnt1, tnp0);
1273 }
1274 break;
1275 case eInterpEndPtDir0:
1276 {
1277 int nfaces = m_interpNtraces[dir][i];
1278 for (int k = 0; k < fnp0; ++k)
1279 {
1280 Vmath::Vcopy(nfaces * fnp1, locfaces.data() + cnt + k,
1281 fnp0, tmp.data() + cnt1 + k, tnp0);
1282 }
1283
1284 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1285 Blas::Dgemv('T', fnp0, tnp1 * nfaces, 1.0,
1286 tmp.data() + cnt1, tnp0, I0.data(), 1, 0.0,
1287 tmp.data() + cnt1 + tnp0 - 1, tnp0);
1288 }
1289 break;
1290 case eInterpDir1:
1291 {
1292 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1293 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1294 {
1295 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1296 locfaces.data() + cnt + j * fnp0 * fnp1,
1297 tnp0, I1->GetPtr().data(), tnp1, 0.0,
1298 tmp.data() + cnt1 + j * tnp0 * tnp1, tnp0);
1299 }
1300 }
1301 break;
1302 case eInterpEndPtDir1:
1303 {
1304 Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1305 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1306 {
1307 // copy all points
1308 Vmath::Vcopy(fnp0 * fnp1,
1309 locfaces.data() + cnt + j * fnp0 * fnp1, 1,
1310 tmp.data() + cnt1 + j * tnp0 * tnp1, 1);
1311
1312 // interpolate end points
1313 for (int k = 0; k < tnp0; ++k)
1314 {
1315 tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1316 Blas::Ddot(fnp1,
1317 locfaces.data() + cnt +
1318 j * fnp0 * fnp1 + k,
1319 fnp0, &I1[0], 1);
1320 }
1321 }
1322 }
1323 break;
1324 case eInterpBothDirs:
1325 {
1326 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1327 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1328 Array<OneD, NekDouble> wsp(m_interpNtraces[dir][i] * fnp0 *
1329 tnp1);
1330
1331 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1332 {
1333 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1334 locfaces.data() + cnt + j * fnp0 * fnp1,
1335 fnp0, I1->GetPtr().data(), tnp1, 0.0,
1336 wsp.data() + j * fnp0 * tnp1, fnp0);
1337 }
1338 Blas::Dgemm('N', 'N', tnp0, tnp1 * m_interpNtraces[dir][i],
1339 fnp0, 1.0, I0->GetPtr().data(), tnp0,
1340 wsp.data(), fnp0, 0.0, tmp.data() + cnt1, tnp0);
1341 }
1342 break;
1343 case eInterpEndPtDir0InterpDir1:
1344 {
1345 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1346 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1347
1348 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1349 {
1350 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1351 locfaces.data() + cnt + j * fnp0 * fnp1,
1352 fnp0, I1->GetPtr().data(), tnp1, 0.0,
1353 tmp.data() + cnt1 + j * tnp0 * tnp1, tnp0);
1354 }
1355
1356 Blas::Dgemv('T', fnp0, tnp1 * m_interpNtraces[dir][i], 1.0,
1357 tmp.data() + cnt1, tnp0, I0.data(), 1, 0.0,
1358 tmp.data() + cnt1 + tnp0 - 1, tnp0);
1359 }
1360 break;
1361 default:
1362 ASSERTL0(false, "Interplation case needs implementing");
1363 break;
1364 }
1365 cnt += nfromfacepts;
1366 cnt1 += m_interpNtraces[dir][i] * tnp0 * tnp1;
1367 }
1368 }
1369
1370 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(), tmp.data(),
1371 m_locInterpTraceToTraceMap[dir].data(), faces.data());
1372}
1373
1374void LocTraceToTraceMap::InterpLocTracesToTraceTranspose(
1375 const int dir, const Array<OneD, const NekDouble> &trace,
1376 Array<OneD, NekDouble> &loctrace)
1377{
1378 switch (m_expdim)
1379 {
1380 case 2:
1381 InterpLocEdgesToTraceTranspose(dir, trace, loctrace);
1382 break;
1383 case 3:
1384 InterpLocFacesToTraceTranspose(dir, trace, loctrace);
1385 break;
1386 default:
1387 NEKERROR(ErrorUtil::efatal, "Not set up");
1388 break;
1389 }
1390}
1391
1392/**
1393 * @brief Transpose of Interp local edges to trace
1394 *
1395 * @param dir Selects forwards (0) or backwards (1) direction.
1396 * @param edges Global trace edge
1397 * @param locedges Local trace edge
1398 */
1399void LocTraceToTraceMap::InterpLocEdgesToTraceTranspose(
1400 const int dir, const Array<OneD, const NekDouble> &edges,
1401 Array<OneD, NekDouble> &locedges)
1402{
1403 ASSERTL1(dir < 2, "option dir out of range, "
1404 " dir=0 is fwd, dir=1 is bwd");
1405
1406 int cnt = 0;
1407 int cnt1 = 0;
1408
1409 // tmp space assuming forward map is of size of trace
1410 Array<OneD, NekDouble> tmp{size_t(m_nTracePts)};
1411 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(), edges.data(),
1412 m_locInterpTraceToTraceMap[dir].data(), tmp.data());
1413
1414 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1415 {
1416 // Check if there are edges to interpolate
1417 if (m_interpNtraces[dir][i])
1418 {
1419 // Get to/from points
1420 LibUtilities::PointsKey fromPointsKey0 =
1421 std::get<0>(m_interpPoints[dir][i]);
1422 LibUtilities::PointsKey toPointsKey0 =
1423 std::get<2>(m_interpPoints[dir][i]);
1424
1425 int fnp = fromPointsKey0.GetNumPoints();
1426 int tnp = toPointsKey0.GetNumPoints();
1427 int nedges = m_interpNtraces[dir][i];
1428
1429 // Do interpolation here if required
1430 switch (m_interpTrace[dir][i])
1431 {
1432 case eNoInterp: // Just copy
1433 {
1434 Vmath::Vcopy(nedges * fnp, tmp.data() + cnt1, 1,
1435 locedges.data() + cnt, 1);
1436 }
1437 break;
1438 case eInterpDir0:
1439 {
1440 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1441 Blas::Dgemm('T', 'N', fnp, nedges, tnp, 1.0,
1442 I0->GetPtr().data(), tnp, tmp.data() + cnt1,
1443 tnp, 0.0, locedges.data() + cnt, fnp);
1444 }
1445 break;
1446 case eInterpEndPtDir0:
1447 {
1448 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1449
1450 for (int k = 0; k < nedges; ++k)
1451 {
1452 Vmath::Vcopy(fnp, &tmp[cnt1 + k * tnp], 1,
1453 &locedges[cnt + k * fnp], 1);
1454
1455 Vmath::Svtvp(fnp, tmp[cnt1 + k * tnp + tnp - 1], &I0[0],
1456 1, locedges.data() + cnt + k * fnp, 1,
1457 locedges.data() + cnt + k * fnp, 1);
1458 }
1459 }
1460 break;
1461 default:
1462 NEKERROR(ErrorUtil::efatal,
1463 "Invalid interpolation type for 2D elements");
1464 break;
1465 }
1466
1467 cnt += nedges * fnp;
1468 cnt1 += nedges * tnp;
1469 }
1470 }
1471}
1472
1473/**
1474 * @brief transpose of interp local faces to trace
1475 *
1476 * @param dir Selects forwards (0) or backwards (1) direction.
1477 * @param loctraces Local trace
1478 * @param traces trace .
1479 */
1480void LocTraceToTraceMap::InterpLocFacesToTraceTranspose(
1481 const int dir, const Array<OneD, const NekDouble> &traces,
1482 Array<OneD, NekDouble> &loctraces)
1483{
1484 ASSERTL1(dir < 2, "option dir out of range, "
1485 " dir=0 is fwd, dir=1 is bwd");
1486
1487 int cnt = 0;
1488 int cnt1 = 0;
1489
1490 // tmp space assuming forward map is of size of trace
1491 Array<OneD, NekDouble> tmp{size_t(m_nTracePts)};
1492 // The static cast is necessary because m_locInterpTraceToTraceMap should be
1493 // Array<OneD, size_t> ... or at least the same type as
1494 // m_locInterpTraceToTraceMap.size() ...
1495 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1496 traces.data(), m_locInterpTraceToTraceMap[dir].data(),
1497 tmp.data());
1498
1499 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1500 {
1501 // Check if there are elementboundaries to interpolate
1502 if (m_interpNtraces[dir][i])
1503 {
1504 // Get to/from points
1505 LibUtilities::PointsKey fromPointsKey0 =
1506 std::get<0>(m_interpPoints[dir][i]);
1507 LibUtilities::PointsKey fromPointsKey1 =
1508 std::get<1>(m_interpPoints[dir][i]);
1509 LibUtilities::PointsKey toPointsKey0 =
1510 std::get<2>(m_interpPoints[dir][i]);
1511 LibUtilities::PointsKey toPointsKey1 =
1512 std::get<3>(m_interpPoints[dir][i]);
1513 // Here the f(from) and t(to) are chosen to be consistent with
1514 // InterpLocFacesToTrace
1515 int fnp0 = fromPointsKey0.GetNumPoints();
1516 int fnp1 = fromPointsKey1.GetNumPoints();
1517 int tnp0 = toPointsKey0.GetNumPoints();
1518 int tnp1 = toPointsKey1.GetNumPoints();
1519 int nfromfacepts = m_interpNtraces[dir][i] * fnp0 * fnp1;
1520
1521 // Do transpose interpolation here if required
1522 switch (m_interpTrace[dir][i])
1523 {
1524 case eNoInterp: // Just copy
1525 {
1526 Vmath::Vcopy(nfromfacepts, tmp.data() + cnt1, 1,
1527 loctraces.data() + cnt, 1);
1528 }
1529 break;
1530 case eInterpDir0:
1531 {
1532 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1533 Blas::Dgemm('T', 'N', fnp0, tnp1, tnp0, 1.0,
1534 I0->GetPtr().data(), tnp0, tmp.data() + cnt1,
1535 tnp0, 0.0, loctraces.data() + cnt, fnp0);
1536 }
1537 break;
1538 case eInterpEndPtDir0:
1539 {
1540 int nfaces = m_interpNtraces[dir][i];
1541 for (int k = 0; k < fnp0; ++k)
1542 {
1543 Vmath::Vcopy(nfaces * fnp1, tmp.data() + cnt1 + k, tnp0,
1544 loctraces.data() + cnt + k, fnp0);
1545 }
1546
1547 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1548 for (int k = 0; k < tnp1 * nfaces; k++)
1549 {
1550 Vmath::Svtvp(fnp0, tmp[cnt1 + tnp0 - 1 + k * tnp0],
1551 &I0[0], 1,
1552 loctraces.data() + cnt + k * fnp0, 1,
1553 loctraces.data() + cnt + k * fnp0, 1);
1554 }
1555 }
1556 break;
1557 case eInterpDir1:
1558 {
1559 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1560
1561 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1562 {
1563 Blas::Dgemm('N', 'N', tnp0, fnp1, tnp1, 1.0,
1564 tmp.data() + cnt1 + j * tnp0 * tnp1, tnp0,
1565 I1->GetPtr().data(), tnp1, 0.0,
1566 loctraces.data() + cnt + j * fnp0 * fnp1,
1567 tnp0);
1568 }
1569 }
1570 break;
1571 case eInterpEndPtDir1:
1572 {
1573 Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1574 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1575 {
1576 Vmath::Vcopy(
1577 fnp0 * fnp1, tmp.data() + cnt1 + j * tnp0 * tnp1, 1,
1578 loctraces.data() + cnt + j * fnp0 * fnp1, 1);
1579
1580 for (int k = 0; k < fnp1; k++)
1581 {
1582 Vmath::Svtvp(
1583 fnp0, I1[k],
1584 &tmp[cnt1 + (j + 1) * tnp0 * tnp1 - tnp0], 1,
1585 &loctraces[cnt + j * fnp0 * fnp1 + k * fnp0], 1,
1586 &loctraces[cnt + j * fnp0 * fnp1 + k * fnp0],
1587 1);
1588 }
1589 }
1590 }
1591 break;
1592 case eInterpBothDirs:
1593 {
1594 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1595 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1596
1597 Array<OneD, NekDouble> wsp{
1598 size_t(m_interpNtraces[dir][i] * fnp0 * tnp1)};
1599
1600 Blas::Dgemm('T', 'N', fnp0, tnp1 * m_interpNtraces[dir][i],
1601 tnp0, 1.0, I0->GetPtr().data(), tnp0,
1602 tmp.data() + cnt1, tnp0, 0.0, wsp.data(), fnp0);
1603
1604 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1605 {
1606 Blas::Dgemm('N', 'N', fnp0, fnp1, tnp1, 1.0,
1607 wsp.data() + j * fnp0 * tnp1, fnp0,
1608 I1->GetPtr().data(), tnp1, 0.0,
1609 loctraces.data() + cnt + j * fnp0 * fnp1,
1610 fnp0);
1611 }
1612 }
1613 break;
1614 case eInterpEndPtDir0InterpDir1:
1615 {
1616 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1617 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1618
1619 Array<OneD, NekDouble> wsp{
1620 size_t(m_interpNtraces[dir][i] * fnp0 * tnp1)};
1621
1622 for (int k = 0; k < tnp1 * m_interpNtraces[dir][i]; k++)
1623 {
1624 Vmath::Svtvp(fnp0, tmp[cnt1 + tnp0 - 1 + k * tnp0],
1625 &I0[0], 1, tmp.data() + cnt1 + k * tnp0, 1,
1626 wsp.data() + k * fnp0, 1);
1627 }
1628
1629 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1630 {
1631 Blas::Dgemm('N', 'N', fnp0, fnp1, tnp1, 1.0,
1632 wsp.data() + j * fnp0 * tnp1, fnp0,
1633 I1->GetPtr().data(), tnp1, 0.0,
1634 loctraces.data() + cnt + j * fnp0 * fnp1,
1635 fnp0);
1636 }
1637 }
1638 break;
1639 }
1640 cnt += nfromfacepts;
1641 cnt1 += m_interpNtraces[dir][i] * tnp0 * tnp1;
1642 }
1643 }
1644}
1645
1646void LocTraceToTraceMap::InterpTraceToLocTrace(
1647 const int dir, const Array<OneD, NekDouble> &traces,
1648 Array<OneD, NekDouble> &loctraces)
1649{
1650 switch (m_expdim)
1651 {
1652 case 1: // Essentially do copy
1653 Vmath::Gathr(
1654 static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1655 traces.data(), m_locInterpTraceToTraceMap[dir].data(),
1656 loctraces.data());
1657 break;
1658 case 2:
1659 InterpTraceToLocEdges(dir, traces, loctraces);
1660 break;
1661 case 3:
1662 InterpTraceToLocFaces(dir, traces, loctraces);
1663 break;
1664 default:
1665 ASSERTL0(false, "Not set up");
1666 break;
1667 }
1668}
1669
1670/**
1671 * @brief Interpolate global trace edge to local trace edges point
1672 * distributions where required.
1673 *
1674 * @param dir Selects forwards (0) or backwards (1) direction.
1675 * @param locfaces Local trace edge storage.
1676 * @param faces Global trace edge storage
1677 */
1678void LocTraceToTraceMap::InterpTraceToLocEdges(
1679 const int dir, const Array<OneD, const NekDouble> &edges,
1680 Array<OneD, NekDouble> &locedges)
1681{
1682 ASSERTL1(dir < 2, "option dir out of range, "
1683 " dir=0 is fwd, dir=1 is bwd");
1684
1685 int cnt = 0;
1686 int cnt1 = 0;
1687
1688 Array<OneD, NekDouble> tmp(m_nTracePts);
1689
1690 // unshuffles trace into lcoally orientated format.
1691 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1692 edges.data(), m_locInterpTraceToTraceMap[dir].data(),
1693 tmp.data());
1694
1695 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1696 {
1697 // Check if there are edges to interpolate
1698 if (m_interpNtraces[dir][i])
1699 {
1700 // Get to/from points
1701 LibUtilities::PointsKey fromPointsKey0 =
1702 std::get<2>(m_interpPoints[dir][i]);
1703 LibUtilities::PointsKey toPointsKey0 =
1704 std::get<0>(m_interpPoints[dir][i]);
1705
1706 int fnp = fromPointsKey0.GetNumPoints();
1707 int tnp = toPointsKey0.GetNumPoints();
1708 int nedges = m_interpNtraces[dir][i];
1709
1710 // Do interpolation here if required
1711 switch (m_interpTrace[dir][i])
1712 {
1713 case eNoInterp: // Just copy
1714 {
1715 Vmath::Vcopy(nedges * fnp, tmp.data() + cnt, 1,
1716 locedges.data() + cnt1, 1);
1717 }
1718 break;
1719 case eInterpDir0:
1720 {
1721 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
1722 Blas::Dgemm('N', 'N', tnp, nedges, fnp, 1.0,
1723 I0->GetPtr().data(), tnp, tmp.data() + cnt, fnp,
1724 0.0, locedges.data() + cnt1, tnp);
1725 }
1726 break;
1727 case eInterpEndPtDir0:
1728 {
1729 // Just copy points back
1730 for (int k = 0; k < nedges; ++k)
1731 {
1732 // Should be tnp rather than fnp on first argument.
1733 Vmath::Vcopy(tnp, &tmp[cnt + k * fnp], 1,
1734 &locedges[cnt1 + k * tnp], 1);
1735 }
1736 }
1737 break;
1738 default:
1739 ASSERTL0(false,
1740 "Invalid interpolation type for 2D elements");
1741 break;
1742 }
1743
1744 cnt += nedges * fnp;
1745 cnt1 += nedges * tnp;
1746 }
1747 }
1748}
1749
1750/**
1751 * @brief Interpolate global trace edge to local trace edges point
1752 * distributions where required.
1753 *
1754 * @param dir Selects forwards (0) or backwards (1) direction.
1755 * @param locfaces Local trace edge storage.
1756 * @param faces Global trace edge storage
1757 */
1758void LocTraceToTraceMap::InterpTraceToLocFaces(
1759 const int dir, const Array<OneD, const NekDouble> &faces,
1760 Array<OneD, NekDouble> &locfaces)
1761{
1762 ASSERTL1(dir < 2, "option dir out of range, "
1763 " dir=0 is fwd, dir=1 is bwd");
1764
1765 int cnt = 0;
1766 int cnt1 = 0;
1767
1768 Array<OneD, NekDouble> tmp(m_nTracePts);
1769
1770 // unshuffles trace into lcoally orientated format.
1771 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1772 faces.data(), m_locInterpTraceToTraceMap[dir].data(),
1773 tmp.data());
1774
1775 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1776 {
1777 // Check if there are faces to interpolate
1778 if (m_interpNtraces[dir][i])
1779 {
1780 // Get to/from points
1781 LibUtilities::PointsKey fromPointsKey0 =
1782 std::get<2>(m_interpPoints[dir][i]);
1783 LibUtilities::PointsKey fromPointsKey1 =
1784 std::get<3>(m_interpPoints[dir][i]);
1785 LibUtilities::PointsKey toPointsKey0 =
1786 std::get<0>(m_interpPoints[dir][i]);
1787 LibUtilities::PointsKey toPointsKey1 =
1788 std::get<1>(m_interpPoints[dir][i]);
1789
1790 int fnp0 = fromPointsKey0.GetNumPoints();
1791 int fnp1 = fromPointsKey1.GetNumPoints();
1792 int tnp0 = toPointsKey0.GetNumPoints();
1793 int tnp1 = toPointsKey1.GetNumPoints();
1794 int nfaces = m_interpNtraces[dir][i];
1795 int nfromfacepts = nfaces * fnp0 * fnp1;
1796
1797 // Do interpolation here if required
1798 switch (m_interpTrace[dir][i])
1799 {
1800 case eNoInterp: // Just copy
1801 {
1802 Vmath::Vcopy(nfromfacepts, tmp.data() + cnt, 1,
1803 locfaces.data() + cnt1, 1);
1804 }
1805 break;
1806 case eInterpDir0:
1807 {
1808 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
1809 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1810 I0->GetPtr().data(), tnp0, tmp.data() + cnt,
1811 fnp0, 0.0, locfaces.data() + cnt1, tnp0);
1812 }
1813 break;
1814 case eInterpDir1:
1815 {
1816 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
1817 for (int j = 0; j < nfaces; ++j)
1818 {
1819 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1820 tmp.data() + cnt + j * fnp0 * fnp1, tnp0,
1821 I1->GetPtr().data(), tnp1, 0.0,
1822 locfaces.data() + cnt1 + j * tnp0 * tnp1,
1823 tnp0);
1824 }
1825 }
1826 break;
1827 case eInterpEndPtDir1:
1828 {
1829 for (int j = 0; j < nfaces; ++j)
1830 {
1831 // copy all points missing off top verex in dir 1
1832 Vmath::Vcopy(
1833 tnp0 * tnp1, tmp.data() + cnt + j * fnp0 * fnp1, 1,
1834 locfaces.data() + cnt1 + j * tnp0 * tnp1, 1);
1835 }
1836 }
1837 break;
1838 case eInterpBothDirs:
1839 {
1840 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
1841 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
1842 Array<OneD, NekDouble> wsp(nfaces * fnp0 * tnp1 * fnp0);
1843
1844 for (int j = 0; j < nfaces; ++j)
1845 {
1846 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1847 tmp.data() + cnt + j * fnp0 * fnp1, fnp0,
1848 I1->GetPtr().data(), tnp1, 0.0,
1849 wsp.data() + j * fnp0 * tnp1, fnp0);
1850 }
1851
1852 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1853 I0->GetPtr().data(), tnp0, wsp.data(), fnp0,
1854 0.0, locfaces.data() + cnt1, tnp0);
1855 }
1856 break;
1857 case eInterpEndPtDir0InterpDir1:
1858 {
1859 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
1860 for (int j = 0; j < nfaces; ++j)
1861 {
1862 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1863 tmp.data() + cnt + j * fnp0 * fnp1, fnp0,
1864 I1->GetPtr().data(), tnp1, 0.0,
1865 locfaces.data() + cnt1 + j * tnp0 * tnp1,
1866 tnp0);
1867 }
1868 }
1869 break;
1870 default:
1871 ASSERTL0(false, "Interpolation case not implemneted (yet)");
1872 break;
1873 }
1874 cnt += nfromfacepts;
1875 cnt1 += m_interpNtraces[dir][i] * tnp0 * tnp1;
1876 }
1877 }
1878}
1879
1880/**
1881 * @brief Add contributions from trace coefficients to the elemental field
1882 * storage.
1883 *
1884 * @param trace Array of global trace coefficients.
1885 * @param field Array containing field coefficients storage.
1886 */
1887void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1888 const Array<OneD, const NekDouble> &trace, Array<OneD, NekDouble> &field)
1889{
1890 int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
1891 for (int i = 0; i < nvals; ++i)
1892 {
1893 field[m_traceCoeffsToElmtMap[0][i]] +=
1894 m_traceCoeffsToElmtSign[0][i] *
1895 trace[m_traceCoeffsToElmtTrace[0][i]];
1896 }
1897}
1898
1899/**
1900 * @brief Add contributions from backwards or forwards oriented trace
1901 * coefficients to the elemental field storage.
1902 *
1903 * @param dir Selects forwards (0) or backwards (1) direction
1904 * @param trace Array of global trace coefficients.
1905 * @param field Array containing field coefficients storage.
1906 */
1907void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1908 const int dir, const Array<OneD, const NekDouble> &trace,
1909 Array<OneD, NekDouble> &field)
1910{
1911 int nvals = m_nTraceCoeffs[dir];
1912 for (int i = 0; i < nvals; ++i)
1913 {
1914 field[m_traceCoeffsToElmtMap[dir][i]] +=
1915 m_traceCoeffsToElmtSign[dir][i] *
1916 trace[m_traceCoeffsToElmtTrace[dir][i]];
1917 }
1918}
1919
1920} // namespace Nektar::MultiRegions
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:50
Nektar::LibUtilities::PointsKey::GetPointsType
PointsType GetPointsType() const
Definition: Points.h:90
Nektar::LibUtilities::PointsKey::GetNumPoints
size_t GetNumPoints() const
Definition: Points.h:85
Nektar::MultiRegions::ExpList
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:99
Nektar::MultiRegions::ExpList::GetCoeffsToElmt
const Array< OneD, const std::pair< int, int > > & GetCoeffsToElmt() const
Get m_coeffs to elemental value map.
Definition: ExpList.h:2118
Nektar::MultiRegions::ExpList::GetCoeff_Offset
int GetCoeff_Offset(int n) const
Get the start offset position for a local contiguous list of coeffs correspoinding to element n.
Definition: ExpList.h:2078
Nektar::MultiRegions::ExpList::GetExp
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2070
Nektar::MultiRegions::ExpList::GetPhys_Offset
int GetPhys_Offset(int n) const
Get the start offset position for a local contiguous list of quadrature points in a full array corres...
Definition: ExpList.h:2085
Nektar::MultiRegions::LocTraceToTraceMap::m_traceCoeffToLeftRightExpCoeffMap
Array< OneD, Array< OneD, Array< OneD, int > > > m_traceCoeffToLeftRightExpCoeffMap
The map of every coeff from current trace to the left & right adjacent expasion coeffs.
Definition: LocTraceToTraceMap.h:389
Nektar::MultiRegions::LocTraceToTraceMap::m_traceCoeffsToElmtSign
Array< OneD, Array< OneD, int > > m_traceCoeffsToElmtSign
Sign array for mapping from forwards/backwards trace coefficients to local trace storage.
Definition: LocTraceToTraceMap.h:380
Nektar::MultiRegions::LocTraceToTraceMap::m_interpNtraces
Array< OneD, Array< OneD, int > > m_interpNtraces
Number of edges/faces on a 2D/3D element that require interpolation.
Definition: LocTraceToTraceMap.h:369
Nektar::MultiRegions::LocTraceToTraceMap::m_locTraceToElmtTraceMap
Array< OneD, Array< OneD, int > > m_locTraceToElmtTraceMap
A mapping from the local elemental trace points, arranged as all forwards traces followed by backward...
Definition: LocTraceToTraceMap.h:337
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocEdgesToTraceTranspose
void InterpLocEdgesToTraceTranspose(const int dir, const Array< OneD, const NekDouble > &edges, Array< OneD, NekDouble > &locedges)
Transpose of interp local edges to Trace methods.
Definition: LocTraceToTraceMap.cpp:1399
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocTracesToTraceTranspose
void InterpLocTracesToTraceTranspose(const int dir, const Array< OneD, const NekDouble > &traces, Array< OneD, NekDouble > &loctraces)
Definition: LocTraceToTraceMap.cpp:1374
Nektar::MultiRegions::LocTraceToTraceMap::m_traceCoeffsToElmtTrace
Array< OneD, Array< OneD, int > > m_traceCoeffsToElmtTrace
Mapping from forwards/backwards trace coefficients to the position of the trace element in global sto...
Definition: LocTraceToTraceMap.h:377
Nektar::MultiRegions::LocTraceToTraceMap::m_leftRightAdjacentExpId
Array< OneD, Array< OneD, int > > m_leftRightAdjacentExpId
The expansion ID that are the left & right adjacent to current trace.
Definition: LocTraceToTraceMap.h:385
Nektar::MultiRegions::LocTraceToTraceMap::LocTraceToTraceMap
LocTraceToTraceMap(const ExpList &locExp, const ExpListSharedPtr &trace, const Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > &elmtToTrace, const std::vector< bool > &LeftAdjacents)
Set up trace to trace mapping components.
Definition: LocTraceToTraceMap.cpp:60
Nektar::MultiRegions::LocTraceToTraceMap::m_interpFromTraceI1
Array< OneD, Array< OneD, DNekMatSharedPtr > > m_interpFromTraceI1
Interpolation matrices for either 2D edges or first coordinate of 3D face using going "from' to 'to' ...
Definition: LocTraceToTraceMap.h:357
Nektar::MultiRegions::LocTraceToTraceMap::m_interpEndPtI0
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_interpEndPtI0
Mapping to hold first coordinate direction endpoint interpolation, which can be more optimal if using...
Definition: LocTraceToTraceMap.h:363
Nektar::MultiRegions::LocTraceToTraceMap::LocTracesFromField
void LocTracesFromField(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > faces)
Gather the local elemental traces in physical space from field using m_locTraceToFieldMap....
Definition: LocTraceToTraceMap.cpp:1028
Nektar::MultiRegions::LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_3D
void CalcLocTracePhysToTraceIDMap_3D(const ExpListSharedPtr &tracelist)
Definition: LocTraceToTraceMap.cpp:755
Nektar::MultiRegions::LocTraceToTraceMap::m_nLocTracePts
int m_nLocTracePts
The number of local trace points.
Definition: LocTraceToTraceMap.h:327
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocEdgesToTrace
void InterpLocEdgesToTrace(const int dir, const Array< OneD, const NekDouble > &locfaces, Array< OneD, NekDouble > &edges)
Interpolate local trace edges to global trace edge point distributions where required.
Definition: LocTraceToTraceMap.cpp:1140
Nektar::MultiRegions::LocTraceToTraceMap::m_nTracePts
int m_nTracePts
The number of global trace points.
Definition: LocTraceToTraceMap.h:329
Nektar::MultiRegions::LocTraceToTraceMap::UnshuffleLocTraces
void UnshuffleLocTraces(const int dir, const Array< OneD, const NekDouble > &loctraces, Array< OneD, NekDouble > unshuffle)
Unshuffle local elemental traces in physical space from similar faces points are blocked together to ...
Definition: LocTraceToTraceMap.cpp:1095
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocTracesToTrace
void InterpLocTracesToTrace(const int dir, const Array< OneD, const NekDouble > &loctraces, Array< OneD, NekDouble > &traces)
Definition: LocTraceToTraceMap.cpp:1109
Nektar::MultiRegions::LocTraceToTraceMap::m_leftRightAdjacentExpFlag
Array< OneD, Array< OneD, bool > > m_leftRightAdjacentExpFlag
Flag indicates whether the expansion that are the left & right adjacent to current trace exists.
Definition: LocTraceToTraceMap.h:383
Nektar::MultiRegions::LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap
void CalcLocTracePhysToTraceIDMap(const ExpListSharedPtr &tracelist, const int ndim)
Definition: LocTraceToTraceMap.cpp:660
Nektar::MultiRegions::LocTraceToTraceMap::AddLocTracesToField
void AddLocTracesToField(const Array< OneD, const NekDouble > &faces, Array< OneD, NekDouble > &field)
Reverse process of LocTracesFromField() Add the local traces in physical space to field using m_locTr...
Definition: LocTraceToTraceMap.cpp:1046
Nektar::MultiRegions::LocTraceToTraceMap::m_nTraceCoeffs
int m_nTraceCoeffs[2]
Number of forwards/backwards trace coefficients.
Definition: LocTraceToTraceMap.h:371
Nektar::MultiRegions::LocTraceToTraceMap::m_nFwdLocTracePts
int m_nFwdLocTracePts
The number of forward trace points. A local trace element is ‘forward’ if it is the side selected for...
Definition: LocTraceToTraceMap.h:325
Nektar::MultiRegions::LocTraceToTraceMap::m_traceCoeffToLeftRightExpCoeffSign
Array< OneD, Array< OneD, Array< OneD, int > > > m_traceCoeffToLeftRightExpCoeffSign
The sign of every coeff from current trace to the left & right adjacent expasion coeffs.
Definition: LocTraceToTraceMap.h:393
Nektar::MultiRegions::LocTraceToTraceMap::InterpTraceToLocFaces
void InterpTraceToLocFaces(const int dir, const Array< OneD, const NekDouble > &faces, Array< OneD, NekDouble > &locfaces)
Interpolate global trace edge to local trace edges point distributions where required.
Definition: LocTraceToTraceMap.cpp:1758
Nektar::MultiRegions::LocTraceToTraceMap::m_interpEndPtI1
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_interpEndPtI1
Mapping to hold second coordinate direction endpoint interpolation, which can be more optimal if usin...
Definition: LocTraceToTraceMap.h:367
Nektar::MultiRegions::LocTraceToTraceMap::TraceLocToElmtLocCoeffMap
void TraceLocToElmtLocCoeffMap(const ExpList &locExp, const ExpListSharedPtr &trace)
Set up maps between coefficients on trace and in cells.
Definition: LocTraceToTraceMap.cpp:847
Nektar::MultiRegions::LocTraceToTraceMap::m_traceCoeffsToElmtMap
Array< OneD, Array< OneD, int > > m_traceCoeffsToElmtMap
Mapping from forwards/backwards trace coefficients to elemental coefficient storage.
Definition: LocTraceToTraceMap.h:374
Nektar::MultiRegions::LocTraceToTraceMap::m_expdim
int m_expdim
Expansion Dimension we have setup for trace mapping.
Definition: LocTraceToTraceMap.h:322
Nektar::MultiRegions::LocTraceToTraceMap::m_interpTrace
Array< OneD, Array< OneD, InterpLocTraceToTrace > > m_interpTrace
A mapping holding the type of interpolation needed for each local trace. Dimension 0 holds forward tr...
Definition: LocTraceToTraceMap.h:343
Nektar::MultiRegions::LocTraceToTraceMap::InterpTraceToLocEdges
void InterpTraceToLocEdges(const int dir, const Array< OneD, const NekDouble > &locfaces, Array< OneD, NekDouble > &edges)
Interpolate global trace edge to local trace edges point distributions where required.
Definition: LocTraceToTraceMap.cpp:1678
Nektar::MultiRegions::LocTraceToTraceMap::FwdLocTracesFromField
void FwdLocTracesFromField(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > faces)
Gather the forwards-oriented local traces in physical space from field using m_locTraceToFieldMap.
Definition: LocTraceToTraceMap.cpp:1063
Nektar::MultiRegions::LocTraceToTraceMap::Setup
void Setup(const ExpList &locExp, const ExpListSharedPtr &trace, const Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > &elmtToTrace, const std::vector< bool > &LeftAdjacents)
Set up member variables for a two-dimensional problem.
Definition: LocTraceToTraceMap.cpp:87
Nektar::MultiRegions::LocTraceToTraceMap::m_locInterpTraceToTraceMap
Array< OneD, Array< OneD, int > > m_locInterpTraceToTraceMap
A mapping from local trace points to the global trace. Dimension 0 holds forward traces,...
Definition: LocTraceToTraceMap.h:340
Nektar::MultiRegions::LocTraceToTraceMap::m_ElemNeighbsId
Array< OneD, Array< OneD, int > > m_ElemNeighbsId
Definition: LocTraceToTraceMap.h:404
Nektar::MultiRegions::LocTraceToTraceMap::m_interpPoints
Array< OneD, Array< OneD, TraceInterpPoints > > m_interpPoints
Interpolation points key distributions to each of the local to global mappings.
Definition: LocTraceToTraceMap.h:360
Nektar::MultiRegions::LocTraceToTraceMap::m_locTraceToFieldMap
Array< OneD, int > m_locTraceToFieldMap
A mapping from the local elemental trace points, arranged as all forwards traces followed by backward...
Definition: LocTraceToTraceMap.h:333
Nektar::MultiRegions::LocTraceToTraceMap::m_interpFromTraceI0
Array< OneD, Array< OneD, DNekMatSharedPtr > > m_interpFromTraceI0
Interpolation matrices for either 2D edges or first coordinate of 3D face using going "from' to 'to' ...
Definition: LocTraceToTraceMap.h:353
Nektar::MultiRegions::LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs
void AddTraceCoeffsToFieldCoeffs(const Array< OneD, const NekDouble > &trace, Array< OneD, NekDouble > &field)
Add contributions from trace coefficients to the elemental field storage.
Definition: LocTraceToTraceMap.cpp:1887
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocFacesToTraceTranspose
void InterpLocFacesToTraceTranspose(const int dir, const Array< OneD, const NekDouble > &traces, Array< OneD, NekDouble > &loctraces)
Transpose of interp local faces to Trace methods.
Definition: LocTraceToTraceMap.cpp:1480
Nektar::MultiRegions::LocTraceToTraceMap::m_interpTraceI1
Array< OneD, Array< OneD, DNekMatSharedPtr > > m_interpTraceI1
Interpolation matrices for the second coordinate of 3D face, not used in 2D.
Definition: LocTraceToTraceMap.h:349
Nektar::MultiRegions::LocTraceToTraceMap::m_interpTraceI0
Array< OneD, Array< OneD, DNekMatSharedPtr > > m_interpTraceI0
Interpolation matrices for either 2D edges or first coordinate of 3D face.
Definition: LocTraceToTraceMap.h:346
Nektar::MultiRegions::LocTraceToTraceMap::FindElmtNeighbors
void FindElmtNeighbors(const ExpList &locExp, const ExpListSharedPtr &trace)
Definition: LocTraceToTraceMap.cpp:907
Nektar::MultiRegions::LocTraceToTraceMap::InterpTraceToLocTrace
void InterpTraceToLocTrace(const int dir, const Array< OneD, NekDouble > &traces, Array< OneD, NekDouble > &loctraces)
Definition: LocTraceToTraceMap.cpp:1646
Nektar::MultiRegions::LocTraceToTraceMap::ReshuffleLocTracesForInterp
void ReshuffleLocTracesForInterp(const int dir, const Array< OneD, const NekDouble > &loctraces, Array< OneD, NekDouble > reshuffle)
Reshuffle local elemental traces in physical space so that similar faces points are blocked together ...
Definition: LocTraceToTraceMap.cpp:1077
Nektar::MultiRegions::LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_2D
void CalcLocTracePhysToTraceIDMap_2D(const ExpListSharedPtr &tracelist)
Definition: LocTraceToTraceMap.cpp:677
Nektar::MultiRegions::LocTraceToTraceMap::m_LocTracephysToTraceIDMap
Array< OneD, Array< OneD, int > > m_LocTracephysToTraceIDMap
Definition: LocTraceToTraceMap.h:406
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocFacesToTrace
void InterpLocFacesToTrace(const int dir, const Array< OneD, const NekDouble > &locfaces, Array< OneD, NekDouble > faces)
Interpolate local faces to trace face point distributions where required.
Definition: LocTraceToTraceMap.cpp:1222
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = alpha A x plus beta y where A[m x n].
Definition: Blas.hpp:211
Blas::Ddot
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
FilterPython_Function.field
field
Definition: FilterPython_Function.py:4
Nektar::LibUtilities::PointsManager
PointsManagerT & PointsManager(void)
Definition: ManagerAccess.cpp:40
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
Nektar::LibUtilities::P
@ P
Monomial polynomials .
Definition: BasisType.h:62
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
Nektar::LocalRegions::ExpansionVector
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:68
Nektar::MultiRegions::TraceInterpPoints
std::tuple< LibUtilities::PointsKey, LibUtilities::PointsKey, LibUtilities::PointsKey, LibUtilities::PointsKey > TraceInterpPoints
Map holding points distributions required for interpolation of local traces onto global trace in two ...
Definition: LocTraceToTraceMap.h:80
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:49
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1794
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.hpp:539
Vmath::Assmb
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.hpp:577
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::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
std
STL namespace.
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
Generated on Sun Jun 15 2025 07:52:31 for Nektar++ by doxygen 1.9.4