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 .get() +
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 .get() +
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 .get() +
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 .get() +
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.get(), m_locInterpTraceToTraceMap[dir].get(),
714 tmp.get());
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.get(), m_locInterpTraceToTraceMap[dir].get(),
793 tmp.get());
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.get(), m_locInterpTraceToTraceMap[dir].get(),
1118 traces.get());
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.get() + cnt, 1,
1174 tmp.get() + 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().get(), tnp, locedges.get() + cnt,
1182 fnp, 0.0, tmp.get() + 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.get() + 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.get(),
1211 m_locInterpTraceToTraceMap[dir].get(), edges.get());
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.get() + cnt, 1,
1263 tmp.get() + 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().get(), tnp0, locfaces.get() + cnt,
1271 fnp0, 0.0, tmp.get() + cnt1, tnp0);
1272 }
1273 break;
1274 case eInterpEndPtDir0:
1275 {
1276 int nfaces = m_interpNtraces[dir][i];
1277 for (int k = 0; k < fnp0; ++k)
1278 {
1279 Vmath::Vcopy(nfaces * fnp1, locfaces.get() + cnt + k,
1280 fnp0, tmp.get() + cnt1 + k, tnp0);
1281 }
1282
1283 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1284 Blas::Dgemv('T', fnp0, tnp1 * nfaces, 1.0, tmp.get() + cnt1,
1285 tnp0, I0.get(), 1, 0.0,
1286 tmp.get() + cnt1 + tnp0 - 1, tnp0);
1287 }
1288 break;
1289 case eInterpDir1:
1290 {
1291 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1292 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1293 {
1294 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1295 locfaces.get() + cnt + j * fnp0 * fnp1,
1296 tnp0, I1->GetPtr().get(), tnp1, 0.0,
1297 tmp.get() + cnt1 + j * tnp0 * tnp1, tnp0);
1298 }
1299 }
1300 break;
1301 case eInterpEndPtDir1:
1302 {
1303 Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1304 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1305 {
1306 // copy all points
1307 Vmath::Vcopy(fnp0 * fnp1,
1308 locfaces.get() + cnt + j * fnp0 * fnp1, 1,
1309 tmp.get() + cnt1 + j * tnp0 * tnp1, 1);
1310
1311 // interpolate end points
1312 for (int k = 0; k < tnp0; ++k)
1313 {
1314 tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1315 Blas::Ddot(fnp1,
1316 locfaces.get() + cnt +
1317 j * fnp0 * fnp1 + k,
1318 fnp0, &I1[0], 1);
1319 }
1320 }
1321 }
1322 break;
1323 case eInterpBothDirs:
1324 {
1325 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1326 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1327 Array<OneD, NekDouble> wsp(m_interpNtraces[dir][i] * fnp0 *
1328 tnp1);
1329
1330 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1331 {
1332 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1333 locfaces.get() + cnt + j * fnp0 * fnp1,
1334 fnp0, I1->GetPtr().get(), tnp1, 0.0,
1335 wsp.get() + j * fnp0 * tnp1, fnp0);
1336 }
1337 Blas::Dgemm('N', 'N', tnp0, tnp1 * m_interpNtraces[dir][i],
1338 fnp0, 1.0, I0->GetPtr().get(), tnp0, wsp.get(),
1339 fnp0, 0.0, tmp.get() + cnt1, tnp0);
1340 }
1341 break;
1342 case eInterpEndPtDir0InterpDir1:
1343 {
1344 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1345 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1346
1347 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1348 {
1349 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1350 locfaces.get() + cnt + j * fnp0 * fnp1,
1351 fnp0, I1->GetPtr().get(), tnp1, 0.0,
1352 tmp.get() + cnt1 + j * tnp0 * tnp1, tnp0);
1353 }
1354
1355 Blas::Dgemv('T', fnp0, tnp1 * m_interpNtraces[dir][i], 1.0,
1356 tmp.get() + cnt1, tnp0, I0.get(), 1, 0.0,
1357 tmp.get() + cnt1 + tnp0 - 1, tnp0);
1358 }
1359 break;
1360 default:
1361 ASSERTL0(false, "Interplation case needs implementing");
1362 break;
1363 }
1364 cnt += nfromfacepts;
1365 cnt1 += m_interpNtraces[dir][i] * tnp0 * tnp1;
1366 }
1367 }
1368
1369 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(), tmp.get(),
1370 m_locInterpTraceToTraceMap[dir].get(), faces.get());
1371}
1372
1373void LocTraceToTraceMap::InterpLocTracesToTraceTranspose(
1374 const int dir, const Array<OneD, const NekDouble> &trace,
1375 Array<OneD, NekDouble> &loctrace)
1376{
1377 switch (m_expdim)
1378 {
1379 case 2:
1380 InterpLocEdgesToTraceTranspose(dir, trace, loctrace);
1381 break;
1382 case 3:
1383 InterpLocFacesToTraceTranspose(dir, trace, loctrace);
1384 break;
1385 default:
1386 NEKERROR(ErrorUtil::efatal, "Not set up");
1387 break;
1388 }
1389}
1390
1391/**
1392 * @brief Transpose of Interp local edges to trace
1393 *
1394 * @param dir Selects forwards (0) or backwards (1) direction.
1395 * @param edges Global trace edge
1396 * @param locedges Local trace edge
1397 */
1398void LocTraceToTraceMap::InterpLocEdgesToTraceTranspose(
1399 const int dir, const Array<OneD, const NekDouble> &edges,
1400 Array<OneD, NekDouble> &locedges)
1401{
1402 ASSERTL1(dir < 2, "option dir out of range, "
1403 " dir=0 is fwd, dir=1 is bwd");
1404
1405 int cnt = 0;
1406 int cnt1 = 0;
1407
1408 // tmp space assuming forward map is of size of trace
1409 Array<OneD, NekDouble> tmp{size_t(m_nTracePts)};
1410 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(), edges.get(),
1411 m_locInterpTraceToTraceMap[dir].get(), tmp.get());
1412
1413 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1414 {
1415 // Check if there are edges to interpolate
1416 if (m_interpNtraces[dir][i])
1417 {
1418 // Get to/from points
1419 LibUtilities::PointsKey fromPointsKey0 =
1420 std::get<0>(m_interpPoints[dir][i]);
1421 LibUtilities::PointsKey toPointsKey0 =
1422 std::get<2>(m_interpPoints[dir][i]);
1423
1424 int fnp = fromPointsKey0.GetNumPoints();
1425 int tnp = toPointsKey0.GetNumPoints();
1426 int nedges = m_interpNtraces[dir][i];
1427
1428 // Do interpolation here if required
1429 switch (m_interpTrace[dir][i])
1430 {
1431 case eNoInterp: // Just copy
1432 {
1433 Vmath::Vcopy(nedges * fnp, tmp.get() + cnt1, 1,
1434 locedges.get() + cnt, 1);
1435 }
1436 break;
1437 case eInterpDir0:
1438 {
1439 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1440 Blas::Dgemm('T', 'N', fnp, nedges, tnp, 1.0,
1441 I0->GetPtr().get(), tnp, tmp.get() + cnt1, tnp,
1442 0.0, locedges.get() + cnt, fnp);
1443 }
1444 break;
1445 case eInterpEndPtDir0:
1446 {
1447 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1448
1449 for (int k = 0; k < nedges; ++k)
1450 {
1451 Vmath::Vcopy(fnp, &tmp[cnt1 + k * tnp], 1,
1452 &locedges[cnt + k * fnp], 1);
1453
1454 Vmath::Svtvp(fnp, tmp[cnt1 + k * tnp + tnp - 1], &I0[0],
1455 1, locedges.get() + cnt + k * fnp, 1,
1456 locedges.get() + cnt + k * fnp, 1);
1457 }
1458 }
1459 break;
1460 default:
1461 NEKERROR(ErrorUtil::efatal,
1462 "Invalid interpolation type for 2D elements");
1463 break;
1464 }
1465
1466 cnt += nedges * fnp;
1467 cnt1 += nedges * tnp;
1468 }
1469 }
1470}
1471
1472/**
1473 * @brief transpose of interp local faces to trace
1474 *
1475 * @param dir Selects forwards (0) or backwards (1) direction.
1476 * @param loctraces Local trace
1477 * @param traces trace .
1478 */
1479void LocTraceToTraceMap::InterpLocFacesToTraceTranspose(
1480 const int dir, const Array<OneD, const NekDouble> &traces,
1481 Array<OneD, NekDouble> &loctraces)
1482{
1483 ASSERTL1(dir < 2, "option dir out of range, "
1484 " dir=0 is fwd, dir=1 is bwd");
1485
1486 int cnt = 0;
1487 int cnt1 = 0;
1488
1489 // tmp space assuming forward map is of size of trace
1490 Array<OneD, NekDouble> tmp{size_t(m_nTracePts)};
1491 // The static cast is necessary because m_locInterpTraceToTraceMap should be
1492 // Array<OneD, size_t> ... or at least the same type as
1493 // m_locInterpTraceToTraceMap.size() ...
1494 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1495 traces.data(), m_locInterpTraceToTraceMap[dir].data(),
1496 tmp.data());
1497
1498 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1499 {
1500 // Check if there are elementboundaries to interpolate
1501 if (m_interpNtraces[dir][i])
1502 {
1503 // Get to/from points
1504 LibUtilities::PointsKey fromPointsKey0 =
1505 std::get<0>(m_interpPoints[dir][i]);
1506 LibUtilities::PointsKey fromPointsKey1 =
1507 std::get<1>(m_interpPoints[dir][i]);
1508 LibUtilities::PointsKey toPointsKey0 =
1509 std::get<2>(m_interpPoints[dir][i]);
1510 LibUtilities::PointsKey toPointsKey1 =
1511 std::get<3>(m_interpPoints[dir][i]);
1512 // Here the f(from) and t(to) are chosen to be consistent with
1513 // InterpLocFacesToTrace
1514 int fnp0 = fromPointsKey0.GetNumPoints();
1515 int fnp1 = fromPointsKey1.GetNumPoints();
1516 int tnp0 = toPointsKey0.GetNumPoints();
1517 int tnp1 = toPointsKey1.GetNumPoints();
1518 int nfromfacepts = m_interpNtraces[dir][i] * fnp0 * fnp1;
1519
1520 // Do transpose interpolation here if required
1521 switch (m_interpTrace[dir][i])
1522 {
1523 case eNoInterp: // Just copy
1524 {
1525 Vmath::Vcopy(nfromfacepts, tmp.get() + cnt1, 1,
1526 loctraces.get() + cnt, 1);
1527 }
1528 break;
1529 case eInterpDir0:
1530 {
1531 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1532 Blas::Dgemm('T', 'N', fnp0, tnp1, tnp0, 1.0,
1533 I0->GetPtr().get(), tnp0, tmp.get() + cnt1,
1534 tnp0, 0.0, loctraces.get() + cnt, fnp0);
1535 }
1536 break;
1537 case eInterpEndPtDir0:
1538 {
1539 int nfaces = m_interpNtraces[dir][i];
1540 for (int k = 0; k < fnp0; ++k)
1541 {
1542 Vmath::Vcopy(nfaces * fnp1, tmp.get() + cnt1 + k, tnp0,
1543 loctraces.get() + cnt + k, fnp0);
1544 }
1545
1546 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1547 for (int k = 0; k < tnp1 * nfaces; k++)
1548 {
1549 Vmath::Svtvp(fnp0, tmp[cnt1 + tnp0 - 1 + k * tnp0],
1550 &I0[0], 1,
1551 loctraces.get() + cnt + k * fnp0, 1,
1552 loctraces.get() + cnt + k * fnp0, 1);
1553 }
1554 }
1555 break;
1556 case eInterpDir1:
1557 {
1558 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1559
1560 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1561 {
1562 Blas::Dgemm('N', 'N', tnp0, fnp1, tnp1, 1.0,
1563 tmp.get() + cnt1 + j * tnp0 * tnp1, tnp0,
1564 I1->GetPtr().get(), tnp1, 0.0,
1565 loctraces.get() + cnt + j * fnp0 * fnp1,
1566 tnp0);
1567 }
1568 }
1569 break;
1570 case eInterpEndPtDir1:
1571 {
1572 Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1573 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1574 {
1575 Vmath::Vcopy(
1576 fnp0 * fnp1, tmp.get() + cnt1 + j * tnp0 * tnp1, 1,
1577 loctraces.get() + cnt + j * fnp0 * fnp1, 1);
1578
1579 for (int k = 0; k < fnp1; k++)
1580 {
1581 Vmath::Svtvp(
1582 fnp0, I1[k],
1583 &tmp[cnt1 + (j + 1) * tnp0 * tnp1 - tnp0], 1,
1584 &loctraces[cnt + j * fnp0 * fnp1 + k * fnp0], 1,
1585 &loctraces[cnt + j * fnp0 * fnp1 + k * fnp0],
1586 1);
1587 }
1588 }
1589 }
1590 break;
1591 case eInterpBothDirs:
1592 {
1593 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1594 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1595
1596 Array<OneD, NekDouble> wsp{
1597 size_t(m_interpNtraces[dir][i] * fnp0 * tnp1)};
1598
1599 Blas::Dgemm('T', 'N', fnp0, tnp1 * m_interpNtraces[dir][i],
1600 tnp0, 1.0, I0->GetPtr().get(), tnp0,
1601 tmp.get() + cnt1, tnp0, 0.0, wsp.get(), fnp0);
1602
1603 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1604 {
1605 Blas::Dgemm('N', 'N', fnp0, fnp1, tnp1, 1.0,
1606 wsp.get() + j * fnp0 * tnp1, fnp0,
1607 I1->GetPtr().get(), tnp1, 0.0,
1608 loctraces.get() + cnt + j * fnp0 * fnp1,
1609 fnp0);
1610 }
1611 }
1612 break;
1613 case eInterpEndPtDir0InterpDir1:
1614 {
1615 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1616 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1617
1618 Array<OneD, NekDouble> wsp{
1619 size_t(m_interpNtraces[dir][i] * fnp0 * tnp1)};
1620
1621 for (int k = 0; k < tnp1 * m_interpNtraces[dir][i]; k++)
1622 {
1623 Vmath::Svtvp(fnp0, tmp[cnt1 + tnp0 - 1 + k * tnp0],
1624 &I0[0], 1, tmp.get() + cnt1 + k * tnp0, 1,
1625 wsp.get() + k * fnp0, 1);
1626 }
1627
1628 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1629 {
1630 Blas::Dgemm('N', 'N', fnp0, fnp1, tnp1, 1.0,
1631 wsp.get() + j * fnp0 * tnp1, fnp0,
1632 I1->GetPtr().get(), tnp1, 0.0,
1633 loctraces.get() + cnt + j * fnp0 * fnp1,
1634 fnp0);
1635 }
1636 }
1637 break;
1638 }
1639 cnt += nfromfacepts;
1640 cnt1 += m_interpNtraces[dir][i] * tnp0 * tnp1;
1641 }
1642 }
1643}
1644
1645void LocTraceToTraceMap::InterpTraceToLocTrace(
1646 const int dir, const Array<OneD, NekDouble> &traces,
1647 Array<OneD, NekDouble> &loctraces)
1648{
1649 switch (m_expdim)
1650 {
1651 case 1: // Essentially do copy
1652 Vmath::Gathr(
1653 static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1654 traces.get(), m_locInterpTraceToTraceMap[dir].get(),
1655 loctraces.get());
1656 break;
1657 case 2:
1658 InterpTraceToLocEdges(dir, traces, loctraces);
1659 break;
1660 case 3:
1661 InterpTraceToLocFaces(dir, traces, loctraces);
1662 break;
1663 default:
1664 ASSERTL0(false, "Not set up");
1665 break;
1666 }
1667}
1668
1669/**
1670 * @brief Interpolate global trace edge to local trace edges point
1671 * distributions where required.
1672 *
1673 * @param dir Selects forwards (0) or backwards (1) direction.
1674 * @param locfaces Local trace edge storage.
1675 * @param faces Global trace edge storage
1676 */
1677void LocTraceToTraceMap::InterpTraceToLocEdges(
1678 const int dir, const Array<OneD, const NekDouble> &edges,
1679 Array<OneD, NekDouble> &locedges)
1680{
1681 ASSERTL1(dir < 2, "option dir out of range, "
1682 " dir=0 is fwd, dir=1 is bwd");
1683
1684 int cnt = 0;
1685 int cnt1 = 0;
1686
1687 Array<OneD, NekDouble> tmp(m_nTracePts);
1688
1689 // unshuffles trace into lcoally orientated format.
1690 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1691 edges.get(), m_locInterpTraceToTraceMap[dir].get(), tmp.get());
1692
1693 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1694 {
1695 // Check if there are edges to interpolate
1696 if (m_interpNtraces[dir][i])
1697 {
1698 // Get to/from points
1699 LibUtilities::PointsKey fromPointsKey0 =
1700 std::get<2>(m_interpPoints[dir][i]);
1701 LibUtilities::PointsKey toPointsKey0 =
1702 std::get<0>(m_interpPoints[dir][i]);
1703
1704 int fnp = fromPointsKey0.GetNumPoints();
1705 int tnp = toPointsKey0.GetNumPoints();
1706 int nedges = m_interpNtraces[dir][i];
1707
1708 // Do interpolation here if required
1709 switch (m_interpTrace[dir][i])
1710 {
1711 case eNoInterp: // Just copy
1712 {
1713 Vmath::Vcopy(nedges * fnp, tmp.get() + cnt, 1,
1714 locedges.get() + cnt1, 1);
1715 }
1716 break;
1717 case eInterpDir0:
1718 {
1719 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
1720 Blas::Dgemm('N', 'N', tnp, nedges, fnp, 1.0,
1721 I0->GetPtr().get(), tnp, tmp.get() + cnt, fnp,
1722 0.0, locedges.get() + cnt1, tnp);
1723 }
1724 break;
1725 case eInterpEndPtDir0:
1726 {
1727 // Just copy points back
1728 for (int k = 0; k < nedges; ++k)
1729 {
1730 // Should be tnp rather than fnp on first argument.
1731 Vmath::Vcopy(tnp, &tmp[cnt + k * fnp], 1,
1732 &locedges[cnt1 + k * tnp], 1);
1733 }
1734 }
1735 break;
1736 default:
1737 ASSERTL0(false,
1738 "Invalid interpolation type for 2D elements");
1739 break;
1740 }
1741
1742 cnt += nedges * fnp;
1743 cnt1 += nedges * tnp;
1744 }
1745 }
1746}
1747
1748/**
1749 * @brief Interpolate global trace edge to local trace edges point
1750 * distributions where required.
1751 *
1752 * @param dir Selects forwards (0) or backwards (1) direction.
1753 * @param locfaces Local trace edge storage.
1754 * @param faces Global trace edge storage
1755 */
1756void LocTraceToTraceMap::InterpTraceToLocFaces(
1757 const int dir, const Array<OneD, const NekDouble> &faces,
1758 Array<OneD, NekDouble> &locfaces)
1759{
1760 ASSERTL1(dir < 2, "option dir out of range, "
1761 " dir=0 is fwd, dir=1 is bwd");
1762
1763 int cnt = 0;
1764 int cnt1 = 0;
1765
1766 Array<OneD, NekDouble> tmp(m_nTracePts);
1767
1768 // unshuffles trace into lcoally orientated format.
1769 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1770 faces.get(), m_locInterpTraceToTraceMap[dir].get(), tmp.get());
1771
1772 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1773 {
1774 // Check if there are faces to interpolate
1775 if (m_interpNtraces[dir][i])
1776 {
1777 // Get to/from points
1778 LibUtilities::PointsKey fromPointsKey0 =
1779 std::get<2>(m_interpPoints[dir][i]);
1780 LibUtilities::PointsKey fromPointsKey1 =
1781 std::get<3>(m_interpPoints[dir][i]);
1782 LibUtilities::PointsKey toPointsKey0 =
1783 std::get<0>(m_interpPoints[dir][i]);
1784 LibUtilities::PointsKey toPointsKey1 =
1785 std::get<1>(m_interpPoints[dir][i]);
1786
1787 int fnp0 = fromPointsKey0.GetNumPoints();
1788 int fnp1 = fromPointsKey1.GetNumPoints();
1789 int tnp0 = toPointsKey0.GetNumPoints();
1790 int tnp1 = toPointsKey1.GetNumPoints();
1791 int nfaces = m_interpNtraces[dir][i];
1792 int nfromfacepts = nfaces * fnp0 * fnp1;
1793
1794 // Do interpolation here if required
1795 switch (m_interpTrace[dir][i])
1796 {
1797 case eNoInterp: // Just copy
1798 {
1799 Vmath::Vcopy(nfromfacepts, tmp.get() + cnt, 1,
1800 locfaces.get() + cnt1, 1);
1801 }
1802 break;
1803 case eInterpDir0:
1804 {
1805 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
1806 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1807 I0->GetPtr().get(), tnp0, tmp.get() + cnt, fnp0,
1808 0.0, locfaces.get() + cnt1, tnp0);
1809 }
1810 break;
1811 case eInterpDir1:
1812 {
1813 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
1814 for (int j = 0; j < nfaces; ++j)
1815 {
1816 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1817 tmp.get() + cnt + j * fnp0 * fnp1, tnp0,
1818 I1->GetPtr().get(), tnp1, 0.0,
1819 locfaces.get() + cnt1 + j * tnp0 * tnp1,
1820 tnp0);
1821 }
1822 }
1823 break;
1824 case eInterpEndPtDir1:
1825 {
1826 for (int j = 0; j < nfaces; ++j)
1827 {
1828 // copy all points missing off top verex in dir 1
1829 Vmath::Vcopy(
1830 tnp0 * tnp1, tmp.get() + cnt + j * fnp0 * fnp1, 1,
1831 locfaces.get() + cnt1 + j * tnp0 * tnp1, 1);
1832 }
1833 }
1834 break;
1835 case eInterpBothDirs:
1836 {
1837 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
1838 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
1839 Array<OneD, NekDouble> wsp(nfaces * fnp0 * tnp1 * fnp0);
1840
1841 for (int j = 0; j < nfaces; ++j)
1842 {
1843 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1844 tmp.get() + cnt + j * fnp0 * fnp1, fnp0,
1845 I1->GetPtr().get(), tnp1, 0.0,
1846 wsp.get() + j * fnp0 * tnp1, fnp0);
1847 }
1848
1849 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1850 I0->GetPtr().get(), tnp0, wsp.get(), fnp0, 0.0,
1851 locfaces.get() + cnt1, tnp0);
1852 }
1853 break;
1854 case eInterpEndPtDir0InterpDir1:
1855 {
1856 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
1857 for (int j = 0; j < nfaces; ++j)
1858 {
1859 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1860 tmp.get() + cnt + j * fnp0 * fnp1, fnp0,
1861 I1->GetPtr().get(), tnp1, 0.0,
1862 locfaces.get() + cnt1 + j * tnp0 * tnp1,
1863 tnp0);
1864 }
1865 }
1866 break;
1867 default:
1868 ASSERTL0(false, "Interpolation case not implemneted (yet)");
1869 break;
1870 }
1871 cnt += nfromfacepts;
1872 cnt1 += m_interpNtraces[dir][i] * tnp0 * tnp1;
1873 }
1874 }
1875}
1876
1877/**
1878 * @brief Add contributions from trace coefficients to the elemental field
1879 * storage.
1880 *
1881 * @param trace Array of global trace coefficients.
1882 * @param field Array containing field coefficients storage.
1883 */
1884void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1885 const Array<OneD, const NekDouble> &trace, Array<OneD, NekDouble> &field)
1886{
1887 int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
1888 for (int i = 0; i < nvals; ++i)
1889 {
1890 field[m_traceCoeffsToElmtMap[0][i]] +=
1891 m_traceCoeffsToElmtSign[0][i] *
1892 trace[m_traceCoeffsToElmtTrace[0][i]];
1893 }
1894}
1895
1896/**
1897 * @brief Add contributions from backwards or forwards oriented trace
1898 * coefficients to the elemental field storage.
1899 *
1900 * @param dir Selects forwards (0) or backwards (1) direction
1901 * @param trace Array of global trace coefficients.
1902 * @param field Array containing field coefficients storage.
1903 */
1904void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1905 const int dir, const Array<OneD, const NekDouble> &trace,
1906 Array<OneD, NekDouble> &field)
1907{
1908 int nvals = m_nTraceCoeffs[dir];
1909 for (int i = 0; i < nvals; ++i)
1910 {
1911 field[m_traceCoeffsToElmtMap[dir][i]] +=
1912 m_traceCoeffsToElmtSign[dir][i] *
1913 trace[m_traceCoeffsToElmtTrace[dir][i]];
1914 }
1915}
1916
1917} // 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:1398
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocTracesToTraceTranspose
void InterpLocTracesToTraceTranspose(const int dir, const Array< OneD, const NekDouble > &traces, Array< OneD, NekDouble > &loctraces)
Definition: LocTraceToTraceMap.cpp:1373
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:1756
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:1677
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:1884
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:1479
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:1645
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:1768
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:298
Generated on Wed Oct 30 2024 16:07:29 for Nektar++ by doxygen 1.9.4