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 
36 #include <MultiRegions/ExpList.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,
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  int m_expdim = locExpVector[0]->GetShapeDimension();
74 
75  // Switch between 1D, 2D and 3D
76  switch (m_expdim)
77  {
78  case 1:
79  break;
80  case 2:
81  Setup2D(locExp, trace, elmtToTrace, LeftAdjacents);
82  break;
83  case 3:
84  Setup3D(locExp, trace, elmtToTrace, LeftAdjacents);
85  break;
86  default:
87  ASSERTL0(false, "Number of dimensions greater than 3")
88  break;
89  }
90 }
91 
92 LocTraceToTraceMap::~LocTraceToTraceMap()
93 {
94 }
95 
96 /**
97  * @brief Set up member variables for a two-dimensional problem.
98  *
99  * @param locExp Expansion list of 2D elements
100  * @param trace Expansion list of the one-dimensional trace.
101  * @param elmtToTrace Mapping from elemental edges to trace.
102  * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
103  */
104 void LocTraceToTraceMap::Setup2D(
105  const ExpList &locExp,
106  const ExpListSharedPtr &trace,
108  &elmtToTrace,
109  const vector<bool> &LeftAdjacents)
110 {
111  m_LocTraceToTraceMap = Array<OneD, Array<OneD, int> >(2);
113  m_interpTraceI0 = Array<OneD, Array<OneD, DNekMatSharedPtr> >(2);
114  m_interpEndPtI0 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(2);
115  m_interpPoints = Array<OneD, Array<OneD, TraceInterpPoints> >(2);
116  m_interpNfaces = Array<OneD, Array<OneD, int> >(2);
117 
118  m_traceCoeffsToElmtMap = Array<OneD, Array<OneD, int> >(2);
119  m_traceCoeffsToElmtTrace = Array<OneD, Array<OneD, int> >(2);
120  m_traceCoeffsToElmtSign = Array<OneD, Array<OneD, int> >(2);
121 
123  const std::shared_ptr<LocalRegions::ExpansionVector> exp =
124  locExp.GetExp();
125 
126  int cnt, n, e, phys_offset;
127 
128  int nexp = exp->size();
129  m_nTracePts = trace->GetTotPoints();
130 
131  // Count number of edges and points required for maps
132  int nFwdPts = 0;
133  int nBwdPts = 0;
134  int nFwdCoeffs = 0;
135  int nBwdCoeffs = 0;
136  m_nFwdLocTracePts = 0;
137  m_nLocTracePts = 0;
138 
139  for (cnt = n = 0; n < nexp; ++n)
140  {
141  exp2d = (*exp)[n]->as<LocalRegions::Expansion2D>();
142 
143  for (int i = 0; i < exp2d->GetNedges(); ++i, ++cnt)
144  {
145  int nLocPts = exp2d->GetEdgeNumPoints(i);
146  m_nLocTracePts += nLocPts;
147 
148  if (LeftAdjacents[cnt])
149  {
150  nFwdPts += elmtToTrace[n][i]->GetTotPoints();
151  nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
152  m_nFwdLocTracePts += nLocPts;
153  }
154  else
155  {
156  nBwdPts += elmtToTrace[n][i]->GetTotPoints();
157  nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
158  }
159  }
160  }
161 
162  m_fieldToLocTraceMap = Array<OneD, int>(m_nLocTracePts);
163 
164  m_LocTraceToTraceMap[0] = Array<OneD, int>(nFwdPts);
165  m_LocTraceToTraceMap[1] = Array<OneD, int>(nBwdPts);
166 
167  m_nTraceCoeffs[0] = nFwdCoeffs;
168  m_nTraceCoeffs[1] = nBwdCoeffs;
169 
170  m_traceCoeffsToElmtMap[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
171  m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
172  m_traceCoeffsToElmtTrace[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
173  m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
174  m_traceCoeffsToElmtSign[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
175  m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
176 
177  // Gather information about trace interpolations
178  map<TraceInterpPoints, vector<pair<int, int> >, cmpop> TraceInterpMap;
179 
180  vector<vector<int> > TraceOrder;
181  TraceOrder.resize(nexp);
182  int nedge;
183  int fwdcnt = 0;
184  int bwdcnt = 0;
185 
186  // Generate a map of similar traces with the same
187  // interpolation requirements
188  for (cnt = n = 0; n < nexp; ++n)
189  {
190  exp2d = (*exp)[n]->as<LocalRegions::Expansion2D>();
191  nedge = exp2d->GetNedges();
192  TraceOrder[n].resize(nedge);
193 
194  int coeffoffset = locExp.GetCoeff_Offset(n);
195  for (e = 0; e < nedge; ++e, ++cnt)
196  {
197  StdRegions::StdExpansionSharedPtr edge = elmtToTrace[n][e];
198  StdRegions::Orientation orient = exp2d->GetEorient(e);
199 
200  LibUtilities::PointsKey fromPointsKey0, fromPointsKey1;
201  LibUtilities::PointsKey toPointsKey0, toPointsKey1;
202 
203  // For Spencer's data structure
204  int basisDir = 0;
205  if (exp2d->DetShapeType() == LibUtilities::eTriangle)
206  {
207  basisDir = e == 0 ? 0 : 1;
208  }
209  else
210  {
211  basisDir = e % 2;
212  }
213 
214  fromPointsKey0 = exp2d->GetBasis(basisDir)->GetPointsKey();
215  fromPointsKey1 =
217  toPointsKey0 = edge->GetBasis(0)->GetPointsKey();
218  toPointsKey1 =
220 
221  // Spencer's data structure
222  TraceInterpPoints fpoint(
223  fromPointsKey0, fromPointsKey1, toPointsKey0, toPointsKey1);
224 
225  pair<int, int> epf(n, e);
226  TraceInterpMap[fpoint].push_back(epf);
227  TraceOrder[n][e] = cnt;
228 
229  // Setup for coefficient mapping from trace normal flux
230  // to elements
233  // Not sure about the call GetBasisNumModes passing (0)!
234  exp2d->GetEdgeToElementMap(
235  e, orient, map, sign, edge->GetBasisNumModes(0));
236 
237  int order_f = edge->GetNcoeffs();
238  int foffset = trace->GetCoeff_Offset(edge->GetElmtId());
239 
240  double fac = (*exp)[n]->EdgeNormalNegated(e) ? -1.0 : 1.0;
241 
243  elmtToTrace[n][e]->as<LocalRegions::Expansion1D>();
244 
245  if (exp2d->GetEdgeExp(e)->GetRightAdjacentElementExp())
246  {
247  if (locExp1d->GetRightAdjacentElementExp()
248  ->GetGeom()
249  ->GetGlobalID() == exp2d->GetGeom()->GetGlobalID())
250  {
251  fac = -1.0;
252  }
253  }
254 
255  if (LeftAdjacents[cnt])
256  {
257  for (int i = 0; i < order_f; ++i)
258  {
259  m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
260  m_traceCoeffsToElmtTrace[0][fwdcnt] = foffset + i;
261  m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
262  }
263  }
264  else
265  {
266  for (int i = 0; i < order_f; ++i)
267  {
268  m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
269  m_traceCoeffsToElmtTrace[1][bwdcnt] = foffset + i;
270  m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
271  }
272  }
273  }
274  }
275 
276  int nInterpType = TraceInterpMap.size();
277 
278  for (int i = 0; i < 2; ++i)
279  {
280  m_interpTrace[i] = Array<OneD, InterpLocTraceToTrace>(nInterpType);
281  m_interpTraceI0[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
282  m_interpEndPtI0[i] = Array<OneD, Array<OneD, NekDouble> >(nInterpType);
283  m_interpPoints[i] = Array<OneD, TraceInterpPoints>(nInterpType);
284  m_interpNfaces[i] = Array<OneD, int>(nInterpType, 0);
285  }
286 
287  int nedgepts, nedgepts1;
288  int cnt1 = 0;
289  int cnt2 = 0;
290  int cntFwd = 0;
291  int cntBwd = 0;
292  int cntFwd1 = 0;
293  int cntBwd1 = 0;
294  int set;
295  Array<OneD, int> edgeids;
296  Array<OneD, int> locTraceToTraceMap;
297  cnt = 0;
298 
299  for (auto it = TraceInterpMap.begin(); it != TraceInterpMap.end();
300  ++it, ++cnt1)
301  {
302  LibUtilities::PointsKey fromPointsKey0 = std::get<0>(it->first);
303  LibUtilities::PointsKey toPointsKey0 = std::get<2>(it->first);
304 
305  bool fwdSet = false;
306  bool bwdSet = false;
307 
308  for (int f = 0; f < it->second.size(); ++f, ++cnt2)
309  {
310  n = it->second[f].first;
311  e = it->second[f].second;
312 
313  StdRegions::StdExpansionSharedPtr edge = elmtToTrace[n][e];
314 
315  exp2d = (*exp)[n]->as<LocalRegions::Expansion2D>();
316  phys_offset = locExp.GetPhys_Offset(n);
317 
318  // Mapping of new edge order to one that loops over elmts
319  // then set up mapping of faces in standard cartesian order
320  exp2d->GetEdgePhysMap(e, edgeids);
321 
322  nedgepts = exp2d->GetEdgeNumPoints(e);
323  nedgepts1 = edge->GetTotPoints();
324 
325  StdRegions::Orientation orient = exp2d->GetEorient(e);
326 
327  // Account for eBackwards orientation
328  exp2d->ReOrientEdgePhysMap(elmtToTrace[n][e]->GetNverts(),
329  orient,
330  toPointsKey0.GetNumPoints(),
331  locTraceToTraceMap);
332 
333  int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
334 
335  if (LeftAdjacents[TraceOrder[n][e]])
336  {
337  for (int i = 0; i < nedgepts; ++i)
338  {
339  m_fieldToLocTraceMap[cntFwd + i] = phys_offset + edgeids[i];
340  }
341 
342  for (int i = 0; i < nedgepts1; ++i)
343  {
344  m_LocTraceToTraceMap[0][cntFwd1 + i] =
345  offset + locTraceToTraceMap[i];
346  }
347 
348  cntFwd += nedgepts;
349  cntFwd1 += nedgepts1;
350  set = 0;
351  }
352  else
353  {
354  for (int i = 0; i < nedgepts; ++i)
355  {
356  m_fieldToLocTraceMap[m_nFwdLocTracePts + cntBwd + i] =
357  phys_offset + edgeids[i];
358  }
359 
360  for (int i = 0; i < nedgepts1; ++i)
361  {
362  m_LocTraceToTraceMap[1][cntBwd1 + i] =
363  offset + locTraceToTraceMap[i];
364  }
365 
366  cntBwd += nedgepts;
367  cntBwd1 += nedgepts1;
368  set = 1;
369  }
370 
371  m_interpNfaces[set][cnt1] += 1;
372 
373  if ((fwdSet == false && set == 0) || (bwdSet == false && set == 1))
374  {
375  m_interpPoints[set][cnt1] = it->first;
376 
377  if (fromPointsKey0 == toPointsKey0)
378  {
379  m_interpTrace[set][cnt1] = eNoInterp;
380  }
381  else
382  {
383  m_interpTrace[set][cnt1] = eInterpDir0;
384  m_interpTraceI0[set][cnt1] =
385  LibUtilities::PointsManager()[fromPointsKey0]->GetI(
386  toPointsKey0);
387 
388  // Check to see if we can
389  // just interpolate endpoint
390  if ((fromPointsKey0.GetPointsType() ==
392  (toPointsKey0.GetPointsType() ==
394  {
395  if (fromPointsKey0.GetNumPoints() + 1 ==
396  toPointsKey0.GetNumPoints())
397  {
398  m_interpTrace[set][cnt1] = eInterpEndPtDir0;
399 
400  int fnp0 = fromPointsKey0.GetNumPoints();
401 
402  int tnp0 = toPointsKey0.GetNumPoints();
403 
404  m_interpEndPtI0[set][cnt1] =
406 
407  Vmath::Vcopy(
408  fnp0,
409  m_interpTraceI0[set][cnt1]->GetPtr().get() +
410  tnp0 - 1,
411  tnp0,
412  &m_interpEndPtI0[set][cnt1][0],
413  1);
414  }
415  }
416  }
417 
418  if (set == 0)
419  {
420  fwdSet = true;
421  }
422  else
423  {
424  bwdSet = true;
425  }
426  }
427  }
428  }
429 }
430 
431 /**
432  * @brief Set up member variables for a three-dimensional problem.
433  *
434  * @param locExp Expansion list of 3D elements
435  * @param trace Expansion list of the two-dimensional trace.
436  * @param elmtToTrace Mapping from elemental faces to trace.
437  * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
438  */
439 void LocTraceToTraceMap::Setup3D(
440  const ExpList &locExp,
441  const ExpListSharedPtr &trace,
443  &elmtToTrace,
444  const vector<bool> &LeftAdjacents)
445 {
446  m_LocTraceToTraceMap = Array<OneD, Array<OneD, int> >(2);
447 
449  m_interpTraceI0 = Array<OneD, Array<OneD, DNekMatSharedPtr> >(2);
450  m_interpTraceI1 = Array<OneD, Array<OneD, DNekMatSharedPtr> >(2);
451  m_interpEndPtI0 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(2);
452  m_interpEndPtI1 = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(2);
453  m_interpPoints = Array<OneD, Array<OneD, TraceInterpPoints> >(2);
454  m_interpNfaces = Array<OneD, Array<OneD, int> >(2);
455 
456  m_traceCoeffsToElmtMap = Array<OneD, Array<OneD, int> >(2);
457  m_traceCoeffsToElmtTrace = Array<OneD, Array<OneD, int> >(2);
458  m_traceCoeffsToElmtSign = Array<OneD, Array<OneD, int> >(2);
459 
461  const std::shared_ptr<LocalRegions::ExpansionVector> exp =
462  locExp.GetExp();
463 
464  int cnt, n, e, phys_offset;
465 
466  int nexp = exp->size();
467  m_nTracePts = trace->GetTotPoints();
468 
469  // Count number of faces and points required for maps
470  int nFwdPts = 0;
471  int nBwdPts = 0;
472  int nFwdCoeffs = 0;
473  int nBwdCoeffs = 0;
474  m_nFwdLocTracePts = 0;
475  m_nLocTracePts = 0;
476 
477  for (cnt = n = 0; n < nexp; ++n)
478  {
479  exp3d = (*exp)[n]->as<LocalRegions::Expansion3D>();
480 
481  for (int i = 0; i < exp3d->GetNfaces(); ++i, ++cnt)
482  {
483  int nLocPts = exp3d->GetFaceNumPoints(i);
484  m_nLocTracePts += nLocPts;
485 
486  if (LeftAdjacents[cnt])
487  {
488  nFwdPts += elmtToTrace[n][i]->GetTotPoints();
489  nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
490  m_nFwdLocTracePts += nLocPts;
491  }
492  else
493  {
494  nBwdPts += elmtToTrace[n][i]->GetTotPoints();
495  nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
496  }
497  }
498  }
499 
500  m_fieldToLocTraceMap = Array<OneD, int>(m_nLocTracePts);
501 
502  m_LocTraceToTraceMap[0] = Array<OneD, int>(nFwdPts);
503  m_LocTraceToTraceMap[1] = Array<OneD, int>(nBwdPts);
504 
505  m_nTraceCoeffs[0] = nFwdCoeffs;
506  m_nTraceCoeffs[1] = nBwdCoeffs;
507 
508  m_traceCoeffsToElmtMap[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
509  m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
510  m_traceCoeffsToElmtTrace[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
511  m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
512  m_traceCoeffsToElmtSign[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
513  m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
514 
515  // Gather information about trace interpolations
516  map<TraceInterpPoints, vector<pair<int, int> >, cmpop> TraceInterpMap;
517 
518  vector<vector<int> > TraceOrder;
519  TraceOrder.resize(nexp);
520  int nface;
521  int fwdcnt = 0;
522  int bwdcnt = 0;
523 
524  // Generate a map of similar traces with the same
525  // interpolation or projection requirements
526  for (cnt = n = 0; n < nexp; ++n)
527  {
528  exp3d = (*exp)[n]->as<LocalRegions::Expansion3D>();
529  nface = exp3d->GetNfaces();
530  TraceOrder[n].resize(nface);
531 
532  int coeffoffset = locExp.GetCoeff_Offset(n);
533  for (e = 0; e < nface; ++e, ++cnt)
534  {
535  StdRegions::StdExpansionSharedPtr face = elmtToTrace[n][e];
536  StdRegions::Orientation orient = exp3d->GetForient(e);
537 
538  LibUtilities::PointsKey fromPointsKey0, fromPointsKey1;
539  LibUtilities::PointsKey toPointsKey0, toPointsKey1;
540 
541  // 3D specific
542  int dir0 = exp3d->GetGeom3D()->GetDir(e, 0);
543  int dir1 = exp3d->GetGeom3D()->GetDir(e, 1);
544 
545  fromPointsKey0 = exp3d->GetBasis(dir0)->GetPointsKey();
546  fromPointsKey1 = exp3d->GetBasis(dir1)->GetPointsKey();
547 
549  {
550  toPointsKey0 = face->GetBasis(0)->GetPointsKey();
551  toPointsKey1 = face->GetBasis(1)->GetPointsKey();
552  }
553  else // transpose points key evaluation
554  {
555  toPointsKey0 = face->GetBasis(1)->GetPointsKey();
556  toPointsKey1 = face->GetBasis(0)->GetPointsKey();
557  }
558 
559  TraceInterpPoints fpoint(
560  fromPointsKey0, fromPointsKey1, toPointsKey0, toPointsKey1);
561 
562  pair<int, int> epf(n, e);
563  TraceInterpMap[fpoint].push_back(epf);
564  TraceOrder[n][e] = cnt;
565 
566  // Setup for coefficient mapping from trace normal
567  // flux to elements
570  exp3d->GetFaceToElementMap(e,
571  orient,
572  map,
573  sign,
574  face->GetBasisNumModes(0),
575  face->GetBasisNumModes(1));
576 
577  int order_f = face->GetNcoeffs();
578  int foffset = trace->GetCoeff_Offset(face->GetElmtId());
579 
580  int fac = (*exp)[n]->FaceNormalNegated(e) ? -1.0 : 1.0;
581 
582  if (exp3d->GetFaceExp(e)->GetRightAdjacentElementExp())
583  {
584  if (exp3d->GetFaceExp(e)
585  ->GetRightAdjacentElementExp()
586  ->GetGeom3D()
587  ->GetGlobalID() == exp3d->GetGeom3D()->GetGlobalID())
588  {
589  fac = -1.0;
590  }
591  }
592 
593  if (LeftAdjacents[cnt])
594  {
595  for (int i = 0; i < order_f; ++i)
596  {
597  m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
598  m_traceCoeffsToElmtTrace[0][fwdcnt] = foffset + i;
599  m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
600  }
601  }
602  else
603  {
604  for (int i = 0; i < order_f; ++i)
605  {
606  m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
607  m_traceCoeffsToElmtTrace[1][bwdcnt] = foffset + i;
608  m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
609  }
610  }
611  }
612  }
613 
614  int nInterpType = TraceInterpMap.size();
615  for (int i = 0; i < 2; ++i)
616  {
617  m_interpTrace[i] = Array<OneD, InterpLocTraceToTrace>(nInterpType);
618  m_interpTraceI0[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
619  m_interpTraceI1[i] = Array<OneD, DNekMatSharedPtr>(nInterpType);
620  m_interpEndPtI0[i] = Array<OneD, Array<OneD, NekDouble> >(nInterpType);
621  m_interpEndPtI1[i] = Array<OneD, Array<OneD, NekDouble> >(nInterpType);
622  m_interpPoints[i] = Array<OneD, TraceInterpPoints>(nInterpType);
623  m_interpNfaces[i] = Array<OneD, int>(nInterpType, 0);
624  }
625 
626  int nfacepts, nfacepts1;
627  int cnt1 = 0;
628  int cnt2 = 0;
629  int cntFwd = 0;
630  int cntBwd = 0;
631  int cntFwd1 = 0;
632  int cntBwd1 = 0;
633  int set;
634  Array<OneD, int> faceids;
635  Array<OneD, int> locTraceToTraceMap;
636  cnt = 0;
637  for (auto it = TraceInterpMap.begin(); it != TraceInterpMap.end();
638  ++it, ++cnt1)
639  {
640  LibUtilities::PointsKey fromPointsKey0 = std::get<0>(it->first);
641  LibUtilities::PointsKey fromPointsKey1 = std::get<1>(it->first);
642  LibUtilities::PointsKey toPointsKey0 = std::get<2>(it->first);
643  LibUtilities::PointsKey toPointsKey1 = std::get<3>(it->first);
644 
645  bool fwdSet = false;
646  bool bwdSet = false;
647 
648  for (int f = 0; f < it->second.size(); ++f, ++cnt2)
649  {
650  n = it->second[f].first;
651  e = it->second[f].second;
652 
653  StdRegions::StdExpansionSharedPtr face = elmtToTrace[n][e];
654 
655  exp3d = (*exp)[n]->as<LocalRegions::Expansion3D>();
656  phys_offset = locExp.GetPhys_Offset(n);
657 
658  // mapping of new face order to one that loops over elmts
659  // then faces set up mapping of faces in standard cartesian
660  // order
661  exp3d->GetFacePhysMap(e, faceids);
662  nfacepts = exp3d->GetFaceNumPoints(e);
663  nfacepts1 = face->GetTotPoints();
664 
665  StdRegions::Orientation orient = exp3d->GetForient(e);
666 
667  exp3d->ReOrientFacePhysMap(elmtToTrace[n][e]->GetNverts(),
668  orient,
669  toPointsKey0.GetNumPoints(),
670  toPointsKey1.GetNumPoints(),
671  locTraceToTraceMap);
672 
673  int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
674 
675  if (LeftAdjacents[TraceOrder[n][e]])
676  {
677  for (int i = 0; i < nfacepts; ++i)
678  {
679  m_fieldToLocTraceMap[cntFwd + i] = phys_offset + faceids[i];
680  }
681 
682  for (int i = 0; i < nfacepts1; ++i)
683  {
684  m_LocTraceToTraceMap[0][cntFwd1 + i] =
685  offset + locTraceToTraceMap[i];
686  }
687 
688  cntFwd += nfacepts;
689  cntFwd1 += nfacepts1;
690  set = 0;
691  }
692  else
693  {
694  for (int i = 0; i < nfacepts; ++i)
695  {
696  m_fieldToLocTraceMap[m_nFwdLocTracePts + cntBwd + i] =
697  phys_offset + faceids[i];
698  }
699 
700  for (int i = 0; i < nfacepts1; ++i)
701  {
702  m_LocTraceToTraceMap[1][cntBwd1 + i] =
703  offset + locTraceToTraceMap[i];
704  }
705 
706  cntBwd += nfacepts;
707  cntBwd1 += nfacepts1;
708  set = 1;
709  }
710 
711  m_interpNfaces[set][cnt1] += 1;
712 
713  if (((fwdSet == false) && (set == 0)) ||
714  ((bwdSet == false) && (set == 1)))
715  {
716  m_interpPoints[set][cnt1] = it->first;
717 
718  if (fromPointsKey0 == toPointsKey0)
719  {
720  if (fromPointsKey1 == toPointsKey1)
721  {
722  m_interpTrace[set][cnt1] = eNoInterp;
723  }
724  else
725  {
726  m_interpTrace[set][cnt1] = eInterpDir1;
727  m_interpTraceI1[set][cnt1] =
728  LibUtilities::PointsManager()[fromPointsKey1]->GetI(
729  toPointsKey1);
730 
731  // Check to see if we can just
732  // interpolate endpoint
733  if ((fromPointsKey1.GetPointsType() ==
735  (toPointsKey1.GetPointsType() ==
737  {
738  if (fromPointsKey1.GetNumPoints() + 1 ==
739  toPointsKey1.GetNumPoints())
740  {
741  m_interpTrace[set][cnt1] = eInterpEndPtDir1;
742  int fnp1 = fromPointsKey1.GetNumPoints();
743  int tnp1 = toPointsKey1.GetNumPoints();
744  m_interpEndPtI1[set][cnt1] =
746  Vmath::Vcopy(
747  fnp1,
748  m_interpTraceI1[set][cnt1]->GetPtr().get() +
749  tnp1 - 1,
750  tnp1,
751  &m_interpEndPtI1[set][cnt1][0],
752  1);
753  }
754  }
755  }
756  }
757  else
758  {
759  if (fromPointsKey1 == toPointsKey1)
760  {
761  m_interpTrace[set][cnt1] = eInterpDir0;
762  m_interpTraceI0[set][cnt1] =
763  LibUtilities::PointsManager()[fromPointsKey0]->GetI(
764  toPointsKey0);
765 
766  // Check to see if we can just
767  // interpolate endpoint
768  if ((fromPointsKey0.GetPointsType() ==
770  (toPointsKey0.GetPointsType() ==
772  {
773  if (fromPointsKey0.GetNumPoints() + 1 ==
774  toPointsKey0.GetNumPoints())
775  {
776  m_interpTrace[set][cnt1] = eInterpEndPtDir0;
777  int fnp0 = fromPointsKey0.GetNumPoints();
778  int tnp0 = toPointsKey0.GetNumPoints();
779  m_interpEndPtI0[set][cnt1] =
781  Vmath::Vcopy(
782  fnp0,
783  m_interpTraceI0[set][cnt1]->GetPtr().get() +
784  tnp0 - 1,
785  tnp0,
786  &m_interpEndPtI0[set][cnt1][0],
787  1);
788  }
789  }
790  }
791  else
792  {
793  m_interpTrace[set][cnt1] = eInterpBothDirs;
794  m_interpTraceI0[set][cnt1] =
795  LibUtilities::PointsManager()[fromPointsKey0]->GetI(
796  toPointsKey0);
797  m_interpTraceI1[set][cnt1] =
798  LibUtilities::PointsManager()[fromPointsKey1]->GetI(
799  toPointsKey1);
800 
801  // check to see if we can just
802  // interpolate endpoint
803  if ((fromPointsKey0.GetPointsType() ==
805  (toPointsKey0.GetPointsType() ==
807  {
808  if (fromPointsKey0.GetNumPoints() + 1 ==
809  toPointsKey0.GetNumPoints())
810  {
811  m_interpTrace[set][cnt1] =
813  int fnp0 = fromPointsKey0.GetNumPoints();
814  int tnp0 = toPointsKey0.GetNumPoints();
815  m_interpEndPtI0[set][cnt1] =
817  Vmath::Vcopy(
818  fnp0,
819  m_interpTraceI0[set][cnt1]->GetPtr().get() +
820  tnp0 - 1,
821  tnp0,
822  &m_interpEndPtI0[set][cnt1][0],
823  1);
824  }
825  }
826  }
827  }
828 
829  if (set == 0)
830  {
831  fwdSet = true;
832  }
833  else
834  {
835  bwdSet = true;
836  }
837  }
838  }
839  }
840 }
841 
842 /**
843  * @brief Gather the local traces in physical space from field using
844  * #m_fieldToLocTraceMap.
845  *
846  * @param field Solution field in physical space
847  * @param faces Resulting local traces.
848  */
849 void LocTraceToTraceMap::LocTracesFromField(
851 {
852  Vmath::Gathr(m_fieldToLocTraceMap.num_elements(),
853  field,
854  m_fieldToLocTraceMap,
855  faces);
856 }
857 
858 /**
859  * @brief Gather the forwards-oriented local traces in physical space from field
860  * using #m_fieldToLocTraceMap.
861  *
862  * @param field Solution field in physical space
863  * @param faces Resulting local forwards-oriented traces.
864  */
865 void LocTraceToTraceMap::FwdLocTracesFromField(
867 {
868  Vmath::Gathr(m_nFwdLocTracePts, field, m_fieldToLocTraceMap, faces);
869 }
870 
871 /**
872  * @brief Interpolate local trace edges to global trace edge point distributions
873  * where required.
874  *
875  * @param dir Selects forwards (0) or backwards (1) direction.
876  * @param locfaces Local trace edge storage.
877  * @param faces Global trace edge storage
878  */
879 void LocTraceToTraceMap::InterpLocEdgesToTrace(
880  const int dir,
881  const Array<OneD, const NekDouble> &locedges,
883 {
884  ASSERTL1(dir < 2,
885  "option dir out of range, "
886  " dir=0 is fwd, dir=1 is bwd");
887 
888  int cnt = 0;
889  int cnt1 = 0;
890 
891  // tmp space assuming forward map is of size of trace
892  Array<OneD, NekDouble> tmp(m_nTracePts);
893 
894  for (int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
895  {
896  // Check if there are edges to interpolate
897  if (m_interpNfaces[dir][i])
898  {
899  // Get to/from points
900  LibUtilities::PointsKey fromPointsKey0 =
901  std::get<0>(m_interpPoints[dir][i]);
902  LibUtilities::PointsKey toPointsKey0 =
903  std::get<2>(m_interpPoints[dir][i]);
904 
905  int fnp = fromPointsKey0.GetNumPoints();
906  int tnp = toPointsKey0.GetNumPoints();
907  int nedges = m_interpNfaces[dir][i];
908 
909  // Do interpolation here if required
910  switch (m_interpTrace[dir][i])
911  {
912  case eNoInterp: // Just copy
913  {
914  Vmath::Vcopy(nedges * fnp,
915  locedges.get() + cnt,
916  1,
917  tmp.get() + cnt1,
918  1);
919  }
920  break;
921  case eInterpDir0:
922  {
923  DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
924  Blas::Dgemm('N',
925  'N',
926  tnp,
927  nedges,
928  fnp,
929  1.0,
930  I0->GetPtr().get(),
931  tnp,
932  locedges.get() + cnt,
933  fnp,
934  0.0,
935  tmp.get() + cnt1,
936  tnp);
937  }
938  break;
939  case eInterpEndPtDir0:
940  {
941  Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
942 
943  for (int k = 0; k < nedges; ++k)
944  {
945  Vmath::Vcopy(fnp,
946  &locedges[cnt + k * fnp],
947  1,
948  &tmp[cnt1 + k * tnp],
949  1);
950 
951  tmp[cnt1 + k * tnp + tnp - 1] = Blas::Ddot(
952  fnp, locedges.get() + cnt + k * fnp, 1, &I0[0], 1);
953  }
954  }
955  break;
956  default:
957  ASSERTL0(false,
958  "Invalid interpolation type for 2D elements");
959  break;
960  }
961 
962  cnt += nedges * fnp;
963  cnt1 += nedges * tnp;
964  }
965  }
966 
967  Vmath::Scatr(m_LocTraceToTraceMap[dir].num_elements(),
968  tmp.get(),
969  m_LocTraceToTraceMap[dir].get(),
970  edges.get());
971 }
972 
973 /**
974  * @brief Interpolate local faces to trace face point distributions where
975  * required.
976  *
977  * @param dir Selects forwards (0) or backwards (1) direction.
978  * @param locfaces Local trace face storage.
979  * @param faces Global trace face storage
980  */
981 void LocTraceToTraceMap::InterpLocFacesToTrace(
982  const int dir,
983  const Array<OneD, const NekDouble> &locfaces,
985 {
986  ASSERTL1(dir < 2,
987  "option dir out of range, "
988  " dir=0 is fwd, dir=1 is bwd");
989 
990  int cnt1 = 0;
991  int cnt = 0;
992 
993  // tmp space assuming forward map is of size of trace
994  Array<OneD, NekDouble> tmp(m_nTracePts);
995 
996  for (int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
997  {
998  // Check if there are faces to interpolate
999  if (m_interpNfaces[dir][i])
1000  {
1001  // Get to/from points
1002  LibUtilities::PointsKey fromPointsKey0 =
1003  std::get<0>(m_interpPoints[dir][i]);
1004  LibUtilities::PointsKey fromPointsKey1 =
1005  std::get<1>(m_interpPoints[dir][i]);
1006  LibUtilities::PointsKey toPointsKey0 =
1007  std::get<2>(m_interpPoints[dir][i]);
1008  LibUtilities::PointsKey toPointsKey1 =
1009  std::get<3>(m_interpPoints[dir][i]);
1010 
1011  int fnp0 = fromPointsKey0.GetNumPoints();
1012  int fnp1 = fromPointsKey1.GetNumPoints();
1013  int tnp0 = toPointsKey0.GetNumPoints();
1014  int tnp1 = toPointsKey1.GetNumPoints();
1015  int nfromfacepts = m_interpNfaces[dir][i] * fnp0 * fnp1;
1016  ;
1017 
1018  // Do interpolation here if required
1019  switch (m_interpTrace[dir][i])
1020  {
1021  case eNoInterp: // Just copy
1022  {
1023  Vmath::Vcopy(nfromfacepts,
1024  locfaces.get() + cnt,
1025  1,
1026  tmp.get() + cnt1,
1027  1);
1028  }
1029  break;
1030  case eInterpDir0:
1031  {
1032  DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1033  Blas::Dgemm('N',
1034  'N',
1035  tnp0,
1036  tnp1,
1037  fnp0,
1038  1.0,
1039  I0->GetPtr().get(),
1040  tnp0,
1041  locfaces.get() + cnt,
1042  fnp0,
1043  0.0,
1044  tmp.get() + cnt1,
1045  tnp0);
1046  }
1047  break;
1048  case eInterpEndPtDir0:
1049  {
1050  int nfaces = m_interpNfaces[dir][i];
1051  for (int k = 0; k < fnp0; ++k)
1052  {
1053  Vmath::Vcopy(nfaces * fnp1,
1054  locfaces.get() + cnt + k,
1055  fnp0,
1056  tmp.get() + cnt1 + k,
1057  tnp0);
1058  }
1059  Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1060  Blas::Dgemv('T',
1061  fnp0,
1062  tnp1 * m_interpNfaces[dir][i],
1063  1.0,
1064  tmp.get() + cnt1,
1065  tnp0,
1066  I0.get(),
1067  1,
1068  0.0,
1069  tmp.get() + cnt1 + tnp0 - 1,
1070  tnp0);
1071  }
1072  break;
1073  case eInterpDir1:
1074  {
1075  DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1076  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1077  {
1078  Blas::Dgemm('N',
1079  'T',
1080  tnp0,
1081  tnp1,
1082  fnp1,
1083  1.0,
1084  locfaces.get() + cnt + j * fnp0 * fnp1,
1085  tnp0,
1086  I1->GetPtr().get(),
1087  tnp1,
1088  0.0,
1089  tmp.get() + cnt1 + j * tnp0 * tnp1,
1090  tnp0);
1091  }
1092  }
1093  break;
1094  case eInterpEndPtDir1:
1095  {
1096  Array<OneD, NekDouble> I1 = m_interpEndPtI1[dir][i];
1097  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1098  {
1099  // copy all points
1100  Vmath::Vcopy(fnp0 * fnp1,
1101  locfaces.get() + cnt + j * fnp0 * fnp1,
1102  1,
1103  tmp.get() + cnt1 + j * tnp0 * tnp1,
1104  1);
1105 
1106  // interpolate end points
1107  for (int k = 0; k < tnp0; ++k)
1108  {
1109  tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1110  Blas::Ddot(fnp1,
1111  locfaces.get() + cnt +
1112  j * fnp0 * fnp1 + k,
1113  fnp0,
1114  &I1[0],
1115  1);
1116  }
1117  }
1118  }
1119  break;
1120  case eInterpBothDirs:
1121  {
1122  DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1123  DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1124  Array<OneD, NekDouble> wsp(m_interpNfaces[dir][i] * fnp0 *
1125  tnp1 * fnp0);
1126 
1127  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1128  {
1129  Blas::Dgemm('N',
1130  'T',
1131  fnp0,
1132  tnp1,
1133  fnp1,
1134  1.0,
1135  locfaces.get() + cnt + j * fnp0 * fnp1,
1136  fnp0,
1137  I1->GetPtr().get(),
1138  tnp1,
1139  0.0,
1140  wsp.get() + j * fnp0 * tnp1,
1141  fnp0);
1142  }
1143  Blas::Dgemm('N',
1144  'N',
1145  tnp0,
1146  tnp1 * m_interpNfaces[dir][i],
1147  fnp0,
1148  1.0,
1149  I0->GetPtr().get(),
1150  tnp0,
1151  wsp.get(),
1152  fnp0,
1153  0.0,
1154  tmp.get() + cnt1,
1155  tnp0);
1156  }
1157  break;
1159  {
1160  DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1161 
1162  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1163  {
1164  Blas::Dgemm('N',
1165  'T',
1166  fnp0,
1167  tnp1,
1168  fnp1,
1169  1.0,
1170  locfaces.get() + cnt + j * fnp0 * fnp1,
1171  fnp0,
1172  I1->GetPtr().get(),
1173  tnp1,
1174  0.0,
1175  tmp.get() + cnt1 + j * tnp0 * tnp1,
1176  tnp0);
1177  }
1178 
1179  Array<OneD, NekDouble> I0 = m_interpEndPtI0[dir][i];
1180  Blas::Dgemv('T',
1181  fnp0,
1182  tnp1 * m_interpNfaces[dir][i],
1183  1.0,
1184  tmp.get() + cnt1,
1185  tnp0,
1186  I0.get(),
1187  1,
1188  0.0,
1189  tmp.get() + cnt1 + tnp0 - 1,
1190  tnp0);
1191  }
1192  break;
1193  }
1194  cnt += nfromfacepts;
1195  cnt1 += m_interpNfaces[dir][i] * tnp0 * tnp1;
1196  }
1197  }
1198 
1199  Vmath::Scatr(m_LocTraceToTraceMap[dir].num_elements(),
1200  tmp.get(),
1201  m_LocTraceToTraceMap[dir].get(),
1202  faces.get());
1203 }
1204 
1205 /**
1206  * @brief Add contributions from trace coefficients to the elemental field
1207  * storage.
1208  *
1209  * @param trace Array of global trace coefficients.
1210  * @param field Array containing field coefficients storage.
1211  */
1212 void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1214 {
1215  int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
1216  for (int i = 0; i < nvals; ++i)
1217  {
1218  field[m_traceCoeffsToElmtMap[0][i]] +=
1219  m_traceCoeffsToElmtSign[0][i] *
1220  trace[m_traceCoeffsToElmtTrace[0][i]];
1221  }
1222 }
1223 
1224 /**
1225  * @brief Add contributions from backwards or forwards oriented trace
1226  * coefficients to the elemental field storage.
1227  *
1228  * @param dir Selects forwards (0) or backwards (1) direction
1229  * @param trace Array of global trace coefficients.
1230  * @param field Array containing field coefficients storage.
1231  */
1232 void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1233  const int dir,
1234  const Array<OneD, const NekDouble> &trace,
1235  Array<OneD, NekDouble> &field)
1236 {
1237  int nvals = m_nTraceCoeffs[dir];
1238  for (int i = 0; i < nvals; ++i)
1239  {
1240  field[m_traceCoeffsToElmtMap[dir][i]] +=
1241  m_traceCoeffsToElmtSign[dir][i] *
1242  trace[m_traceCoeffsToElmtTrace[dir][i]];
1243  }
1244 }
1245 
1246 }
1247 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
PointsType GetPointsType() const
Definition: Points.h:112
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:647
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:16
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
STL namespace.
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
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 A[m x n], B[n x k], C[m x k].
Definition: Blas.hpp:213
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2191
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:67
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:49
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:103
PointsManagerT & PointsManager(void)
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:676
Defines a specification for a set of points.
Definition: Points.h:59
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:168
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:140
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 ...
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
unsigned int GetNumPoints() const
Definition: Points.h:107
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:2200
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
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:2208
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51