Nektar++
Loading...
Searching...
No Matches
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
53
54/**
55 * @brief Set up trace to trace mapping components.
56 *
57 * @param locExp Expansion list of full dimension problem.
58 * @param trace Expansion list of one dimension lower trace.
59 * @param elmtToTrace Mapping from elemental facets to trace.
60 * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
61 *
62 * @todo Add 1D support
63 */
64LocTraceToTraceMap::LocTraceToTraceMap(
65 const ExpList &locExp, const ExpListSharedPtr &trace,
66 const Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr>>
67 &elmtToTrace,
68 const vector<bool> &LeftAdjacents)
69{
70 const std::shared_ptr<LocalRegions::ExpansionVector> exp = locExp.GetExp();
71 // Assume that all the elements have same dimension
72 m_expdim = (*exp)[0]->GetShapeDimension();
73
74 m_locInterpTraceToTraceMap = Array<OneD, Array<OneD, int>>(2);
75 m_locTraceToElmtTraceMap = Array<OneD, Array<OneD, int>>(2);
76 m_interpTrace = Array<OneD, Array<OneD, InterpLocTraceToTrace>>(2);
77 m_interpTraceI0 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
78 m_interpEndPtI0 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
79 m_interpFromTraceI0 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
80 m_interpPoints = Array<OneD, Array<OneD, TraceInterpPoints>>(2);
81 m_interpNtraces = Array<OneD, Array<OneD, int>>(2);
82
83 if (m_expdim == 3)
84 {
85 m_interpTraceI1 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
86 m_interpFromTraceI1 = Array<OneD, Array<OneD, DNekMatSharedPtr>>(2);
87 m_interpEndPtI1 = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(2);
88 }
89
90 m_traceCoeffsToElmtMap = Array<OneD, Array<OneD, int>>(2);
91 m_traceCoeffsToElmtTrace = Array<OneD, Array<OneD, int>>(2);
92 m_traceCoeffsToElmtSign = Array<OneD, Array<OneD, int>>(2);
93
94 LocalRegions::ExpansionSharedPtr elmt;
95
96 int cnt, n, e, phys_offset;
97
98 int nexp = exp->size();
99 m_nTracePts = trace->GetTotPoints();
100
101 // Count number of traces and points required for maps
102 int nFwdPts = 0;
103 int nBwdPts = 0;
104 int nFwdCoeffs = 0;
105 int nBwdCoeffs = 0;
106 int nLocTraces = 0; // num of total local traces
107 m_nFwdLocTracePts = 0;
108 m_nLocTracePts = 0;
109
110 for (cnt = n = 0; n < nexp; ++n)
111 {
112 elmt = (*exp)[n];
113
114 for (int i = 0; i < elmt->GetNtraces(); ++i, ++cnt)
115 {
116 int nLocPts = elmt->GetTraceNumPoints(i);
117 m_nLocTracePts += nLocPts;
118 nLocTraces++;
119
120 if (LeftAdjacents[cnt])
121 {
122 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
123 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
124 m_nFwdLocTracePts += nLocPts;
125 }
126 else
127 {
128 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
129 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
130 }
131 }
132 }
133
134 m_locTraceToFieldMap = Array<OneD, int>(m_nLocTracePts);
135
136 m_locTraceToElmtTraceMap[0] = Array<OneD, int>(m_nFwdLocTracePts);
137 m_locTraceToElmtTraceMap[1] =
138 Array<OneD, int>(m_nLocTracePts - m_nFwdLocTracePts);
139
140 m_locInterpTraceToTraceMap[0] = Array<OneD, int>(nFwdPts);
141 m_locInterpTraceToTraceMap[1] = Array<OneD, int>(nBwdPts);
142
143 m_nTraceCoeffs[0] = nFwdCoeffs;
144 m_nTraceCoeffs[1] = nBwdCoeffs;
145
146 m_traceCoeffsToElmtMap[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
147 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
148 m_traceCoeffsToElmtTrace[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
149 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
150 m_traceCoeffsToElmtSign[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
151 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
152
153 //------ new entry variables ----------
154 m_locTracePtsEntry = Array<OneD, Array<OneD, int>>(nexp);
155 m_interpTracePtsEntry = Array<OneD, Array<OneD, int>>(nexp);
156 m_interpTraceIndex = Array<OneD, Array<OneD, int>>(nexp);
157 m_traceCoeffsEntry = Array<OneD, Array<OneD, int>>(nexp);
158
159 // size of total loc traces
160 m_locToTracePhysOffset = Array<OneD, int>(nLocTraces);
161
162 // Gather information about trace interpolations
163 map<TraceInterpPoints, vector<pair<int, int>>, cmpop> TraceInterpMap;
164
165 vector<vector<int>> TraceOrder;
166 TraceOrder.resize(nexp);
167 vector<vector<int>> ElmtPhysTraceOffset;
168 ElmtPhysTraceOffset.resize(nexp);
169 int ntrace;
170 int fwdcnt = 0;
171 int bwdcnt = 0;
172 int neoffset = 0;
173 // Generate a map of similar traces with the same
174 // interpolation requirements
175
176 for (cnt = n = 0; n < nexp; ++n)
177 {
178 elmt = (*exp)[n];
179 ntrace = elmt->GetNtraces();
180 TraceOrder[n].resize(ntrace);
181 ElmtPhysTraceOffset[n].resize(ntrace);
182
183 m_locTracePtsEntry[n] = Array<OneD, int>(ntrace);
184 m_interpTracePtsEntry[n] = Array<OneD, int>(ntrace);
185 m_interpTraceIndex[n] = Array<OneD, int>(ntrace);
186 m_traceCoeffsEntry[n] = Array<OneD, int>(ntrace);
187
188 int coeffoffset = locExp.GetCoeff_Offset(n);
189 for (e = 0; e < ntrace; ++e, ++cnt)
190 {
191 LocalRegions::ExpansionSharedPtr elmttrace = elmtToTrace[n][e];
192 StdRegions::Orientation orient = elmt->GetTraceOrient(e);
193
194 LibUtilities::PointsKey fromPointsKey0, fromPointsKey1;
195 LibUtilities::PointsKey toPointsKey0, toPointsKey1;
196 Array<OneD, int> P(2, -1);
197
198 switch (m_expdim)
199 {
200 case 1:
201 {
202 fromPointsKey0 = elmt->GetBasis(0)->GetPointsKey();
203 fromPointsKey1 =
204 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
205 // dummy info since no interpolation is required in this
206 // case.
207 toPointsKey0 =
208 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
209 toPointsKey1 =
210 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
211 }
212 break;
213 case 2:
214 {
215 int dir0 = elmt->GetGeom()->GetDir(e, 0);
216
217 fromPointsKey0 = elmt->GetBasis(dir0)->GetPointsKey();
218 fromPointsKey1 =
219 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
220
221 toPointsKey0 = elmttrace->GetBasis(0)->GetPointsKey();
222 toPointsKey1 =
223 LibUtilities::PointsKey(0, LibUtilities::eNoPointsType);
224
225 P[0] = elmttrace->GetBasisNumModes(0);
226 }
227 break;
228 case 3:
229 {
230 int dir0 = elmt->GetGeom()->GetDir(e, 0);
231 int dir1 = elmt->GetGeom()->GetDir(e, 1);
232
233 fromPointsKey0 = elmt->GetBasis(dir0)->GetPointsKey();
234 fromPointsKey1 = elmt->GetBasis(dir1)->GetPointsKey();
235
236 if (orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
237 {
238 toPointsKey0 = elmttrace->GetBasis(0)->GetPointsKey();
239 toPointsKey1 = elmttrace->GetBasis(1)->GetPointsKey();
240 P[0] = elmttrace->GetBasisNumModes(0);
241 P[1] = elmttrace->GetBasisNumModes(1);
242 }
243 else // transpose points key evaluation
244 {
245 toPointsKey0 = elmttrace->GetBasis(1)->GetPointsKey();
246 toPointsKey1 = elmttrace->GetBasis(0)->GetPointsKey();
247 P[0] = elmttrace->GetBasisNumModes(1);
248 P[1] = elmttrace->GetBasisNumModes(0);
249 }
250 }
251 break;
252 }
253
254 TraceInterpPoints fpoint(fromPointsKey0, fromPointsKey1,
255 toPointsKey0, toPointsKey1);
256
257 pair<int, int> epf(n, e);
258 TraceInterpMap[fpoint].push_back(epf);
259 TraceOrder[n][e] = cnt;
260
261 ElmtPhysTraceOffset[n][e] = neoffset;
262 neoffset += elmt->GetTraceNumPoints(e);
263
264 // Setup for coefficient mapping from trace normal flux
265 // to elements
266 Array<OneD, unsigned int> map;
267 Array<OneD, int> sign;
268 // Test shows we should swap P0 and P1 before calling
269 // GetTraceToElementMap if orientation is transposed.
270 // This is contrary to ReOrientTracePhysMap
271 elmt->GetTraceToElementMap(e, map, sign, orient, P[0], P[1]);
272
273 int order_t = elmttrace->GetNcoeffs();
274 int t_offset = trace->GetCoeff_Offset(elmttrace->GetElmtId());
275
276 double fac = 1.0;
277
278 if (elmt->GetTraceExp(e)->GetRightAdjacentElementExp())
279 {
280 if (elmttrace->GetRightAdjacentElementExp()
281 ->GetGeom()
282 ->GetGlobalID() == elmt->GetGeom()->GetGlobalID())
283 {
284 fac = -1.0;
285 }
286 }
287
288 if (LeftAdjacents[cnt])
289 {
290 m_traceCoeffsEntry[n][e] = fwdcnt;
291 for (int i = 0; i < order_t; ++i)
292 {
293 m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
294 m_traceCoeffsToElmtTrace[0][fwdcnt] = t_offset + i;
295 m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
296 }
297 }
298 else
299 {
300 m_traceCoeffsEntry[n][e] = nFwdCoeffs + bwdcnt;
301 for (int i = 0; i < order_t; ++i)
302 {
303 m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
304 m_traceCoeffsToElmtTrace[1][bwdcnt] = t_offset + i;
305 m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
306 }
307 }
308
309 // Get the start pos of this local trace in the trace phys space
310 m_locToTracePhysOffset[cnt] =
311 trace->GetPhys_Offset(elmttrace->GetElmtId());
312
313 // //-----debug------
314 // if (orient>=9 || P[0]!=P[1])
315 // {
316 // int dir0 = elmt->GetGeom()->GetDir(e, 0);
317 // int dir1 = elmt->GetGeom()->GetDir(e, 1);
318 // std::cout << "n = " << n << " e = " << e
319 // << " orient =" << orient << std::endl;
320 // std::cout << " glotrace nm0, nm1 = " << P[0] << ", "<< P[1]
321 // << std::endl; std::cout << " locTrace nm0, nm1 = " <<
322 // elmt->GetBasisNumModes(dir0)
323 // << ", "<< elmt->GetBasisNumModes(dir1) << std::endl;
324 // }
325 }
326 }
327
328 int nInterpType = TraceInterpMap.size();
329
330 // need to decide on 1D case here !!!!!
331 for (int i = 0; i < 2; ++i)
332 {
333 m_interpTrace[i] = Array<OneD, InterpLocTraceToTrace>(nInterpType);
334 m_interpTraceI0[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
335 m_interpEndPtI0[i] = Array<OneD, Array<OneD, NekDouble>>(nInterpType);
336 m_interpFromTraceI0[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
337 m_interpPoints[i] = Array<OneD, TraceInterpPoints>(nInterpType);
338 m_interpNtraces[i] = Array<OneD, int>(nInterpType, 0);
339 }
340
341 if (m_expdim > 2)
342 {
343 for (int i = 0; i < 2; ++i)
344 {
345 m_interpTraceI1[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
346 m_interpFromTraceI1[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
347 m_interpEndPtI1[i] =
348 Array<OneD, Array<OneD, NekDouble>>(nInterpType);
349 }
350 }
351
352 int ntracepts, ntracepts1;
353 int cnt1 = 0;
354 int cnt2 = 0; // counter for traces
355 int cntFwd = 0;
356 int cntBwd = 0; // counter for loc trace points
357 int cntFwd1 = 0;
358 int cntBwd1 = 0; // counter for trace points
359 int set;
360 Array<OneD, int> traceids;
361 Array<OneD, int> locTraceToTraceMap;
362 cnt = 0;
363
364 for (auto it = TraceInterpMap.begin(); it != TraceInterpMap.end();
365 ++it, ++cnt1)
366 {
367 LibUtilities::PointsKey fromPointsKey0 = std::get<0>(it->first);
368 LibUtilities::PointsKey fromPointsKey1 = std::get<1>(it->first);
369 LibUtilities::PointsKey toPointsKey0 = std::get<2>(it->first);
370 LibUtilities::PointsKey toPointsKey1 = std::get<3>(it->first);
371
372 bool fwdSet = false;
373 bool bwdSet = false;
374
375 for (int f = 0; f < it->second.size(); ++f, ++cnt2)
376 {
377 n = it->second[f].first;
378 e = it->second[f].second;
379
380 StdRegions::StdExpansionSharedPtr elmttrace = elmtToTrace[n][e];
381
382 elmt = (*exp)[n];
383 phys_offset = locExp.GetPhys_Offset(n);
384
385 // Mapping of new edge order to one that loops over elmts
386 // then set up mapping of faces in standard cartesian order
387 elmt->GetTracePhysMap(e, traceids);
388
389 ntracepts = elmt->GetTraceNumPoints(e);
390 ntracepts1 = elmttrace->GetTotPoints();
391
392 StdRegions::Orientation orient = elmt->GetTraceOrient(e);
393
394 elmt->ReOrientTracePhysMap(orient, locTraceToTraceMap,
395 toPointsKey0.GetNumPoints(),
396 toPointsKey1.GetNumPoints(), false);
397
398 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
399
400 if (LeftAdjacents[TraceOrder[n][e]])
401 {
402 m_locTracePtsEntry[n][e] = cntFwd;
403 for (int i = 0; i < ntracepts; ++i)
404 {
405 m_locTraceToFieldMap[cntFwd + i] =
406 phys_offset + traceids[i];
407 }
408
409 for (int i = 0; i < ntracepts; ++i)
410 {
411 m_locTraceToElmtTraceMap[0][cntFwd + i] =
412 ElmtPhysTraceOffset[n][e] + i;
413 }
414
415 m_interpTracePtsEntry[n][e] = cntFwd1;
416 for (int i = 0; i < ntracepts1; ++i)
417 {
418 m_locInterpTraceToTraceMap[0][cntFwd1 + i] =
419 offset + locTraceToTraceMap[i];
420 }
421
422 cntFwd += ntracepts;
423 cntFwd1 += ntracepts1;
424 set = 0;
425 }
426 else
427 {
428 m_locTracePtsEntry[n][e] = m_nFwdLocTracePts + cntBwd;
429 for (int i = 0; i < ntracepts; ++i)
430 {
431 m_locTraceToFieldMap[m_nFwdLocTracePts + cntBwd + i] =
432 phys_offset + traceids[i];
433 }
434
435 for (int i = 0; i < ntracepts; ++i)
436 {
437 m_locTraceToElmtTraceMap[1][cntBwd + i] =
438 ElmtPhysTraceOffset[n][e] + i;
439 }
440
441 // it doesn't matter if it is fwd or bwd
442 m_interpTracePtsEntry[n][e] = cntBwd1;
443 for (int i = 0; i < ntracepts1; ++i)
444 {
445 m_locInterpTraceToTraceMap[1][cntBwd1 + i] =
446 offset + locTraceToTraceMap[i];
447 }
448
449 cntBwd += ntracepts;
450 cntBwd1 += ntracepts1;
451 set = 1;
452 }
453
454 m_interpNtraces[set][cnt1] += 1;
455 m_interpTraceIndex[n][e] = cnt1;
456 // For each unique interpolation type, set the interp matrix only
457 // once either for fwdSet or bwdSet
458 if ((fwdSet == false && set == 0) || (bwdSet == false && set == 1))
459 {
460 m_interpPoints[set][cnt1] = it->first;
461
462 switch (m_expdim)
463 {
464 case 1:
465 {
466 // Always no interplation in this case
467 m_interpTrace[set][cnt1] = eNoInterp;
468 }
469 break;
470 case 2:
471 {
472 if (fromPointsKey0 == toPointsKey0)
473 {
474 m_interpTrace[set][cnt1] = eNoInterp;
475 }
476 else
477 {
478 m_interpTrace[set][cnt1] = eInterpDir0;
479 m_interpTraceI0[set][cnt1] =
480 LibUtilities::PointsManager()[fromPointsKey0]
481 ->GetI(toPointsKey0);
482 m_interpFromTraceI0[set][cnt1] =
483 LibUtilities::PointsManager()[toPointsKey0]
484 ->GetI(fromPointsKey0);
485 // Check to see if we can
486 // just interpolate endpoint
487 if (fromPointsKey0.GetNumPoints() + 1 ==
488 toPointsKey0.GetNumPoints())
489 {
490 if (((fromPointsKey0.GetPointsType() ==
491 LibUtilities::eGaussRadauMAlpha1Beta0) &&
492 (toPointsKey0.GetPointsType() ==
493 LibUtilities::eGaussLobattoLegendre)) ||
494 ((fromPointsKey0.GetPointsType() ==
495 LibUtilities::eGaussLegendreWithM) &&
496 (toPointsKey0.GetPointsType() ==
497 LibUtilities::eGaussLegendreWithMP)))
498 {
499 m_interpTrace[set][cnt1] = eInterpEndPtDir0;
500
501 int fnp0 = fromPointsKey0.GetNumPoints();
502 int tnp0 = toPointsKey0.GetNumPoints();
503
504 m_interpEndPtI0[set][cnt1] =
505 Array<OneD, NekDouble>(fnp0);
506
507 Vmath::Vcopy(fnp0,
508 m_interpTraceI0[set][cnt1]
509 ->GetPtr()
510 .data() +
511 tnp0 - 1,
512 tnp0,
513 &m_interpEndPtI0[set][cnt1][0],
514 1);
515 }
516 }
517 }
518 }
519 break;
520 case 3:
521 {
522 if (fromPointsKey0 == toPointsKey0)
523 {
524 if (fromPointsKey1 == toPointsKey1)
525 {
526 m_interpTrace[set][cnt1] = eNoInterp;
527 }
528 else
529 {
530 m_interpTrace[set][cnt1] = eInterpDir1;
531 m_interpTraceI1[set][cnt1] =
532 LibUtilities::PointsManager()
533 [fromPointsKey1]
534 ->GetI(toPointsKey1);
535 m_interpFromTraceI1[set][cnt1] =
536 LibUtilities::PointsManager()[toPointsKey1]
537 ->GetI(fromPointsKey1);
538
539 // Check to see if we can just
540 // interpolate endpoint
541 if (fromPointsKey1.GetNumPoints() + 1 ==
542 toPointsKey1.GetNumPoints())
543 {
544 if (((fromPointsKey1.GetPointsType() ==
545 LibUtilities::
546 eGaussRadauMAlpha1Beta0) &&
547 (toPointsKey1.GetPointsType() ==
548 LibUtilities::
549 eGaussLobattoLegendre)) ||
550 ((fromPointsKey1.GetPointsType() ==
551 LibUtilities::eGaussLegendreWithM) &&
552 (toPointsKey1.GetPointsType() ==
553 LibUtilities::eGaussLegendreWithMP)))
554 {
555 m_interpTrace[set][cnt1] =
556 eInterpEndPtDir1;
557 int fnp1 =
558 fromPointsKey1.GetNumPoints();
559 int tnp1 = toPointsKey1.GetNumPoints();
560 m_interpEndPtI1[set][cnt1] =
561 Array<OneD, NekDouble>(fnp1);
562 Vmath::Vcopy(
563 fnp1,
564 m_interpTraceI1[set][cnt1]
565 ->GetPtr()
566 .data() +
567 tnp1 - 1,
568 tnp1,
569 &m_interpEndPtI1[set][cnt1][0], 1);
570 }
571 }
572 }
573 }
574 else
575 {
576 if (fromPointsKey1 == toPointsKey1)
577 {
578 m_interpTrace[set][cnt1] = eInterpDir0;
579 m_interpTraceI0[set][cnt1] =
580 LibUtilities::PointsManager()
581 [fromPointsKey0]
582 ->GetI(toPointsKey0);
583 m_interpFromTraceI0[set][cnt1] =
584 LibUtilities::PointsManager()[toPointsKey0]
585 ->GetI(fromPointsKey0);
586
587 // Check to see if we can just
588 // interpolate endpoint
589 if (fromPointsKey0.GetNumPoints() + 1 ==
590 toPointsKey0.GetNumPoints())
591 {
592 if (((fromPointsKey0.GetPointsType() ==
593 LibUtilities::
594 eGaussRadauMAlpha1Beta0) &&
595 (toPointsKey0.GetPointsType() ==
596 LibUtilities::
597 eGaussLobattoLegendre)) ||
598 ((fromPointsKey0.GetPointsType() ==
599 LibUtilities::eGaussLegendreWithM) &&
600 (toPointsKey0.GetPointsType() ==
601 LibUtilities::eGaussLegendreWithMP)))
602 {
603 m_interpTrace[set][cnt1] =
604 eInterpEndPtDir0;
605 int fnp0 =
606 fromPointsKey0.GetNumPoints();
607 int tnp0 = toPointsKey0.GetNumPoints();
608 m_interpEndPtI0[set][cnt1] =
609 Array<OneD, NekDouble>(fnp0);
610 Vmath::Vcopy(
611 fnp0,
612 m_interpTraceI0[set][cnt1]
613 ->GetPtr()
614 .data() +
615 tnp0 - 1,
616 tnp0,
617 &m_interpEndPtI0[set][cnt1][0], 1);
618 }
619 }
620 }
621 else
622 {
623 m_interpTrace[set][cnt1] = eInterpBothDirs;
624 m_interpTraceI0[set][cnt1] =
625 LibUtilities::PointsManager()
626 [fromPointsKey0]
627 ->GetI(toPointsKey0);
628 m_interpFromTraceI0[set][cnt1] =
629 LibUtilities::PointsManager()[toPointsKey0]
630 ->GetI(fromPointsKey0);
631 m_interpTraceI1[set][cnt1] =
632 LibUtilities::PointsManager()
633 [fromPointsKey1]
634 ->GetI(toPointsKey1);
635 m_interpFromTraceI1[set][cnt1] =
636 LibUtilities::PointsManager()[toPointsKey1]
637 ->GetI(fromPointsKey1);
638
639 // check to see if we can just
640 // interpolate endpoint
641 if (fromPointsKey0.GetNumPoints() + 1 ==
642 toPointsKey0.GetNumPoints())
643 {
644 if (((fromPointsKey0.GetPointsType() ==
645 LibUtilities::
646 eGaussRadauMAlpha1Beta0) &&
647 (toPointsKey0.GetPointsType() ==
648 LibUtilities::
649 eGaussLobattoLegendre)) ||
650 ((fromPointsKey0.GetPointsType() ==
651 LibUtilities::eGaussLegendreWithM) &&
652 (toPointsKey0.GetPointsType() ==
653 LibUtilities::eGaussLegendreWithMP)))
654 {
655 m_interpTrace[set][cnt1] =
656 eInterpEndPtDir0InterpDir1;
657 int fnp0 =
658 fromPointsKey0.GetNumPoints();
659 int tnp0 = toPointsKey0.GetNumPoints();
660 m_interpEndPtI0[set][cnt1] =
661 Array<OneD, NekDouble>(fnp0);
662 Vmath::Vcopy(
663 fnp0,
664 m_interpTraceI0[set][cnt1]
665 ->GetPtr()
666 .data() +
667 tnp0 - 1,
668 tnp0,
669 &m_interpEndPtI0[set][cnt1][0], 1);
670 }
671 }
672 }
673 }
674 }
675 }
676
677 if (set == 0)
678 {
679 fwdSet = true;
680 }
681 else
682 {
683 bwdSet = true;
684 }
685 }
686 }
687 }
688
689 TraceLocToElmtLocCoeffMap(locExp, trace);
690 FindElmtNeighbors(locExp, trace);
691
692 //---- Below is to construct traceFieldMap and traceInterp essential -----
693
694 auto collections = locExp.GetCollections();
695
696 m_collExpOffset = Array<OneD, int>(collections.size(), 0);
697
698 for (int cid = 1; cid < collections.size(); ++cid)
699 {
700 m_collExpOffset[cid] = m_collExpOffset[cid - 1] +
701 collections[cid - 1].GetExpVector().size();
702 }
703
704 m_traceFieldMap = Array<OneD, TraceFieldMapEssential>(collections.size());
705 m_traceInterp = Array<OneD, TraceInterpEssential>(collections.size());
706
707 // map<{shape, orient, dir, typid}, orientationid>
708 std::map<std::tuple<int, int, int, int>, int> map_locTraceToTrace;
709
710 // temporary
711 std::vector<Array<OneD, int>> orientationMaps;
712
713 int cntLocTrace = 0;
714 int elmtId = 0;
715 for (int cid = 0; cid < collections.size(); ++cid)
716 {
717 auto pCollExp = collections[cid].GetExpVector();
718 auto pGeomData = collections[cid].GetGeomSharedPtr();
719
720 auto exp =
721 std::dynamic_pointer_cast<LocalRegions::Expansion>(pCollExp[0]);
722
723 // directly assign the same name variables
724 m_traceInterp[cid].m_interpTrace = m_interpTrace;
725 m_traceInterp[cid].m_interpTraceI0 = m_interpTraceI0;
726 m_traceInterp[cid].m_interpTraceI1 = m_interpTraceI1;
727 m_traceInterp[cid].m_interpFromTraceI0 = m_interpFromTraceI0;
728 m_traceInterp[cid].m_interpFromTraceI1 = m_interpFromTraceI1;
729 m_traceInterp[cid].m_interpPoints = m_interpPoints;
730 m_traceInterp[cid].m_interpEndPtI0 = m_interpEndPtI0;
731 m_traceInterp[cid].m_interpEndPtI1 = m_interpEndPtI1;
732 m_traceInterp[cid].m_interpNtraces = m_interpNtraces;
733
734 m_traceInterp[cid].m_maxTraceSize = 1;
735 for (int dir = 0; dir < 2; ++dir)
736 {
737 for (int typid = 0; typid < m_interpPoints[dir].size(); ++typid)
738 {
739 if (m_expdim == 3)
740 {
741 int fnp0 =
742 std::get<0>(m_interpPoints[dir][typid]).GetNumPoints();
743 int fnp1 =
744 std::get<1>(m_interpPoints[dir][typid]).GetNumPoints();
745 int tnp0 =
746 std::get<2>(m_interpPoints[dir][typid]).GetNumPoints();
747 int tnp1 =
748 std::get<3>(m_interpPoints[dir][typid]).GetNumPoints();
749
750 int tmp = std::max(fnp0 * fnp1, tnp0 * tnp1);
751 m_traceInterp[cid].m_maxTraceSize =
752 std::max(tmp, m_traceInterp[cid].m_maxTraceSize);
753 }
754 else
755 {
756 int fnp0 =
757 std::get<0>(m_interpPoints[dir][typid]).GetNumPoints();
758 int tnp0 =
759 std::get<2>(m_interpPoints[dir][typid]).GetNumPoints();
760 int tmp = std::max(fnp0, tnp0);
761 m_traceInterp[cid].m_maxTraceSize =
762 std::max(tmp, m_traceInterp[cid].m_maxTraceSize);
763 }
764 }
765 }
766
767 // reshape m_interpTraceIndex of this collection from (n,e) to (e,n)
768 // and stored in m_traceInterp
769 m_traceInterp[cid].m_interpTraceIndex =
770 Array<OneD, Array<OneD, int>>(exp->GetNtraces());
771 for (int e = 0; e < exp->GetNtraces(); e++)
772 {
773 m_traceInterp[cid].m_interpTraceIndex[e] =
774 Array<OneD, int>(pCollExp.size(), 0);
775 }
776
777 // Get trace to elmt map from first exp and stored
778 m_traceFieldMap[cid].m_locTracePhysToElmtMaps =
779 Array<OneD, Array<OneD, int>>(exp->GetNtraces());
780
781 for (int e = 0; e < exp->GetNtraces(); e++)
782 {
783 exp->GetTracePhysMap(
784 e, m_traceFieldMap[cid].m_locTracePhysToElmtMaps[e]);
785 }
786
787 // Allocate memory for orientations
788 m_traceFieldMap[cid].m_orientationIds =
789 Array<OneD, Array<OneD, int>>(exp->GetNtraces());
790 for (int e = 0; e < exp->GetNtraces(); e++)
791 {
792 m_traceFieldMap[cid].m_orientationIds[e] =
793 Array<OneD, int>(pCollExp.size(), 0);
794 }
795
796 // allocate memory for locToTraceId
797 m_traceFieldMap[cid].m_locToTraceId =
798 Array<OneD, Array<OneD, int>>(exp->GetNtraces());
799 for (int e = 0; e < exp->GetNtraces(); e++)
800 {
801 m_traceFieldMap[cid].m_locToTraceId[e] =
802 Array<OneD, int>(pCollExp.size(), 0);
803 }
804 // clear the maps : so each collection has its own maps
805 // However, even if we disable this, things won't go wrong
806 // Just have bigger arrays storing redundant data
807 orientationMaps.clear();
808 map_locTraceToTrace.clear();
809
810 int cntCollLocTrace = 0;
811 for (int n = 0; n < pCollExp.size(); n++, elmtId++)
812 {
813 auto locExp =
814 std::dynamic_pointer_cast<LocalRegions::Expansion>(pCollExp[n]);
815 for (int e = 0; e < locExp->GetNtraces(); e++, cntCollLocTrace++)
816 {
817 int dir = 0;
818 if (!LeftAdjacents[cntLocTrace + cntCollLocTrace])
819 {
820 dir = 1;
821 }
822 int typid = m_interpTraceIndex[elmtId][e];
823 int orient = static_cast<int>(locExp->GetTraceOrient(e));
824 int shape = 0;
825 switch (locExp->GetShapeDimension())
826 {
827 case 1:
828 shape = static_cast<int>(LibUtilities::ePoint);
829 break;
830 case 2:
831 shape = static_cast<int>(
832 locExp->GetGeom()->GetEdge(e)->GetShapeType());
833 break;
834 case 3:
835 shape = static_cast<int>(
836 locExp->GetGeom()->GetFace(e)->GetShapeType());
837 break;
838 }
839 auto thisKey =
840 std::tuple<int, int, int, int>{shape, orient, dir, typid};
841 if (map_locTraceToTrace.find(thisKey) ==
842 map_locTraceToTrace.end())
843 {
844 // insert new item
845 Array<OneD, int> thisMapArray;
846 auto tnp0 =
847 std::get<2>(m_interpPoints[dir][typid]).GetNumPoints();
848 auto tnp1 =
849 std::get<3>(m_interpPoints[dir][typid]).GetNumPoints();
850 exp->ReOrientTracePhysMap(locExp->GetTraceOrient(e),
851 thisMapArray, tnp0, tnp1, false);
852 orientationMaps.push_back(thisMapArray);
853 m_traceFieldMap[cid].m_orientationIds[e][n] =
854 orientationMaps.size() - 1;
855 map_locTraceToTrace[thisKey] = orientationMaps.size() - 1;
856 }
857 else // already exist
858 {
859 m_traceFieldMap[cid].m_orientationIds[e][n] =
860 map_locTraceToTrace[thisKey];
861 }
862
863 std::array<LibUtilities::BasisKey, 3> thisExpKeys{
864 LibUtilities::NullBasisKey, LibUtilities::NullBasisKey,
865 LibUtilities::NullBasisKey};
866
867 m_traceInterp[cid].m_interpTraceIndex[e][n] =
868 m_interpTraceIndex[elmtId][e];
869 m_traceFieldMap[cid].m_locToTraceId[e][n] =
870 elmtToTrace[elmtId][e]->GetElmtId();
871 }
872 }
873
874 // copy orientationMaps to m_traceFieldMap[cid] because orientationMaps
875 // are std::vector
876 m_traceFieldMap[cid].m_orientationMaps =
877 Array<OneD, Array<OneD, int>>(orientationMaps.size());
878 for (int i = 0; i < orientationMaps.size(); ++i)
879 {
880 m_traceFieldMap[cid].m_orientationMaps[i] = orientationMaps[i];
881 }
882
883 // reshape LeftAdjacents from (n,e) to (e,n) and copy to
884 // m_traceFieldMap[cid].m_isLocTraceLeftAdjacent
885 m_traceFieldMap[cid].m_isLocTraceLeftAdjacent =
886 Array<OneD, Array<OneD, bool>>(exp->GetNtraces());
887 for (int e = 0; e < exp->GetNtraces(); e++)
888 {
889 m_traceFieldMap[cid].m_isLocTraceLeftAdjacent[e] =
890 Array<OneD, bool>(pCollExp.size(), false);
891 for (int n = 0; n < pCollExp.size(); n++)
892 {
893 m_traceFieldMap[cid].m_isLocTraceLeftAdjacent[e][n] =
894 LeftAdjacents[cntLocTrace + n * exp->GetNtraces() + e];
895 }
896 }
897
898 // copy the m_locToTracePhysOffset with offset to a new array
899 // reshape from (n,e) to (e,n) and assgin to m_traceFieldMap[cid]
900 m_traceFieldMap[cid].m_locToTracePhysOffset =
901 Array<OneD, Array<OneD, int>>(exp->GetNtraces());
902 for (int e = 0; e < exp->GetNtraces(); e++)
903 {
904 m_traceFieldMap[cid].m_locToTracePhysOffset[e] =
905 Array<OneD, int>(pCollExp.size(), 0);
906 for (int n = 0; n < pCollExp.size(); n++)
907 {
908 m_traceFieldMap[cid].m_locToTracePhysOffset[e][n] =
909 m_locToTracePhysOffset[cntLocTrace + n * exp->GetNtraces() +
910 e];
911 }
912 }
913
914 // update counters
915 cntLocTrace += cntCollLocTrace;
916 }
917}
918
919const TraceInterpEssential &LocTraceToTraceMap::GetTraceInterpEssential(
920 const int cid)
921{
922 return m_traceInterp[cid];
923}
924
930
931void LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap(
932 const ExpListSharedPtr &tracelist, const int ndim)
933{
934 switch (ndim)
935 {
936 case 2:
937 CalcLocTracePhysToTraceIDMap_2D(tracelist);
938 break;
939 case 3:
940 CalcLocTracePhysToTraceIDMap_3D(tracelist);
941 break;
942 default:
943 NEKERROR(ErrorUtil::efatal,
944 "CalcLocTracePhysToTraceIDMap not coded");
945 }
946}
947
948void LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_2D(
949 const ExpListSharedPtr &tracelist)
950{
951 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
952 tracelist->GetExp();
953 int ntotTrace = (*traceExp).size();
954 int ntPnts, noffset;
955
956 m_LocTracephysToTraceIDMap = Array<OneD, Array<OneD, int>>(2);
957 m_LocTracephysToTraceIDMap[0] = Array<OneD, int>(m_nFwdLocTracePts, -1);
958 m_LocTracephysToTraceIDMap[1] =
959 Array<OneD, int>(m_nLocTracePts - m_nFwdLocTracePts, -1);
960
961 Array<OneD, NekDouble> tracePnts(m_nTracePts, 0.0);
962 for (int nt = 0; nt < ntotTrace; nt++)
963 {
964 ntPnts = tracelist->GetTotPoints(nt);
965 noffset = tracelist->GetPhys_Offset(nt);
966 for (int i = 0; i < ntPnts; i++)
967 {
968 tracePnts[noffset + i] = NekDouble(nt);
969 }
970 }
971
972 Array<OneD, Array<OneD, NekDouble>> loctracePntsLR(2);
973 loctracePntsLR[0] = Array<OneD, NekDouble>(m_nFwdLocTracePts, 0.0);
974 loctracePntsLR[1] =
975 Array<OneD, NekDouble>(m_nLocTracePts - m_nFwdLocTracePts, 0.0);
976
977 for (int dir = 0; dir < 2; dir++)
978 {
979 int cnt = 0;
980 int cnt1 = 0;
981
982 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
983 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(),
984 tracePnts.data(), m_locInterpTraceToTraceMap[dir].data(),
985 tmp.data());
986
987 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
988 {
989 if (m_interpNtraces[dir][i])
990 {
991 LibUtilities::PointsKey fromPointsKey0 =
992 std::get<0>(m_interpPoints[dir][i]);
993 LibUtilities::PointsKey toPointsKey0 =
994 std::get<2>(m_interpPoints[dir][i]);
995
996 int fnp = fromPointsKey0.GetNumPoints();
997 int tnp = toPointsKey0.GetNumPoints();
998 int nedges = m_interpNtraces[dir][i];
999
1000 for (int ne = 0; ne < nedges; ne++)
1001 {
1002 Vmath::Fill(fnp, tmp[cnt1], &loctracePntsLR[dir][cnt], 1);
1003 cnt += fnp;
1004 cnt1 += tnp;
1005 }
1006 }
1007 }
1008 }
1009
1010 NekDouble error = 0.0;
1011 for (int nlr = 0; nlr < 2; nlr++)
1012 {
1013 for (int i = 0; i < loctracePntsLR[nlr].size(); i++)
1014 {
1015 m_LocTracephysToTraceIDMap[nlr][i] =
1016 std::round(loctracePntsLR[nlr][i]);
1017 error += abs(loctracePntsLR[nlr][i] -
1018 NekDouble(m_LocTracephysToTraceIDMap[nlr][i]));
1019 }
1020 }
1021 error = error / NekDouble(m_nLocTracePts);
1022 ASSERTL0(error < NekConstants::kNekZeroTol,
1023 "m_LocTracephysToTraceIDMap may not be integer !!");
1024}
1025
1026void LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_3D(
1027 const ExpListSharedPtr &tracelist)
1028{
1029 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
1030 tracelist->GetExp();
1031 int ntotTrace = (*traceExp).size();
1032 int ntPnts, noffset;
1033
1034 m_LocTracephysToTraceIDMap = Array<OneD, Array<OneD, int>>(2);
1035 m_LocTracephysToTraceIDMap[0] = Array<OneD, int>(m_nFwdLocTracePts, -1);
1036 m_LocTracephysToTraceIDMap[1] =
1037 Array<OneD, int>(m_nLocTracePts - m_nFwdLocTracePts, -1);
1038
1039 Array<OneD, NekDouble> tracePnts(m_nTracePts, 0.0);
1040 for (int nt = 0; nt < ntotTrace; nt++)
1041 {
1042 ntPnts = tracelist->GetTotPoints(nt);
1043 noffset = tracelist->GetPhys_Offset(nt);
1044 for (int i = 0; i < ntPnts; i++)
1045 {
1046 tracePnts[noffset + i] = NekDouble(nt);
1047 }
1048 }
1049
1050 Array<OneD, Array<OneD, NekDouble>> loctracePntsLR(2);
1051 loctracePntsLR[0] = Array<OneD, NekDouble>(m_nFwdLocTracePts, 0.0);
1052 loctracePntsLR[1] =
1053 Array<OneD, NekDouble>(m_nLocTracePts - m_nFwdLocTracePts, 0.0);
1054
1055 for (int dir = 0; dir < 2; dir++)
1056 {
1057 int cnt = 0;
1058 int cnt1 = 0;
1059
1060 // tmp space assuming forward map is of size of trace
1061 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
1062 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(),
1063 tracePnts.data(), m_locInterpTraceToTraceMap[dir].data(),
1064 tmp.data());
1065
1066 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1067 {
1068 if (m_interpNtraces[dir][i])
1069 {
1070 LibUtilities::PointsKey fromPointsKey0 =
1071 std::get<0>(m_interpPoints[dir][i]);
1072 LibUtilities::PointsKey fromPointsKey1 =
1073 std::get<1>(m_interpPoints[dir][i]);
1074 LibUtilities::PointsKey toPointsKey0 =
1075 std::get<2>(m_interpPoints[dir][i]);
1076 LibUtilities::PointsKey toPointsKey1 =
1077 std::get<3>(m_interpPoints[dir][i]);
1078
1079 int fnp0 = fromPointsKey0.GetNumPoints();
1080 int fnp1 = fromPointsKey1.GetNumPoints();
1081 int tnp0 = toPointsKey0.GetNumPoints();
1082 int tnp1 = toPointsKey1.GetNumPoints();
1083
1084 int nfttl = fnp0 * fnp1;
1085
1086 for (int ne = 0; ne < m_interpNtraces[dir][i]; ne++)
1087 {
1088 Vmath::Fill(nfttl, tmp[cnt1], &loctracePntsLR[dir][cnt], 1);
1089 cnt += nfttl;
1090 cnt1 += tnp0 * tnp1;
1091 }
1092 }
1093 }
1094 }
1095
1096 NekDouble error = 0.0;
1097 for (int nlr = 0; nlr < 2; nlr++)
1098 {
1099 for (int i = 0; i < loctracePntsLR[nlr].size(); i++)
1100 {
1101 m_LocTracephysToTraceIDMap[nlr][i] =
1102 std::round(loctracePntsLR[nlr][i]);
1103 error += abs(loctracePntsLR[nlr][i] -
1104 NekDouble(m_LocTracephysToTraceIDMap[nlr][i]));
1105 }
1106 }
1107 error = error / NekDouble(m_nLocTracePts);
1108 ASSERTL0(error < NekConstants::kNekZeroTol,
1109 "m_LocTracephysToTraceIDMap may not be integer !!");
1110}
1111
1112/**
1113 * @brief Set up maps between coefficients on trace and in cells.
1114 *
1115 * @param locExp Expansion list in elements
1116 * @param trace Expansion list on traces.
1117 */
1118void LocTraceToTraceMap::TraceLocToElmtLocCoeffMap(
1119 const ExpList &locExp, const ExpListSharedPtr &trace)
1120{
1121 const std::shared_ptr<LocalRegions::ExpansionVector> exptrac =
1122 trace->GetExp();
1123 size_t ntrace = exptrac->size();
1124
1125 Array<OneD, Array<OneD, int>> LRAdjExpid{2};
1126 Array<OneD, Array<OneD, bool>> LRAdjflag{2};
1127
1128 TensorOfArray3D<int> elmtLRMap{2};
1129 TensorOfArray3D<int> elmtLRSign{2};
1130
1131 for (int lr = 0; lr < 2; ++lr)
1132 {
1133 LRAdjExpid[lr] = Array<OneD, int>{ntrace, 0};
1134 LRAdjflag[lr] = Array<OneD, bool>{ntrace, false};
1135 elmtLRMap[lr] = Array<OneD, Array<OneD, int>>{ntrace};
1136 elmtLRSign[lr] = Array<OneD, Array<OneD, int>>{ntrace};
1137 for (int i = 0; i < ntrace; ++i)
1138 {
1139 size_t ncoeff = trace->GetNcoeffs(i);
1140 elmtLRMap[lr][i] = Array<OneD, int>{ncoeff, 0};
1141 elmtLRSign[lr][i] = Array<OneD, int>{ncoeff, 0};
1142 }
1143 }
1144
1145 const Array<OneD, const pair<int, int>> field_coeffToElmt =
1146 locExp.GetCoeffsToElmt();
1147 const Array<OneD, const pair<int, int>> trace_coeffToElmt =
1148 trace->GetCoeffsToElmt();
1149
1150 for (int lr = 0; lr < 2; ++lr)
1151 {
1152 int ntotcoeffs = m_nTraceCoeffs[lr];
1153 for (int i = 0; i < ntotcoeffs; ++i)
1154 {
1155 int ncoeffField = m_traceCoeffsToElmtMap[lr][i];
1156 int ncoeffTrace = m_traceCoeffsToElmtTrace[lr][i];
1157 int sign = m_traceCoeffsToElmtSign[lr][i];
1158
1159 int ntraceelmt = trace_coeffToElmt[ncoeffTrace].first;
1160 int ntracelocN = trace_coeffToElmt[ncoeffTrace].second;
1161
1162 int nfieldelmt = field_coeffToElmt[ncoeffField].first;
1163 int nfieldlocN = field_coeffToElmt[ncoeffField].second;
1164
1165 LRAdjflag[lr][ntraceelmt] = true;
1166 LRAdjExpid[lr][ntraceelmt] = nfieldelmt;
1167
1168 elmtLRMap[lr][ntraceelmt][ntracelocN] = nfieldlocN;
1169 elmtLRSign[lr][ntraceelmt][ntracelocN] = sign;
1170 }
1171 }
1172 m_leftRightAdjacentExpId = LRAdjExpid;
1173 m_leftRightAdjacentExpFlag = LRAdjflag;
1174 m_traceCoeffToLeftRightExpCoeffMap = elmtLRMap;
1175 m_traceCoeffToLeftRightExpCoeffSign = elmtLRSign;
1176}
1177
1178void LocTraceToTraceMap::FindElmtNeighbors(const ExpList &locExp,
1179 const ExpListSharedPtr &trace)
1180{
1181 const std::shared_ptr<LocalRegions::ExpansionVector> exptrac =
1182 trace->GetExp();
1183 int ntrace = exptrac->size();
1184
1185 const std::shared_ptr<LocalRegions::ExpansionVector> exp = locExp.GetExp();
1186 int nexp = exp->size();
1187
1188 Array<OneD, Array<OneD, int>> LRAdjExpid(2);
1189 Array<OneD, Array<OneD, bool>> LRAdjflag(2);
1190 LRAdjExpid = m_leftRightAdjacentExpId;
1191 LRAdjflag = m_leftRightAdjacentExpFlag;
1192
1193 std::set<std::pair<int, int>> neighborSet;
1194 int ntmp0, ntmp1;
1195 for (int nt = 0; nt < ntrace; nt++)
1196 {
1197 if (LRAdjflag[0][nt] && LRAdjflag[1][nt])
1198 {
1199 ntmp0 = LRAdjExpid[0][nt];
1200 ntmp1 = LRAdjExpid[1][nt];
1201
1202 ASSERTL0(ntmp0 != ntmp1,
1203 " ntmp0==ntmp1, trace inside a element?? ");
1204
1205 std::set<std::pair<int, int>>::iterator it = neighborSet.begin();
1206 neighborSet.insert(it, std::make_pair(ntmp0, ntmp1));
1207 neighborSet.insert(it, std::make_pair(ntmp1, ntmp0));
1208 }
1209 }
1210
1211 Array<OneD, int> ElemIndex(nexp, 0);
1212 for (std::set<std::pair<int, int>>::iterator it = neighborSet.begin();
1213 it != neighborSet.end(); ++it)
1214 {
1215 int ncurrent = it->first;
1216 ElemIndex[ncurrent]++;
1217 }
1218
1219 Array<OneD, Array<OneD, int>> ElemNeighbsId(nexp);
1220 Array<OneD, Array<OneD, int>> tmpId(nexp);
1221 Array<OneD, int> ElemNeighbsNumb(nexp, -1);
1222 Vmath::Vcopy(nexp, ElemIndex, 1, ElemNeighbsNumb, 1);
1223 for (int ne = 0; ne < nexp; ne++)
1224 {
1225 int neighb = ElemNeighbsNumb[ne];
1226 ElemNeighbsId[ne] = Array<OneD, int>(neighb, -1);
1227 tmpId[ne] = Array<OneD, int>(neighb, -1);
1228 }
1229
1230 for (int ne = 0; ne < nexp; ne++)
1231 {
1232 ElemIndex[ne] = 0;
1233 }
1234 for (std::set<std::pair<int, int>>::iterator it = neighborSet.begin();
1235 it != neighborSet.end(); ++it)
1236 {
1237 int ncurrent = it->first;
1238 int neighbor = it->second;
1239 ElemNeighbsId[ncurrent][ElemIndex[ncurrent]] = neighbor;
1240 ElemIndex[ncurrent]++;
1241 }
1242
1243 // pickout repeated indexes
1244 for (int ne = 0; ne < nexp; ne++)
1245 {
1246 ElemIndex[ne] = 0;
1247 for (int nb = 0; nb < ElemNeighbsNumb[ne]; nb++)
1248 {
1249 int neighbId = ElemNeighbsId[ne][nb];
1250 bool found = false;
1251 for (int nc = 0; nc < ElemIndex[ne]; nc++)
1252 {
1253 if (ElemNeighbsId[ne][nb] == tmpId[ne][nc])
1254 {
1255 found = true;
1256 }
1257 }
1258 if (!found)
1259 {
1260 tmpId[ne][ElemIndex[ne]] = neighbId;
1261 ElemIndex[ne]++;
1262 }
1263 }
1264 }
1265 ElemNeighbsNumb = ElemIndex;
1266 for (int ne = 0; ne < nexp; ne++)
1267 {
1268 int neighb = ElemNeighbsNumb[ne];
1269 if (neighb > 0)
1270 {
1271 ElemNeighbsId[ne] = Array<OneD, int>(neighb, -1);
1272 Vmath::Vcopy(neighb, tmpId[ne], 1, ElemNeighbsId[ne], 1);
1273 }
1274 }
1275
1276 // check errors
1277 for (int ne = 0; ne < nexp; ne++)
1278 {
1279 for (int nb = 0; nb < ElemNeighbsNumb[ne]; nb++)
1280 {
1281 ASSERTL0((ElemNeighbsId[ne][nb] >= 0) &&
1282 (ElemNeighbsId[ne][nb] <= nexp),
1283 "Element id <0 or >number of total elements")
1284 }
1285 }
1286
1287 m_ElemNeighbsNumb = ElemNeighbsNumb;
1288 m_ElemNeighbsId = ElemNeighbsId;
1289}
1290
1291/**
1292 * @brief Gather the local elemental traces in physical space from
1293 * field using #m_locTraceToFieldMap. Note traces are blocked together
1294 * in similar trace point ordering
1295 *
1296 * @param field Solution field in physical space
1297 * @param faces Resulting local traces.
1298 */
1299void LocTraceToTraceMap::LocTracesFromField(
1300 const Array<OneD, const NekDouble> &field, Array<OneD, NekDouble> faces)
1301{
1302 // The static cast is necessary because m_locTraceToFieldMap should be
1303 // Array<OneD, size_t> ... or at least the same type as
1304 // m_locTraceToFieldMap.size() ...
1305 Vmath::Gathr(static_cast<int>(m_locTraceToFieldMap.size()), field,
1306 m_locTraceToFieldMap, faces);
1307}
1308
1309/**
1310 * @brief Reverse process of LocTracesFromField()
1311 * Add the local traces in physical space to field using
1312 * #m_locTraceToFieldMap.
1313 *
1314 * @param field Solution field in physical space
1315 * @param faces local traces.
1316 */
1317void LocTraceToTraceMap::AddLocTracesToField(
1318 const Array<OneD, const NekDouble> &faces, Array<OneD, NekDouble> &field)
1319{
1320 size_t nfield = field.size();
1321 Array<OneD, NekDouble> tmp{nfield, 0.0};
1322 Vmath::Assmb(m_locTraceToFieldMap.size(), faces.data(),
1323 m_locTraceToFieldMap.data(), tmp.data());
1324 Vmath::Vadd(nfield, tmp, 1, field, 1, field, 1);
1325}
1326
1327/**
1328 * @brief Gather the forwards-oriented local traces in physical space from field
1329 * using #m_locTraceToFieldMap.
1330 *
1331 * @param field Solution field in physical space
1332 * @param faces Resulting local forwards-oriented traces.
1333 */
1339
1340/**
1341 * @brief Reshuffle local elemental traces in physical space so that
1342 * similar faces points are blocked together so they can then be
1343 * interpolated with InterpLocTraceToTrace method.
1344 *
1345 * @param loctrace local traces in physical space
1346 * @param reshuffle traces ordered in reshuffled format of similar patterns.
1347 */
1348void LocTraceToTraceMap::ReshuffleLocTracesForInterp(
1349 const int dir, const Array<OneD, const NekDouble> &loctraces,
1350 Array<OneD, NekDouble> reshuffle)
1351{
1352 ASSERTL1(dir < 2, "option dir out of range, "
1353 " dir=0 is fwd, dir=1 is bwd");
1354
1355 Vmath::Gathr(static_cast<int>(m_locTraceToElmtTraceMap[dir].size()),
1356 loctraces, m_locTraceToElmtTraceMap[dir], reshuffle);
1357}
1358
1359/**
1360 * @brief Unshuffle local elemental traces in physical space from
1361 * similar faces points are blocked together to the local elemental trace format
1362 *
1363 * @param loctrace local traces in physical space
1364 * @param reshuffle traces ordered in reshuffled format of similar patterns.
1365 */
1366void LocTraceToTraceMap::UnshuffleLocTraces(
1367 const int dir, const Array<OneD, const NekDouble> &loctraces,
1368 Array<OneD, NekDouble> unshuffle)
1369{
1370 ASSERTL1(dir < 2, "option dir out of range, "
1371 " dir=0 is fwd, dir=1 is bwd");
1372
1373 if (m_locTraceToElmtTraceMap[dir].size()) // single elemt check
1374 {
1375 Vmath::Scatr(m_locTraceToElmtTraceMap[dir].size(), loctraces,
1376 m_locTraceToElmtTraceMap[dir], unshuffle);
1377 }
1378}
1379
1380void LocTraceToTraceMap::InterpLocTracesToTrace(
1381 const int dir, const Array<OneD, const NekDouble> &loctraces,
1382 Array<OneD, NekDouble> &traces)
1383{
1384 switch (m_expdim)
1385 {
1386 case 1: // Essentially do copy
1387 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(),
1388 loctraces.data(),
1389 m_locInterpTraceToTraceMap[dir].data(), traces.data());
1390 break;
1391 case 2:
1392 InterpLocEdgesToTrace(dir, loctraces, traces);
1393 break;
1394 case 3:
1395 InterpLocFacesToTrace(dir, loctraces, traces);
1396 break;
1397 default:
1398 NEKERROR(ErrorUtil::efatal, "Not set up");
1399 break;
1400 }
1401}
1402
1403/**
1404 * @brief Interpolate local trace edges to global trace edge point distributions
1405 * where required.
1406 *
1407 * @param dir Selects forwards (0) or backwards (1) direction.
1408 * @param locfaces Local trace edge storage.
1409 * @param faces Global trace edge storage
1410 */
1411void LocTraceToTraceMap::InterpLocEdgesToTrace(
1412 const int dir, const Array<OneD, const NekDouble> &locedges,
1413 Array<OneD, NekDouble> &edges)
1414{
1415 ASSERTL1(dir < 2, "option dir out of range, "
1416 " dir=0 is fwd, dir=1 is bwd");
1417
1418 int cnt = 0;
1419 int cnt1 = 0;
1420
1421 // tmp space assuming forward map is of size of trace
1422 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
1423
1424 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1425 {
1426 // Check if there are edges to interpolate
1427 if (m_interpNtraces[dir][i])
1428 {
1429 // Get to/from points
1430 LibUtilities::PointsKey fromPointsKey0 =
1431 std::get<0>(m_interpPoints[dir][i]);
1432 LibUtilities::PointsKey toPointsKey0 =
1433 std::get<2>(m_interpPoints[dir][i]);
1434
1435 int fnp = fromPointsKey0.GetNumPoints();
1436 int tnp = toPointsKey0.GetNumPoints();
1437 int nedges = m_interpNtraces[dir][i];
1438
1439 // Do interpolation here if required
1440 switch (m_interpTrace[dir][i])
1441 {
1442 case eNoInterp: // Just copy
1443 {
1444 if (fnp == tnp) // identical points : just copy all points
1445 {
1446 Vmath::Vcopy(nedges * fnp, locedges.data() + cnt, 1,
1447 tmp.data() + cnt1, 1);
1448 }
1449 else // copy interior points and fill in end points (if any)
1450 {
1451 int fbegin, fend, fsize;
1452 int tbegin, tend, tsize;
1453 LibUtilities::PointsKey::GetEffectiveQuadRange(
1454 fromPointsKey0, fbegin, fend);
1455 LibUtilities::PointsKey::GetEffectiveQuadRange(
1456 toPointsKey0, tbegin, tend);
1457 fsize = fend - fbegin;
1458 tsize = tend - tbegin;
1459 ASSERTL0(
1460 fsize == tsize,
1461 "Quad ranges mismatch in InterpLocEdgesToTrace!");
1462 for (int k = 0; k < nedges; ++k)
1463 {
1464 Vmath::Vcopy(
1465 fsize, locedges.data() + cnt + fnp * k + fbegin,
1466 1, tmp.data() + cnt1 + tnp * k + tbegin, 1);
1467 // fill 0 ~ tbegin with tbegin
1468 Vmath::Fill(tbegin, tmp[cnt1 + tnp * k + tbegin],
1469 tmp.data() + cnt1 + tnp * k, 1);
1470 // fill tend ~ tnp with tend
1471 Vmath::Fill(tnp - tend,
1472 tmp[cnt1 + tnp * k + tend - 1],
1473 tmp.data() + cnt1 + tnp * k + tend, 1);
1474 }
1475 }
1476 }
1477 break;
1478 case eInterpDir0:
1479 {
1480 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1481 Blas::Dgemm('N', 'N', tnp, nedges, fnp, 1.0,
1482 I0->GetPtr().data(), tnp, locedges.data() + cnt,
1483 fnp, 0.0, tmp.data() + cnt1, tnp);
1484 }
1485 break;
1486 case eInterpEndPtDir0:
1487 {
1488 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1489
1490 for (int k = 0; k < nedges; ++k)
1491 {
1492 Vmath::Vcopy(fnp, &locedges[cnt + k * fnp], 1,
1493 &tmp[cnt1 + k * tnp], 1);
1494
1495 tmp[cnt1 + k * tnp + tnp - 1] = Vmath::Dot(
1496 fnp, locedges.data() + cnt + k * fnp, 1, &I0[0], 1);
1497 }
1498 }
1499 break;
1500 default:
1501 NEKERROR(ErrorUtil::efatal,
1502 "Invalid interpolation type for 2D elements");
1503 break;
1504 }
1505
1506 cnt += nedges * fnp;
1507 cnt1 += nedges * tnp;
1508 }
1509 }
1510
1511 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(), tmp.data(),
1512 m_locInterpTraceToTraceMap[dir].data(), edges.data());
1513}
1514
1515/**
1516 * @brief Interpolate local faces to trace face point distributions where
1517 * required.
1518 *
1519 * @param dir Selects forwards (0) or backwards (1) direction.
1520 * @param locfaces Local trace face storage.
1521 * @param faces Global trace face storage
1522 */
1523void LocTraceToTraceMap::InterpLocFacesToTrace(
1524 const int dir, const Array<OneD, const NekDouble> &locfaces,
1525 Array<OneD, NekDouble> faces)
1526{
1527 ASSERTL1(dir < 2, "option dir out of range, "
1528 " dir=0 is fwd, dir=1 is bwd");
1529
1530 int cnt1 = 0;
1531 int cnt = 0;
1532
1533 // tmp space assuming forward map is of size of trace
1534 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
1535
1536 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1537 {
1538 // Check if there are faces to interpolate
1539 if (m_interpNtraces[dir][i])
1540 {
1541 // Get to/from points
1542 LibUtilities::PointsKey fromPointsKey0 =
1543 std::get<0>(m_interpPoints[dir][i]);
1544 LibUtilities::PointsKey fromPointsKey1 =
1545 std::get<1>(m_interpPoints[dir][i]);
1546 LibUtilities::PointsKey toPointsKey0 =
1547 std::get<2>(m_interpPoints[dir][i]);
1548 LibUtilities::PointsKey toPointsKey1 =
1549 std::get<3>(m_interpPoints[dir][i]);
1550
1551 int fnp0 = fromPointsKey0.GetNumPoints();
1552 int fnp1 = fromPointsKey1.GetNumPoints();
1553 int tnp0 = toPointsKey0.GetNumPoints();
1554 int tnp1 = toPointsKey1.GetNumPoints();
1555 int nfaces = m_interpNtraces[dir][i];
1556 int nfromfacepts = fnp0 * fnp1;
1557 int ntofacepts = tnp0 * tnp1;
1558
1559 // Do interpolation here if required
1560 switch (m_interpTrace[dir][i])
1561 {
1562 case eNoInterp: // Just copy
1563 {
1564 if (fnp0 == tnp0 && fnp1 == tnp1)
1565 {
1566 Vmath::Vcopy(nfaces * nfromfacepts,
1567 locfaces.data() + cnt, 1,
1568 tmp.data() + cnt1, 1);
1569 }
1570 else
1571 {
1572 int fbegin0, fend0, fbegin1, fend1;
1573 int tbegin0, tend0, tbegin1, tend1;
1574 LibUtilities::PointsKey::GetEffectiveQuadRange(
1575 fromPointsKey0, fbegin0, fend0);
1576 LibUtilities::PointsKey::GetEffectiveQuadRange(
1577 toPointsKey0, tbegin0, tend0);
1578 LibUtilities::PointsKey::GetEffectiveQuadRange(
1579 fromPointsKey1, fbegin1, fend1);
1580 LibUtilities::PointsKey::GetEffectiveQuadRange(
1581 toPointsKey1, tbegin1, tend1);
1582 ASSERTL0(
1583 fend0 - fbegin0 == tend0 - tbegin0 &&
1584 fend1 - fbegin1 == tend1 - tbegin1,
1585 "Quad ranges mismatch in InterpLocFacesToTrace!");
1586 for (int j = 0; j < nfaces; ++j)
1587 {
1588 for (int k = fbegin1, l = tbegin1; k < fend1;
1589 ++k, ++l)
1590 {
1591 Vmath::Vcopy(
1592 fend0 - fbegin0,
1593 locfaces.data() + cnt + j * nfromfacepts +
1594 k * fnp0 + fbegin0,
1595 1,
1596 tmp.data() + cnt1 + j * ntofacepts +
1597 l * tnp0 + tbegin0,
1598 1);
1599 }
1600 }
1601 }
1602 }
1603 break;
1604 case eInterpDir0:
1605 {
1606 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1607 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1608 I0->GetPtr().data(), tnp0,
1609 locfaces.data() + cnt, fnp0, 0.0,
1610 tmp.data() + cnt1, tnp0);
1611 }
1612 break;
1613 case eInterpEndPtDir0:
1614 {
1615 int nfaces = m_interpNtraces[dir][i];
1616 for (int k = 0; k < fnp0; ++k)
1617 {
1618 Vmath::Vcopy(nfaces * fnp1, locfaces.data() + cnt + k,
1619 fnp0, tmp.data() + cnt1 + k, tnp0);
1620 }
1621
1622 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1623 Blas::Dgemv('T', fnp0, tnp1 * nfaces, 1.0,
1624 tmp.data() + cnt1, tnp0, I0.data(), 1, 0.0,
1625 tmp.data() + cnt1 + tnp0 - 1, tnp0);
1626 }
1627 break;
1628 case eInterpDir1:
1629 {
1630 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1631 for (int j = 0; j < nfaces; ++j)
1632 {
1633 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
1634 locfaces.data() + cnt + j * fnp0 * fnp1,
1635 tnp0, I1->GetPtr().data(), tnp1, 0.0,
1636 tmp.data() + cnt1 + j * tnp0 * tnp1, tnp0);
1637 }
1638 }
1639 break;
1640 case eInterpEndPtDir1:
1641 {
1642 Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1643 for (int j = 0; j < nfaces; ++j)
1644 {
1645 // copy all points
1646 Vmath::Vcopy(fnp0 * fnp1,
1647 locfaces.data() + cnt + j * fnp0 * fnp1, 1,
1648 tmp.data() + cnt1 + j * tnp0 * tnp1, 1);
1649
1650 // interpolate end points
1651 for (int k = 0; k < tnp0; ++k)
1652 {
1653 tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1654 Vmath::Dot(fnp1,
1655 locfaces.data() + cnt +
1656 j * fnp0 * fnp1 + k,
1657 fnp0, &I1[0], 1);
1658 }
1659 }
1660 }
1661 break;
1662 case eInterpBothDirs:
1663 {
1664 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1665 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1666 Array<OneD, NekDouble> wsp(nfaces * fnp0 * tnp1);
1667
1668 for (int j = 0; j < nfaces; ++j)
1669 {
1670 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1671 locfaces.data() + cnt + j * fnp0 * fnp1,
1672 fnp0, I1->GetPtr().data(), tnp1, 0.0,
1673 wsp.data() + j * fnp0 * tnp1, fnp0);
1674 }
1675 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
1676 I0->GetPtr().data(), tnp0, wsp.data(), fnp0,
1677 0.0, tmp.data() + cnt1, tnp0);
1678 }
1679 break;
1680 case eInterpEndPtDir0InterpDir1:
1681 {
1682 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1683 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1684
1685 for (int j = 0; j < nfaces; ++j)
1686 {
1687 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
1688 locfaces.data() + cnt + j * fnp0 * fnp1,
1689 fnp0, I1->GetPtr().data(), tnp1, 0.0,
1690 tmp.data() + cnt1 + j * tnp0 * tnp1, tnp0);
1691 }
1692
1693 Blas::Dgemv('T', fnp0, tnp1 * nfaces, 1.0,
1694 tmp.data() + cnt1, tnp0, I0.data(), 1, 0.0,
1695 tmp.data() + cnt1 + tnp0 - 1, tnp0);
1696 }
1697 break;
1698 default:
1699 ASSERTL0(false, "Interplation case needs implementing");
1700 break;
1701 }
1702 cnt += nfaces * nfromfacepts;
1703 cnt1 += nfaces * ntofacepts;
1704 }
1705 }
1706
1707 Vmath::Scatr(m_locInterpTraceToTraceMap[dir].size(), tmp.data(),
1708 m_locInterpTraceToTraceMap[dir].data(), faces.data());
1709}
1710
1711void LocTraceToTraceMap::InterpLocTracesToTraceTranspose(
1712 const int dir, const Array<OneD, const NekDouble> &trace,
1713 Array<OneD, NekDouble> &loctrace)
1714{
1715 switch (m_expdim)
1716 {
1717 case 2:
1718 InterpLocEdgesToTraceTranspose(dir, trace, loctrace);
1719 break;
1720 case 3:
1721 InterpLocFacesToTraceTranspose(dir, trace, loctrace);
1722 break;
1723 default:
1724 NEKERROR(ErrorUtil::efatal, "Not set up");
1725 break;
1726 }
1727}
1728
1729/**
1730 * @brief Transpose of Interp local edges to trace
1731 *
1732 * @param dir Selects forwards (0) or backwards (1) direction.
1733 * @param edges Global trace edge
1734 * @param locedges Local trace edge
1735 */
1736void LocTraceToTraceMap::InterpLocEdgesToTraceTranspose(
1737 const int dir, const Array<OneD, const NekDouble> &edges,
1738 Array<OneD, NekDouble> &locedges)
1739{
1740 ASSERTL1(dir < 2, "option dir out of range, "
1741 " dir=0 is fwd, dir=1 is bwd");
1742
1743 int cnt = 0;
1744 int cnt1 = 0;
1745
1746 // tmp space assuming forward map is of size of trace
1747 Array<OneD, NekDouble> tmp{size_t(m_nTracePts), 0.0};
1748 Vmath::Gathr((int)m_locInterpTraceToTraceMap[dir].size(), edges.data(),
1749 m_locInterpTraceToTraceMap[dir].data(), tmp.data());
1750
1751 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1752 {
1753 // Check if there are edges to interpolate
1754 if (m_interpNtraces[dir][i])
1755 {
1756 // Get to/from points
1757 LibUtilities::PointsKey fromPointsKey0 =
1758 std::get<0>(m_interpPoints[dir][i]);
1759 LibUtilities::PointsKey toPointsKey0 =
1760 std::get<2>(m_interpPoints[dir][i]);
1761
1762 int fnp = fromPointsKey0.GetNumPoints();
1763 int tnp = toPointsKey0.GetNumPoints();
1764 int nedges = m_interpNtraces[dir][i];
1765
1766 // Do interpolation here if required
1767 switch (m_interpTrace[dir][i])
1768 {
1769 case eNoInterp: // Just copy
1770 {
1771 int fbegin, fsize;
1772 int tbegin, tsize;
1773 LibUtilities::PointsKey::GetEffectiveQuadRange(
1774 fromPointsKey0, fbegin, fsize);
1775 LibUtilities::PointsKey::GetEffectiveQuadRange(
1776 toPointsKey0, tbegin, tsize);
1777 fsize = fsize - fbegin;
1778 tsize = tsize - tbegin;
1779 ASSERTL0(fsize == tsize, "Quad ranges mismatch in "
1780 "InterpLocEdgesToTraceTranspose!");
1781 for (int k = 0; k < nedges; ++k)
1782 {
1783 Vmath::Zero(fnp, locedges.data() + cnt + fnp * k, 1);
1784 Vmath::Vcopy(
1785 fsize, tmp.data() + cnt1 + tnp * k + tbegin, 1,
1786 locedges.data() + cnt + fnp * k + fbegin, 1);
1787 }
1788 }
1789 break;
1790 case eInterpDir0:
1791 {
1792 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1793 Blas::Dgemm('T', 'N', fnp, nedges, tnp, 1.0,
1794 I0->GetPtr().data(), tnp, tmp.data() + cnt1,
1795 tnp, 0.0, locedges.data() + cnt, fnp);
1796 }
1797 break;
1798 case eInterpEndPtDir0:
1799 {
1800 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1801
1802 for (int k = 0; k < nedges; ++k)
1803 {
1804 Vmath::Vcopy(fnp, &tmp[cnt1 + k * tnp], 1,
1805 &locedges[cnt + k * fnp], 1);
1806
1807 Vmath::Svtvp(fnp, tmp[cnt1 + k * tnp + tnp - 1], &I0[0],
1808 1, locedges.data() + cnt + k * fnp, 1,
1809 locedges.data() + cnt + k * fnp, 1);
1810 }
1811 }
1812 break;
1813 default:
1814 NEKERROR(ErrorUtil::efatal,
1815 "Invalid interpolation type for 2D elements");
1816 break;
1817 }
1818
1819 cnt += nedges * fnp;
1820 cnt1 += nedges * tnp;
1821 }
1822 }
1823}
1824
1825/**
1826 * @brief transpose of interp local faces to trace
1827 *
1828 * @param dir Selects forwards (0) or backwards (1) direction.
1829 * @param loctraces Local trace
1830 * @param traces trace .
1831 */
1832void LocTraceToTraceMap::InterpLocFacesToTraceTranspose(
1833 const int dir, const Array<OneD, const NekDouble> &traces,
1834 Array<OneD, NekDouble> &loctraces)
1835{
1836 ASSERTL1(dir < 2, "option dir out of range, "
1837 " dir=0 is fwd, dir=1 is bwd");
1838
1839 int cnt = 0;
1840 int cnt1 = 0;
1841
1842 // tmp space assuming forward map is of size of trace
1843 Array<OneD, NekDouble> tmp{size_t(m_nTracePts), 0.0};
1844 // The static cast is necessary because m_locInterpTraceToTraceMap should be
1845 // Array<OneD, size_t> ... or at least the same type as
1846 // m_locInterpTraceToTraceMap.size() ...
1847 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
1848 traces.data(), m_locInterpTraceToTraceMap[dir].data(),
1849 tmp.data());
1850
1851 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
1852 {
1853 // Check if there are elementboundaries to interpolate
1854 if (m_interpNtraces[dir][i])
1855 {
1856 // Get to/from points
1857 LibUtilities::PointsKey fromPointsKey0 =
1858 std::get<0>(m_interpPoints[dir][i]);
1859 LibUtilities::PointsKey fromPointsKey1 =
1860 std::get<1>(m_interpPoints[dir][i]);
1861 LibUtilities::PointsKey toPointsKey0 =
1862 std::get<2>(m_interpPoints[dir][i]);
1863 LibUtilities::PointsKey toPointsKey1 =
1864 std::get<3>(m_interpPoints[dir][i]);
1865 // Here the f(from) and t(to) are chosen to be consistent with
1866 // InterpLocFacesToTrace
1867 int fnp0 = fromPointsKey0.GetNumPoints();
1868 int fnp1 = fromPointsKey1.GetNumPoints();
1869 int tnp0 = toPointsKey0.GetNumPoints();
1870 int tnp1 = toPointsKey1.GetNumPoints();
1871 int nfaces = m_interpNtraces[dir][i];
1872 int nfromfacepts = fnp0 * fnp1;
1873 int ntofacepts = tnp0 * tnp1;
1874
1875 // Do transpose interpolation here if required
1876 switch (m_interpTrace[dir][i])
1877 {
1878 case eNoInterp: // Just copy
1879 {
1880 int fbegin0, fend0, fbegin1, fend1;
1881 int tbegin0, tend0, tbegin1, tend1;
1882 LibUtilities::PointsKey::GetEffectiveQuadRange(
1883 fromPointsKey0, fbegin0, fend0);
1884 LibUtilities::PointsKey::GetEffectiveQuadRange(
1885 toPointsKey0, tbegin0, tend0);
1886 LibUtilities::PointsKey::GetEffectiveQuadRange(
1887 fromPointsKey1, fbegin1, fend1);
1888 LibUtilities::PointsKey::GetEffectiveQuadRange(
1889 toPointsKey1, tbegin1, tend1);
1890 // overwrite end by size
1891 fend0 = fend0 - fbegin0;
1892 tend0 = tend0 - tbegin0;
1893 ASSERTL0(fend0 == tend0 &&
1894 fend1 - fbegin1 == tend1 - tbegin1,
1895 "Quad ranges mismatch in "
1896 "InterpLocFacesToTraceTranspose!");
1897 for (int j = 0; j < nfaces; ++j)
1898 {
1899 for (int k = fbegin1, l = tbegin1; k < fend1; ++k, ++l)
1900 {
1901 Vmath::Vcopy(fend0,
1902 tmp.data() + cnt1 + j * ntofacepts +
1903 l * tnp0 + tbegin0,
1904 1,
1905 loctraces.data() + cnt +
1906 j * nfromfacepts + k * fnp0 +
1907 fbegin0,
1908 1);
1909 }
1910 }
1911 }
1912 break;
1913 case eInterpDir0:
1914 {
1915 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1916 Blas::Dgemm('T', 'N', fnp0, tnp1, tnp0, 1.0,
1917 I0->GetPtr().data(), tnp0, tmp.data() + cnt1,
1918 tnp0, 0.0, loctraces.data() + cnt, fnp0);
1919 }
1920 break;
1921 case eInterpEndPtDir0:
1922 {
1923 for (int k = 0; k < fnp0; ++k)
1924 {
1925 Vmath::Vcopy(nfaces * fnp1, tmp.data() + cnt1 + k, tnp0,
1926 loctraces.data() + cnt + k, fnp0);
1927 }
1928
1929 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1930 for (int k = 0; k < tnp1 * nfaces; k++)
1931 {
1932 Vmath::Svtvp(fnp0, tmp[cnt1 + tnp0 - 1 + k * tnp0],
1933 &I0[0], 1,
1934 loctraces.data() + cnt + k * fnp0, 1,
1935 loctraces.data() + cnt + k * fnp0, 1);
1936 }
1937 }
1938 break;
1939 case eInterpDir1:
1940 {
1941 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1942
1943 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1944 {
1945 Blas::Dgemm('N', 'N', tnp0, fnp1, tnp1, 1.0,
1946 tmp.data() + cnt1 + j * tnp0 * tnp1, tnp0,
1947 I1->GetPtr().data(), tnp1, 0.0,
1948 loctraces.data() + cnt + j * fnp0 * fnp1,
1949 tnp0);
1950 }
1951 }
1952 break;
1953 case eInterpEndPtDir1:
1954 {
1955 Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1956 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1957 {
1958 Vmath::Vcopy(
1959 fnp0 * fnp1, tmp.data() + cnt1 + j * tnp0 * tnp1, 1,
1960 loctraces.data() + cnt + j * fnp0 * fnp1, 1);
1961
1962 for (int k = 0; k < fnp1; k++)
1963 {
1964 Vmath::Svtvp(
1965 fnp0, I1[k],
1966 &tmp[cnt1 + (j + 1) * tnp0 * tnp1 - tnp0], 1,
1967 &loctraces[cnt + j * fnp0 * fnp1 + k * fnp0], 1,
1968 &loctraces[cnt + j * fnp0 * fnp1 + k * fnp0],
1969 1);
1970 }
1971 }
1972 }
1973 break;
1974 case eInterpBothDirs:
1975 {
1976 DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1977 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1978
1979 Array<OneD, NekDouble> wsp{
1980 size_t(m_interpNtraces[dir][i] * fnp0 * tnp1)};
1981
1982 Blas::Dgemm('T', 'N', fnp0, tnp1 * m_interpNtraces[dir][i],
1983 tnp0, 1.0, I0->GetPtr().data(), tnp0,
1984 tmp.data() + cnt1, tnp0, 0.0, wsp.data(), fnp0);
1985
1986 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
1987 {
1988 Blas::Dgemm('N', 'N', fnp0, fnp1, tnp1, 1.0,
1989 wsp.data() + j * fnp0 * tnp1, fnp0,
1990 I1->GetPtr().data(), tnp1, 0.0,
1991 loctraces.data() + cnt + j * fnp0 * fnp1,
1992 fnp0);
1993 }
1994 }
1995 break;
1996 case eInterpEndPtDir0InterpDir1:
1997 {
1998 DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1999 Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
2000
2001 Array<OneD, NekDouble> wsp{
2002 size_t(m_interpNtraces[dir][i] * fnp0 * tnp1)};
2003
2004 for (int k = 0; k < tnp1 * m_interpNtraces[dir][i]; k++)
2005 {
2006 Vmath::Svtvp(fnp0, tmp[cnt1 + tnp0 - 1 + k * tnp0],
2007 &I0[0], 1, tmp.data() + cnt1 + k * tnp0, 1,
2008 wsp.data() + k * fnp0, 1);
2009 }
2010
2011 for (int j = 0; j < m_interpNtraces[dir][i]; ++j)
2012 {
2013 Blas::Dgemm('N', 'N', fnp0, fnp1, tnp1, 1.0,
2014 wsp.data() + j * fnp0 * tnp1, fnp0,
2015 I1->GetPtr().data(), tnp1, 0.0,
2016 loctraces.data() + cnt + j * fnp0 * fnp1,
2017 fnp0);
2018 }
2019 }
2020 break;
2021 }
2022 cnt += nfaces * nfromfacepts;
2023 cnt1 += nfaces * ntofacepts;
2024 }
2025 }
2026}
2027
2028void LocTraceToTraceMap::InterpTraceToLocTrace(
2029 const int dir, const Array<OneD, NekDouble> &traces,
2030 Array<OneD, NekDouble> &loctraces)
2031{
2032 switch (m_expdim)
2033 {
2034 case 1: // Essentially do copy
2035 Vmath::Gathr(
2036 static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
2037 traces.data(), m_locInterpTraceToTraceMap[dir].data(),
2038 loctraces.data());
2039 break;
2040 case 2:
2041 InterpTraceToLocEdges(dir, traces, loctraces);
2042 break;
2043 case 3:
2044 InterpTraceToLocFaces(dir, traces, loctraces);
2045 break;
2046 default:
2047 ASSERTL0(false, "Not set up");
2048 break;
2049 }
2050}
2051
2052/**
2053 * @brief Interpolate global trace edge to local trace edges point
2054 * distributions where required.
2055 *
2056 * @param dir Selects forwards (0) or backwards (1) direction.
2057 * @param locfaces Local trace edge storage.
2058 * @param faces Global trace edge storage
2059 */
2060void LocTraceToTraceMap::InterpTraceToLocEdges(
2061 const int dir, const Array<OneD, const NekDouble> &edges,
2062 Array<OneD, NekDouble> &locedges)
2063{
2064 ASSERTL1(dir < 2, "option dir out of range, "
2065 " dir=0 is fwd, dir=1 is bwd");
2066
2067 int cnt = 0;
2068 int cnt1 = 0;
2069
2070 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
2071
2072 // unshuffles trace into lcoally orientated format.
2073 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
2074 edges.data(), m_locInterpTraceToTraceMap[dir].data(),
2075 tmp.data());
2076
2077 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
2078 {
2079 // Check if there are edges to interpolate
2080 if (m_interpNtraces[dir][i])
2081 {
2082 // Get to/from points
2083 LibUtilities::PointsKey fromPointsKey0 =
2084 std::get<2>(m_interpPoints[dir][i]);
2085 LibUtilities::PointsKey toPointsKey0 =
2086 std::get<0>(m_interpPoints[dir][i]);
2087
2088 int fnp = fromPointsKey0.GetNumPoints();
2089 int tnp = toPointsKey0.GetNumPoints();
2090 int nedges = m_interpNtraces[dir][i];
2091
2092 // Do interpolation here if required
2093 switch (m_interpTrace[dir][i])
2094 {
2095 case eNoInterp: // Just copy
2096 {
2097 int fbegin, fend;
2098 int tbegin, tend;
2099 LibUtilities::PointsKey::GetEffectiveQuadRange(
2100 fromPointsKey0, fbegin, fend);
2101 LibUtilities::PointsKey::GetEffectiveQuadRange(
2102 toPointsKey0, tbegin, tend);
2103 // overwrite fend by size
2104 fend = fend - fbegin;
2105 tend = tend - tbegin;
2106 ASSERTL0(fend == tend,
2107 "Quad ranges mismatch in InterpTraceToLocEdges!");
2108 for (int k = 0; k < nedges; ++k)
2109 {
2110 Vmath::Zero(tnp, locedges.data() + cnt1 + tnp * k, 1);
2111 Vmath::Vcopy(
2112 fend, tmp.data() + cnt + fnp * k + fbegin, 1,
2113 locedges.data() + cnt1 + tnp * k + tbegin, 1);
2114 }
2115 }
2116 break;
2117 case eInterpDir0:
2118 {
2119 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
2120 Blas::Dgemm('N', 'N', tnp, nedges, fnp, 1.0,
2121 I0->GetPtr().data(), tnp, tmp.data() + cnt, fnp,
2122 0.0, locedges.data() + cnt1, tnp);
2123 }
2124 break;
2125 case eInterpEndPtDir0:
2126 {
2127 // Just copy points back
2128 for (int k = 0; k < nedges; ++k)
2129 {
2130 // Should be tnp rather than fnp on first argument.
2131 Vmath::Vcopy(tnp, &tmp[cnt + k * fnp], 1,
2132 &locedges[cnt1 + k * tnp], 1);
2133 }
2134 }
2135 break;
2136 default:
2137 ASSERTL0(false,
2138 "Invalid interpolation type for 2D elements");
2139 break;
2140 }
2141
2142 cnt += nedges * fnp;
2143 cnt1 += nedges * tnp;
2144 }
2145 }
2146}
2147
2148/**
2149 * @brief Interpolate global trace edge to local trace edges point
2150 * distributions where required.
2151 *
2152 * @param dir Selects forwards (0) or backwards (1) direction.
2153 * @param locfaces Local trace edge storage.
2154 * @param faces Global trace edge storage
2155 */
2156void LocTraceToTraceMap::InterpTraceToLocFaces(
2157 const int dir, const Array<OneD, const NekDouble> &faces,
2158 Array<OneD, NekDouble> &locfaces)
2159{
2160 ASSERTL1(dir < 2, "option dir out of range, "
2161 " dir=0 is fwd, dir=1 is bwd");
2162
2163 int cnt = 0;
2164 int cnt1 = 0;
2165
2166 Array<OneD, NekDouble> tmp(m_nTracePts, 0.0);
2167
2168 // unshuffles trace into lcoally orientated format.
2169 Vmath::Gathr(static_cast<int>(m_locInterpTraceToTraceMap[dir].size()),
2170 faces.data(), m_locInterpTraceToTraceMap[dir].data(),
2171 tmp.data());
2172
2173 for (int i = 0; i < m_interpTrace[dir].size(); ++i)
2174 {
2175 // Check if there are faces to interpolate
2176 if (m_interpNtraces[dir][i])
2177 {
2178 // Get to/from points
2179 LibUtilities::PointsKey fromPointsKey0 =
2180 std::get<2>(m_interpPoints[dir][i]);
2181 LibUtilities::PointsKey fromPointsKey1 =
2182 std::get<3>(m_interpPoints[dir][i]);
2183 LibUtilities::PointsKey toPointsKey0 =
2184 std::get<0>(m_interpPoints[dir][i]);
2185 LibUtilities::PointsKey toPointsKey1 =
2186 std::get<1>(m_interpPoints[dir][i]);
2187
2188 int fnp0 = fromPointsKey0.GetNumPoints();
2189 int fnp1 = fromPointsKey1.GetNumPoints();
2190 int tnp0 = toPointsKey0.GetNumPoints();
2191 int tnp1 = toPointsKey1.GetNumPoints();
2192 int nfaces = m_interpNtraces[dir][i];
2193 int nfromfacepts = fnp0 * fnp1;
2194 int ntofacepts = tnp0 * tnp1;
2195
2196 // Do interpolation here if required
2197 switch (m_interpTrace[dir][i])
2198 {
2199 case eNoInterp: // Just copy
2200 {
2201 int fbegin0, fend0, fbegin1, fend1;
2202 int tbegin0, tend0, tbegin1, tend1;
2203 LibUtilities::PointsKey::GetEffectiveQuadRange(
2204 fromPointsKey0, fbegin0, fend0);
2205 LibUtilities::PointsKey::GetEffectiveQuadRange(
2206 toPointsKey0, tbegin0, tend0);
2207 LibUtilities::PointsKey::GetEffectiveQuadRange(
2208 fromPointsKey1, fbegin1, fend1);
2209 LibUtilities::PointsKey::GetEffectiveQuadRange(
2210 toPointsKey1, tbegin1, tend1);
2211 // overwrite end by size
2212 fend0 = fend0 - fbegin0;
2213 tend0 = tend0 - tbegin0;
2214 ASSERTL0(fend0 == tend0 &&
2215 fend1 - fbegin1 == tend1 - tbegin1,
2216 "Quad ranges mismatch in InterpTraceToLocFaces!");
2217 for (int j = 0; j < nfaces; ++j)
2218 {
2219 Vmath::Zero(ntofacepts,
2220 locfaces.data() + cnt1 + ntofacepts * j, 1);
2221 for (int k = fbegin1, l = tbegin1; k < fend1; ++k, ++l)
2222 {
2223 Vmath::Vcopy(fend0,
2224 tmp.data() + cnt + j * nfromfacepts +
2225 k * fnp0 + fbegin0,
2226 1,
2227 locfaces.data() + cnt1 +
2228 ntofacepts * j + l * tnp0 +
2229 tbegin0,
2230 1);
2231 }
2232 }
2233 }
2234 break;
2235 case eInterpDir0:
2236 {
2237 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
2238 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
2239 I0->GetPtr().data(), tnp0, tmp.data() + cnt,
2240 fnp0, 0.0, locfaces.data() + cnt1, tnp0);
2241 }
2242 break;
2243 case eInterpDir1:
2244 {
2245 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
2246 for (int j = 0; j < nfaces; ++j)
2247 {
2248 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
2249 tmp.data() + cnt + j * fnp0 * fnp1, tnp0,
2250 I1->GetPtr().data(), tnp1, 0.0,
2251 locfaces.data() + cnt1 + j * tnp0 * tnp1,
2252 tnp0);
2253 }
2254 }
2255 break;
2256 case eInterpEndPtDir0:
2257 {
2258 for (int j = 0; j < nfaces * tnp1; ++j)
2259 {
2260 Vmath::Vcopy(tnp0, tmp.data() + cnt + j * fnp0, 1,
2261 locfaces.data() + cnt1 + j * tnp0, 1);
2262 }
2263 }
2264 break;
2265 case eInterpEndPtDir1:
2266 {
2267 for (int j = 0; j < nfaces; ++j)
2268 {
2269 // copy all points missing off top verex in dir 1
2270 Vmath::Vcopy(
2271 tnp0 * tnp1, tmp.data() + cnt + j * fnp0 * fnp1, 1,
2272 locfaces.data() + cnt1 + j * tnp0 * tnp1, 1);
2273 }
2274 }
2275 break;
2276 case eInterpBothDirs:
2277 {
2278 DNekMatSharedPtr I0 = m_interpFromTraceI0[dir][i];
2279 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
2280 Array<OneD, NekDouble> wsp(nfaces * fnp0 * tnp1 * fnp0);
2281
2282 for (int j = 0; j < nfaces; ++j)
2283 {
2284 Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
2285 tmp.data() + cnt + j * fnp0 * fnp1, fnp0,
2286 I1->GetPtr().data(), tnp1, 0.0,
2287 wsp.data() + j * fnp0 * tnp1, fnp0);
2288 }
2289
2290 Blas::Dgemm('N', 'N', tnp0, tnp1 * nfaces, fnp0, 1.0,
2291 I0->GetPtr().data(), tnp0, wsp.data(), fnp0,
2292 0.0, locfaces.data() + cnt1, tnp0);
2293 }
2294 break;
2295 case eInterpEndPtDir0InterpDir1:
2296 {
2297 DNekMatSharedPtr I1 = m_interpFromTraceI1[dir][i];
2298 for (int j = 0; j < nfaces; ++j)
2299 {
2300 Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0,
2301 tmp.data() + cnt + j * fnp0 * fnp1, fnp0,
2302 I1->GetPtr().data(), tnp1, 0.0,
2303 locfaces.data() + cnt1 + j * tnp0 * tnp1,
2304 tnp0);
2305 }
2306 }
2307 break;
2308 default:
2309 ASSERTL0(false, "Interpolation case not implemneted (yet)");
2310 break;
2311 }
2312 cnt += nfaces * nfromfacepts;
2313 cnt1 += nfaces * ntofacepts;
2314 }
2315 }
2316}
2317
2318/**
2319 * @brief Add contributions from trace coefficients to the elemental field
2320 * storage.
2321 *
2322 * @param trace Array of global trace coefficients.
2323 * @param field Array containing field coefficients storage.
2324 */
2325void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
2326 const Array<OneD, const NekDouble> &trace, Array<OneD, NekDouble> &field)
2327{
2328 int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
2329 for (int i = 0; i < nvals; ++i)
2330 {
2331 field[m_traceCoeffsToElmtMap[0][i]] +=
2332 m_traceCoeffsToElmtSign[0][i] *
2333 trace[m_traceCoeffsToElmtTrace[0][i]];
2334 }
2335}
2336
2337/**
2338 * @brief Add contributions from backwards or forwards oriented trace
2339 * coefficients to the elemental field storage.
2340 *
2341 * @param dir Selects forwards (0) or backwards (1) direction
2342 * @param trace Array of global trace coefficients.
2343 * @param field Array containing field coefficients storage.
2344 */
2345void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
2346 const int dir, const Array<OneD, const NekDouble> &trace,
2347 Array<OneD, NekDouble> &field)
2348{
2349 int nvals = m_nTraceCoeffs[dir];
2350 for (int i = 0; i < nvals; ++i)
2351 {
2352 field[m_traceCoeffsToElmtMap[dir][i]] +=
2353 m_traceCoeffsToElmtSign[dir][i] *
2354 trace[m_traceCoeffsToElmtTrace[dir][i]];
2355 }
2356}
2357
2358} // namespace Nektar::MultiRegions
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition ErrorUtil.hpp:210
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition ErrorUtil.hpp:250
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::GetEffectiveQuadRange
static void GetEffectiveQuadRange(const LibUtilities::PointsKey &pkey, int &q_begin, int &q_end)
offset the m_nq_begin and m_nq_end according to the points type
Definition Points.h:219
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:98
Nektar::MultiRegions::ExpList::GetCollections
const Collections::CollectionVector & GetCollections() const
This function returns collections.
Definition ExpList.h:1091
Nektar::MultiRegions::ExpList::GetCoeffsToElmt
const Array< OneD, const std::pair< int, int > > & GetCoeffsToElmt() const
Get m_coeffs to elemental value map.
Definition ExpList.h:2237
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:2197
Nektar::MultiRegions::ExpList::GetShapeDimension
int GetShapeDimension()
This function returns the dimension of the shape of the element eid.
Definition ExpList.h:2014
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:2189
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:2204
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:467
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:457
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:446
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:413
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:1736
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocTracesToTraceTranspose
void InterpLocTracesToTraceTranspose(const int dir, const Array< OneD, const NekDouble > &traces, Array< OneD, NekDouble > &loctraces)
Definition LocTraceToTraceMap.cpp:1711
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:454
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:463
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:64
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:433
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:440
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:1299
Nektar::MultiRegions::LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_3D
void CalcLocTracePhysToTraceIDMap_3D(const ExpListSharedPtr &tracelist)
Definition LocTraceToTraceMap.cpp:1026
Nektar::MultiRegions::LocTraceToTraceMap::m_nLocTracePts
int m_nLocTracePts
The number of local trace points.
Definition LocTraceToTraceMap.h:403
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:1411
Nektar::MultiRegions::LocTraceToTraceMap::m_nTracePts
int m_nTracePts
The number of global trace points.
Definition LocTraceToTraceMap.h:405
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:1366
Nektar::MultiRegions::LocTraceToTraceMap::InterpLocTracesToTrace
void InterpLocTracesToTrace(const int dir, const Array< OneD, const NekDouble > &loctraces, Array< OneD, NekDouble > &traces)
Definition LocTraceToTraceMap.cpp:1380
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:460
Nektar::MultiRegions::LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap
void CalcLocTracePhysToTraceIDMap(const ExpListSharedPtr &tracelist, const int ndim)
Definition LocTraceToTraceMap.cpp:931
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:1317
Nektar::MultiRegions::LocTraceToTraceMap::m_nTraceCoeffs
int m_nTraceCoeffs[2]
Number of forwards/backwards trace coefficients.
Definition LocTraceToTraceMap.h:448
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:401
Nektar::MultiRegions::LocTraceToTraceMap::m_traceFieldMap
Array< OneD, TraceFieldMapEssential > m_traceFieldMap
Definition LocTraceToTraceMap.h:510
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:471
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:2156
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:444
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:1118
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:451
Nektar::MultiRegions::LocTraceToTraceMap::m_expdim
int m_expdim
Expansion Dimension we have setup for trace mapping.
Definition LocTraceToTraceMap.h:398
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:419
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:1334
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:416
Nektar::MultiRegions::LocTraceToTraceMap::InterpTraceToLocEdges
void InterpTraceToLocEdges(const int dir, const Array< OneD, const NekDouble > &edges, Array< OneD, NekDouble > &locedges)
Interpolate global trace edge to local trace edges point distributions where required.
Definition LocTraceToTraceMap.cpp:2060
Nektar::MultiRegions::LocTraceToTraceMap::m_ElemNeighbsId
Array< OneD, Array< OneD, int > > m_ElemNeighbsId
Definition LocTraceToTraceMap.h:482
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:436
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:409
Nektar::MultiRegions::LocTraceToTraceMap::m_interpTracePtsEntry
Array< OneD, Array< OneD, int > > m_interpTracePtsEntry
start entry of each global trace in m_locInterpTraceToTraceMap, referenced by element and trace id
Definition LocTraceToTraceMap.h:496
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:429
Nektar::MultiRegions::LocTraceToTraceMap::m_traceInterp
Array< OneD, TraceInterpEssential > m_traceInterp
Definition LocTraceToTraceMap.h:511
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:2325
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:1832
Nektar::MultiRegions::LocTraceToTraceMap::m_locTracePtsEntry
Array< OneD, Array< OneD, int > > m_locTracePtsEntry
start entry of each local trace in m_locTraceToFieldMap, referenced by element and trace id
Definition LocTraceToTraceMap.h:493
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:425
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:422
Nektar::MultiRegions::LocTraceToTraceMap::m_traceCoeffsEntry
Array< OneD, Array< OneD, int > > m_traceCoeffsEntry
start entry of each global trace in m_traceCoeffsToElmtMap, referenced by element and trace id
Definition LocTraceToTraceMap.h:501
Nektar::MultiRegions::LocTraceToTraceMap::GetTraceFieldMapEssential
const TraceFieldMapEssential & GetTraceFieldMapEssential(const int cid)
Definition LocTraceToTraceMap.cpp:925
Nektar::MultiRegions::LocTraceToTraceMap::FindElmtNeighbors
void FindElmtNeighbors(const ExpList &locExp, const ExpListSharedPtr &trace)
Definition LocTraceToTraceMap.cpp:1178
Nektar::MultiRegions::LocTraceToTraceMap::InterpTraceToLocTrace
void InterpTraceToLocTrace(const int dir, const Array< OneD, NekDouble > &traces, Array< OneD, NekDouble > &loctraces)
Definition LocTraceToTraceMap.cpp:2028
Nektar::MultiRegions::LocTraceToTraceMap::m_interpTraceIndex
Array< OneD, Array< OneD, int > > m_interpTraceIndex
subscript of m_interpTrace, referenced by element and trace id
Definition LocTraceToTraceMap.h:498
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:1348
Nektar::MultiRegions::LocTraceToTraceMap::CalcLocTracePhysToTraceIDMap_2D
void CalcLocTracePhysToTraceIDMap_2D(const ExpListSharedPtr &tracelist)
Definition LocTraceToTraceMap.cpp:948
Nektar::MultiRegions::LocTraceToTraceMap::m_LocTracephysToTraceIDMap
Array< OneD, Array< OneD, int > > m_LocTracephysToTraceIDMap
Definition LocTraceToTraceMap.h:484
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:1523
Nektar::MultiRegions::LocTraceToTraceMap::GetTraceInterpEssential
const TraceInterpEssential & GetTraceInterpEssential(const int cid)
Definition LocTraceToTraceMap.cpp:919
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:152
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:324
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::LibUtilities::PointsManager
PointsManagerT & PointsManager(void)
Definition ManagerAccess.cpp:40
Nektar::LibUtilities::eGaussLegendreWithMP
@ eGaussLegendreWithMP
1D Gauss-Legendre quadrature points with additional x=-1 and x=1 end points
Definition PointsType.h:95
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition PointsType.h:51
Nektar::LibUtilities::eGaussLegendreWithM
@ eGaussLegendreWithM
1D Gauss-Legendre quadrature points with additional x=-1 point
Definition PointsType.h:97
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition Expansion.h:66
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:1817
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition NekTypeDefs.hpp:75
Vmath::Gathr
void Gathr(I n, const T *x, const I *y, T *z)
Gather vector z[i] = x[y[i]].
Definition Vmath.hpp:507
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition Vmath.hpp:396
Vmath::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::Dot
T Dot(int n, const T *w, const T *x)
dot product
Definition Vmath.hpp:761
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::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition Vmath.hpp:273
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:295
Generated on Sat Jul 4 2026 01:02:22 for Nektar++ by doxygen 1.9.8