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