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