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