Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Local trace to general trace mapping information
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <MultiRegions/ExpList.h>
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  */
63  const ExpList &locExp,
64  const ExpListSharedPtr &trace,
66  &elmtToTrace,
67  const vector<bool> &LeftAdjacents)
68 {
69  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
70 
71  // Assume that all the elements have same dimension
72  int m_expdim = locExpVector[0]->GetShapeDimension();
73 
74  // Switch between 1D, 2D and 3D
75  switch (m_expdim)
76  {
77  case 1:
78  break;
79  case 2:
80  Setup2D(locExp, trace, elmtToTrace, LeftAdjacents);
81  break;
82  case 3:
83  Setup3D(locExp, trace, elmtToTrace, LeftAdjacents);
84  break;
85  default:
86  ASSERTL0(false, "Number of dimensions greater than 3")
87  break;
88  }
89 }
90 
92 {
93 }
94 
95 /**
96  * @brief Set up member variables for a two-dimensional problem.
97  *
98  * @param locExp Expansion list of 2D elements
99  * @param trace Expansion list of the one-dimensional trace.
100  * @param elmtToTrace Mapping from elemental edges to trace.
101  * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
102  */
104  const ExpList &locExp,
105  const ExpListSharedPtr &trace,
107  &elmtToTrace,
108  const vector<bool> &LeftAdjacents)
109 {
116 
120 
122  const boost::shared_ptr<LocalRegions::ExpansionVector> exp =
123  locExp.GetExp();
124 
125  int cnt, n, e, phys_offset;
126 
127  int nexp = exp->size();
128  m_nTracePts = trace->GetTotPoints();
129 
130  // Count number of edges and points required for maps
131  int nFwdPts = 0;
132  int nBwdPts = 0;
133  int nFwdCoeffs = 0;
134  int nBwdCoeffs = 0;
135  m_nFwdLocTracePts = 0;
136  m_nLocTracePts = 0;
137 
138  for (cnt = n = 0; n < nexp; ++n)
139  {
140  exp2d = (*exp)[n]->as<LocalRegions::Expansion2D>();
141 
142  for (int i = 0; i < exp2d->GetNedges(); ++i, ++cnt)
143  {
144  int nLocPts = exp2d->GetEdgeNumPoints(i);
145  m_nLocTracePts += nLocPts;
146 
147  if (LeftAdjacents[cnt])
148  {
149  nFwdPts += elmtToTrace[n][i]->GetTotPoints();
150  nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
151  m_nFwdLocTracePts += nLocPts;
152  }
153  else
154  {
155  nBwdPts += elmtToTrace[n][i]->GetTotPoints();
156  nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
157  }
158  }
159  }
160 
162 
164  m_LocTraceToTraceMap[1] = Array<OneD, int>(nBwdPts);
165 
166  m_nTraceCoeffs[0] = nFwdCoeffs;
167  m_nTraceCoeffs[1] = nBwdCoeffs;
168 
169  m_traceCoeffsToElmtMap[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
170  m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
171  m_traceCoeffsToElmtTrace[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
172  m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
173  m_traceCoeffsToElmtSign[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
174  m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
175 
176  // Gather information about trace interpolations
177  map<TraceInterpPoints, vector<pair<int, int> >, cmpop> TraceInterpMap;
178  map<TraceInterpPoints, vector<pair<int, int> >, cmpop>::iterator it;
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 (it = TraceInterpMap.begin(); it != TraceInterpMap.end(); ++it, ++cnt1)
300  {
301  LibUtilities::PointsKey fromPointsKey0 = it->first.get<0>();
302  LibUtilities::PointsKey toPointsKey0 = it->first.get<2>();
303 
304  bool fwdSet = false;
305  bool bwdSet = false;
306 
307  for (int f = 0; f < it->second.size(); ++f, ++cnt2)
308  {
309  n = it->second[f].first;
310  e = it->second[f].second;
311 
312  StdRegions::StdExpansionSharedPtr edge = elmtToTrace[n][e];
313 
314  exp2d = (*exp)[n]->as<LocalRegions::Expansion2D>();
315  phys_offset = locExp.GetPhys_Offset(n);
316 
317  // Mapping of new edge order to one that loops over elmts
318  // then set up mapping of faces in standard cartesian order
319  exp2d->GetEdgePhysMap(e, edgeids);
320 
321  nedgepts = exp2d->GetEdgeNumPoints(e);
322  nedgepts1 = edge->GetTotPoints();
323 
324  StdRegions::Orientation orient = exp2d->GetCartesianEorient(e);
325 
326  // Account for eBackwards orientation
327  exp2d->ReOrientEdgePhysMap(elmtToTrace[n][e]->GetNverts(),
328  orient,
329  toPointsKey0.GetNumPoints(),
330  locTraceToTraceMap);
331 
332  int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
333 
334  if (LeftAdjacents[TraceOrder[n][e]])
335  {
336  for (int i = 0; i < nedgepts; ++i)
337  {
338  m_fieldToLocTraceMap[cntFwd + i] = phys_offset + edgeids[i];
339  }
340 
341  for (int i = 0; i < nedgepts1; ++i)
342  {
343  m_LocTraceToTraceMap[0][cntFwd1 + i] =
344  offset + locTraceToTraceMap[i];
345  }
346 
347  cntFwd += nedgepts;
348  cntFwd1 += nedgepts1;
349  set = 0;
350  }
351  else
352  {
353  for (int i = 0; i < nedgepts; ++i)
354  {
356  phys_offset + edgeids[i];
357  }
358 
359  for (int i = 0; i < nedgepts1; ++i)
360  {
361  m_LocTraceToTraceMap[1][cntBwd1 + i] =
362  offset + locTraceToTraceMap[i];
363  }
364 
365  cntBwd += nedgepts;
366  cntBwd1 += nedgepts1;
367  set = 1;
368  }
369 
370  m_interpNfaces[set][cnt1] += 1;
371 
372  if ((fwdSet == false && set == 0) || (bwdSet == false && set == 1))
373  {
374  m_interpPoints[set][cnt1] = it->first;
375 
376  if (fromPointsKey0 == toPointsKey0)
377  {
378  m_interpTrace[set][cnt1] = eNoInterp;
379  }
380  else
381  {
382  m_interpTrace[set][cnt1] = eInterpDir0;
383  m_interpTraceI0[set][cnt1] =
384  LibUtilities::PointsManager()[fromPointsKey0]->GetI(
385  toPointsKey0);
386 
387  // Check to see if we can
388  // just interpolate endpoint
389  if ((fromPointsKey0.GetPointsType() ==
391  (toPointsKey0.GetPointsType() ==
393  {
394  if (fromPointsKey0.GetNumPoints() + 1 ==
395  toPointsKey0.GetNumPoints())
396  {
397  m_interpTrace[set][cnt1] = eInterpEndPtDir0;
398 
399  int fnp0 = fromPointsKey0.GetNumPoints();
400 
401  int tnp0 = toPointsKey0.GetNumPoints();
402 
403  m_interpEndPtI0[set][cnt1] =
405 
406  Vmath::Vcopy(
407  fnp0,
408  m_interpTraceI0[set][cnt1]->GetPtr().get() +
409  tnp0 - 1,
410  tnp0,
411  &m_interpEndPtI0[set][cnt1][0],
412  1);
413  }
414  }
415  }
416 
417  if (set == 0)
418  {
419  fwdSet = true;
420  }
421  else
422  {
423  bwdSet = true;
424  }
425  }
426  }
427  }
428 }
429 
430 /**
431  * @brief Set up member variables for a three-dimensional problem.
432  *
433  * @param locExp Expansion list of 3D elements
434  * @param trace Expansion list of the two-dimensional trace.
435  * @param elmtToTrace Mapping from elemental faces to trace.
436  * @param leftAdjacents Vector of bools denoting forwards-oriented traces.
437  */
439  const ExpList &locExp,
440  const ExpListSharedPtr &trace,
442  &elmtToTrace,
443  const vector<bool> &LeftAdjacents)
444 {
446 
454 
458 
460  const boost::shared_ptr<LocalRegions::ExpansionVector> exp =
461  locExp.GetExp();
462 
463  int cnt, n, e, phys_offset;
464 
465  int nexp = exp->size();
466  m_nTracePts = trace->GetTotPoints();
467 
468  // Count number of faces and points required for maps
469  int nFwdPts = 0;
470  int nBwdPts = 0;
471  int nFwdCoeffs = 0;
472  int nBwdCoeffs = 0;
473  m_nFwdLocTracePts = 0;
474  m_nLocTracePts = 0;
475 
476  for (cnt = n = 0; n < nexp; ++n)
477  {
478  exp3d = (*exp)[n]->as<LocalRegions::Expansion3D>();
479 
480  for (int i = 0; i < exp3d->GetNfaces(); ++i, ++cnt)
481  {
482  int nLocPts = exp3d->GetFaceNumPoints(i);
483  m_nLocTracePts += nLocPts;
484 
485  if (LeftAdjacents[cnt])
486  {
487  nFwdPts += elmtToTrace[n][i]->GetTotPoints();
488  nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
489  m_nFwdLocTracePts += nLocPts;
490  }
491  else
492  {
493  nBwdPts += elmtToTrace[n][i]->GetTotPoints();
494  nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
495  }
496  }
497  }
498 
500 
502  m_LocTraceToTraceMap[1] = Array<OneD, int>(nBwdPts);
503 
504  m_nTraceCoeffs[0] = nFwdCoeffs;
505  m_nTraceCoeffs[1] = nBwdCoeffs;
506 
507  m_traceCoeffsToElmtMap[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
508  m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
509  m_traceCoeffsToElmtTrace[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
510  m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
511  m_traceCoeffsToElmtSign[0] = Array<OneD, int>(nFwdCoeffs + nBwdCoeffs);
512  m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
513 
514  // Gather information about trace interpolations
515  map<TraceInterpPoints, vector<pair<int, int> >, cmpop> TraceInterpMap;
516  map<TraceInterpPoints, vector<pair<int, int> >, cmpop>::iterator it;
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 (it = TraceInterpMap.begin(); it != TraceInterpMap.end(); ++it, ++cnt1)
638  {
639  LibUtilities::PointsKey fromPointsKey0 = it->first.get<0>();
640  LibUtilities::PointsKey fromPointsKey1 = it->first.get<1>();
641  LibUtilities::PointsKey toPointsKey0 = it->first.get<2>();
642  LibUtilities::PointsKey toPointsKey1 = it->first.get<3>();
643 
644  bool fwdSet = false;
645  bool bwdSet = false;
646 
647  for (int f = 0; f < it->second.size(); ++f, ++cnt2)
648  {
649  n = it->second[f].first;
650  e = it->second[f].second;
651 
652  StdRegions::StdExpansionSharedPtr face = elmtToTrace[n][e];
653 
654  exp3d = (*exp)[n]->as<LocalRegions::Expansion3D>();
655  phys_offset = locExp.GetPhys_Offset(n);
656 
657  // mapping of new face order to one that loops over elmts
658  // then faces set up mapping of faces in standard cartesian
659  // order
660  exp3d->GetFacePhysMap(e, faceids);
661  nfacepts = exp3d->GetFaceNumPoints(e);
662  nfacepts1 = face->GetTotPoints();
663 
664  StdRegions::Orientation orient = exp3d->GetForient(e);
665 
666  exp3d->ReOrientFacePhysMap(elmtToTrace[n][e]->GetNverts(),
667  orient,
668  toPointsKey0.GetNumPoints(),
669  toPointsKey1.GetNumPoints(),
670  locTraceToTraceMap);
671 
672  int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
673 
674  if (LeftAdjacents[TraceOrder[n][e]])
675  {
676  for (int i = 0; i < nfacepts; ++i)
677  {
678  m_fieldToLocTraceMap[cntFwd + i] = phys_offset + faceids[i];
679  }
680 
681  for (int i = 0; i < nfacepts1; ++i)
682  {
683  m_LocTraceToTraceMap[0][cntFwd1 + i] =
684  offset + locTraceToTraceMap[i];
685  }
686 
687  cntFwd += nfacepts;
688  cntFwd1 += nfacepts1;
689  set = 0;
690  }
691  else
692  {
693  for (int i = 0; i < nfacepts; ++i)
694  {
696  phys_offset + faceids[i];
697  }
698 
699  for (int i = 0; i < nfacepts1; ++i)
700  {
701  m_LocTraceToTraceMap[1][cntBwd1 + i] =
702  offset + locTraceToTraceMap[i];
703  }
704 
705  cntBwd += nfacepts;
706  cntBwd1 += nfacepts1;
707  set = 1;
708  }
709 
710  m_interpNfaces[set][cnt1] += 1;
711 
712  if (((fwdSet == false) && (set == 0)) ||
713  ((bwdSet == false) && (set == 1)))
714  {
715  m_interpPoints[set][cnt1] = it->first;
716 
717  if (fromPointsKey0 == toPointsKey0)
718  {
719  if (fromPointsKey1 == toPointsKey1)
720  {
721  m_interpTrace[set][cnt1] = eNoInterp;
722  }
723  else
724  {
725  m_interpTrace[set][cnt1] = eInterpDir1;
726  m_interpTraceI1[set][cnt1] =
727  LibUtilities::PointsManager()[fromPointsKey1]->GetI(
728  toPointsKey1);
729 
730  // Check to see if we can just
731  // interpolate endpoint
732  if ((fromPointsKey1.GetPointsType() ==
734  (toPointsKey1.GetPointsType() ==
736  {
737  if (fromPointsKey1.GetNumPoints() + 1 ==
738  toPointsKey1.GetNumPoints())
739  {
740  m_interpTrace[set][cnt1] = eInterpEndPtDir1;
741  int fnp1 = fromPointsKey1.GetNumPoints();
742  int tnp1 = toPointsKey1.GetNumPoints();
743  m_interpEndPtI1[set][cnt1] =
745  Vmath::Vcopy(
746  fnp1,
747  m_interpTraceI1[set][cnt1]->GetPtr().get() +
748  tnp1 - 1,
749  tnp1,
750  &m_interpEndPtI1[set][cnt1][0],
751  1);
752  }
753  }
754  }
755  }
756  else
757  {
758  if (fromPointsKey1 == toPointsKey1)
759  {
760  m_interpTrace[set][cnt1] = eInterpDir0;
761  m_interpTraceI0[set][cnt1] =
762  LibUtilities::PointsManager()[fromPointsKey0]->GetI(
763  toPointsKey0);
764 
765  // Check to see if we can just
766  // interpolate endpoint
767  if ((fromPointsKey0.GetPointsType() ==
769  (toPointsKey0.GetPointsType() ==
771  {
772  if (fromPointsKey0.GetNumPoints() + 1 ==
773  toPointsKey0.GetNumPoints())
774  {
775  m_interpTrace[set][cnt1] = eInterpEndPtDir0;
776  int fnp0 = fromPointsKey0.GetNumPoints();
777  int tnp0 = toPointsKey0.GetNumPoints();
778  m_interpEndPtI0[set][cnt1] =
780  Vmath::Vcopy(
781  fnp0,
782  m_interpTraceI0[set][cnt1]->GetPtr().get() +
783  tnp0 - 1,
784  tnp0,
785  &m_interpEndPtI0[set][cnt1][0],
786  1);
787  }
788  }
789  }
790  else
791  {
792  m_interpTrace[set][cnt1] = eInterpBothDirs;
793  m_interpTraceI0[set][cnt1] =
794  LibUtilities::PointsManager()[fromPointsKey0]->GetI(
795  toPointsKey0);
796  m_interpTraceI1[set][cnt1] =
797  LibUtilities::PointsManager()[fromPointsKey1]->GetI(
798  toPointsKey1);
799 
800  // check to see if we can just
801  // interpolate endpoint
802  if ((fromPointsKey0.GetPointsType() ==
804  (toPointsKey0.GetPointsType() ==
806  {
807  if (fromPointsKey0.GetNumPoints() + 1 ==
808  toPointsKey0.GetNumPoints())
809  {
810  m_interpTrace[set][cnt1] =
812  int fnp0 = fromPointsKey0.GetNumPoints();
813  int tnp0 = toPointsKey0.GetNumPoints();
814  m_interpEndPtI0[set][cnt1] =
816  Vmath::Vcopy(
817  fnp0,
818  m_interpTraceI0[set][cnt1]->GetPtr().get() +
819  tnp0 - 1,
820  tnp0,
821  &m_interpEndPtI0[set][cnt1][0],
822  1);
823  }
824  }
825  }
826  }
827 
828  if (set == 0)
829  {
830  fwdSet = true;
831  }
832  else
833  {
834  bwdSet = true;
835  }
836  }
837  }
838  }
839 }
840 
841 /**
842  * @brief Gather the local traces in physical space from field using
843  * #m_fieldToLocTraceMap.
844  *
845  * @param field Solution field in physical space
846  * @param faces Resulting local traces.
847  */
850 {
851  Vmath::Gathr(m_fieldToLocTraceMap.num_elements(),
852  field,
854  faces);
855 }
856 
857 /**
858  * @brief Gather the forwards-oriented local traces in physical space from field
859  * using #m_fieldToLocTraceMap.
860  *
861  * @param field Solution field in physical space
862  * @param faces Resulting local forwards-oriented traces.
863  */
866 {
868 }
869 
870 /**
871  * @brief Interpolate local trace edges to global trace edge point distributions
872  * where required.
873  *
874  * @param dir Selects forwards (0) or backwards (1) direction.
875  * @param locfaces Local trace edge storage.
876  * @param faces Global trace edge storage
877  */
879  const int dir,
880  const Array<OneD, const NekDouble> &locedges,
882 {
883  ASSERTL1(dir < 2,
884  "option dir out of range, "
885  " dir=0 is fwd, dir=1 is bwd");
886 
887  int cnt = 0;
888  int cnt1 = 0;
889 
890  // tmp space assuming forward map is of size of trace
892 
893  for (int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
894  {
895  // Check if there are edges to interpolate
896  if (m_interpNfaces[dir][i])
897  {
898  // Get to/from points
899  LibUtilities::PointsKey fromPointsKey0 =
900  m_interpPoints[dir][i].get<0>();
901  LibUtilities::PointsKey toPointsKey0 =
902  m_interpPoints[dir][i].get<2>();
903 
904  int fnp = fromPointsKey0.GetNumPoints();
905  int tnp = toPointsKey0.GetNumPoints();
906  int nedges = m_interpNfaces[dir][i];
907 
908  // Do interpolation here if required
909  switch (m_interpTrace[dir][i])
910  {
911  case eNoInterp: // Just copy
912  {
913  Vmath::Vcopy(nedges * fnp,
914  locedges.get() + cnt,
915  1,
916  tmp.get() + cnt1,
917  1);
918  }
919  break;
920  case eInterpDir0:
921  {
922  DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
923  Blas::Dgemm('N',
924  'N',
925  tnp,
926  nedges,
927  fnp,
928  1.0,
929  I0->GetPtr().get(),
930  tnp,
931  locedges.get() + cnt,
932  fnp,
933  0.0,
934  tmp.get() + cnt1,
935  tnp);
936  }
937  break;
938  case eInterpEndPtDir0:
939  {
941 
942  for (int k = 0; k < nedges; ++k)
943  {
944  Vmath::Vcopy(fnp,
945  &locedges[cnt + k * fnp],
946  1,
947  &tmp[cnt1 + k * tnp],
948  1);
949 
950  tmp[cnt1 + k * tnp + tnp - 1] = Blas::Ddot(
951  fnp, locedges.get() + cnt + k * fnp, 1, &I0[0], 1);
952  }
953  }
954  break;
955  default:
956  ASSERTL0(false,
957  "Invalid interpolation type for 2D elements");
958  break;
959  }
960 
961  cnt += nedges * fnp;
962  cnt1 += nedges * tnp;
963  }
964  }
965 
966  Vmath::Scatr(m_LocTraceToTraceMap[dir].num_elements(),
967  tmp.get(),
968  m_LocTraceToTraceMap[dir].get(),
969  edges.get());
970 }
971 
972 /**
973  * @brief Interpolate local faces to trace face point distributions where
974  * required.
975  *
976  * @param dir Selects forwards (0) or backwards (1) direction.
977  * @param locfaces Local trace face storage.
978  * @param faces Global trace face storage
979  */
981  const int dir,
982  const Array<OneD, const NekDouble> &locfaces,
984 {
985  ASSERTL1(dir < 2,
986  "option dir out of range, "
987  " dir=0 is fwd, dir=1 is bwd");
988 
989  int cnt1 = 0;
990  int cnt = 0;
991 
992  // tmp space assuming forward map is of size of trace
994 
995  for (int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
996  {
997  // Check if there are faces to interpolate
998  if (m_interpNfaces[dir][i])
999  {
1000  // Get to/from points
1001  LibUtilities::PointsKey fromPointsKey0 =
1002  m_interpPoints[dir][i].get<0>();
1003  LibUtilities::PointsKey fromPointsKey1 =
1004  m_interpPoints[dir][i].get<1>();
1005  LibUtilities::PointsKey toPointsKey0 =
1006  m_interpPoints[dir][i].get<2>();
1007  LibUtilities::PointsKey toPointsKey1 =
1008  m_interpPoints[dir][i].get<3>();
1009 
1010  int fnp0 = fromPointsKey0.GetNumPoints();
1011  int fnp1 = fromPointsKey1.GetNumPoints();
1012  int tnp0 = toPointsKey0.GetNumPoints();
1013  int tnp1 = toPointsKey1.GetNumPoints();
1014  int nfromfacepts = m_interpNfaces[dir][i] * fnp0 * fnp1;
1015  ;
1016 
1017  // Do interpolation here if required
1018  switch (m_interpTrace[dir][i])
1019  {
1020  case eNoInterp: // Just copy
1021  {
1022  Vmath::Vcopy(nfromfacepts,
1023  locfaces.get() + cnt,
1024  1,
1025  tmp.get() + cnt1,
1026  1);
1027  }
1028  break;
1029  case eInterpDir0:
1030  {
1031  DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1032  Blas::Dgemm('N',
1033  'N',
1034  tnp0,
1035  tnp1,
1036  fnp0,
1037  1.0,
1038  I0->GetPtr().get(),
1039  tnp0,
1040  locfaces.get() + cnt,
1041  fnp0,
1042  0.0,
1043  tmp.get() + cnt1,
1044  tnp0);
1045  }
1046  break;
1047  case eInterpEndPtDir0:
1048  {
1049  int nfaces = m_interpNfaces[dir][i];
1050  for (int k = 0; k < fnp0; ++k)
1051  {
1052  Vmath::Vcopy(nfaces * fnp1,
1053  locfaces.get() + cnt + k,
1054  fnp0,
1055  tmp.get() + cnt1 + k,
1056  tnp0);
1057  }
1059  Blas::Dgemv('T',
1060  fnp0,
1061  tnp1 * m_interpNfaces[dir][i],
1062  1.0,
1063  tmp.get() + cnt1,
1064  tnp0,
1065  I0.get(),
1066  1,
1067  0.0,
1068  tmp.get() + cnt1 + tnp0 - 1,
1069  tnp0);
1070  }
1071  break;
1072  case eInterpDir1:
1073  {
1074  DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1075  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1076  {
1077  Blas::Dgemm('N',
1078  'T',
1079  tnp0,
1080  tnp1,
1081  fnp1,
1082  1.0,
1083  locfaces.get() + cnt + j * fnp0 * fnp1,
1084  tnp0,
1085  I1->GetPtr().get(),
1086  tnp1,
1087  0.0,
1088  tmp.get() + cnt1 + j * tnp0 * tnp1,
1089  tnp0);
1090  }
1091  }
1092  break;
1093  case eInterpEndPtDir1:
1094  {
1096  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1097  {
1098  // copy all points
1099  Vmath::Vcopy(fnp0 * fnp1,
1100  locfaces.get() + cnt + j * fnp0 * fnp1,
1101  1,
1102  tmp.get() + cnt1 + j * tnp0 * tnp1,
1103  1);
1104 
1105  // interpolate end points
1106  for (int k = 0; k < tnp0; ++k)
1107  {
1108  tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1109  Blas::Ddot(fnp1,
1110  locfaces.get() + cnt +
1111  j * fnp0 * fnp1 + k,
1112  fnp0,
1113  &I1[0],
1114  1);
1115  }
1116  }
1117  }
1118  break;
1119  case eInterpBothDirs:
1120  {
1121  DNekMatSharedPtr I0 = m_interpTraceI0[dir][i];
1122  DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1123  Array<OneD, NekDouble> wsp(m_interpNfaces[dir][i] * fnp0 *
1124  tnp1 * fnp0);
1125 
1126  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1127  {
1128  Blas::Dgemm('N',
1129  'T',
1130  fnp0,
1131  tnp1,
1132  fnp1,
1133  1.0,
1134  locfaces.get() + cnt + j * fnp0 * fnp1,
1135  fnp0,
1136  I1->GetPtr().get(),
1137  tnp1,
1138  0.0,
1139  wsp.get() + j * fnp0 * tnp1,
1140  fnp0);
1141  }
1142  Blas::Dgemm('N',
1143  'N',
1144  tnp0,
1145  tnp1 * m_interpNfaces[dir][i],
1146  fnp0,
1147  1.0,
1148  I0->GetPtr().get(),
1149  tnp0,
1150  wsp.get(),
1151  fnp0,
1152  0.0,
1153  tmp.get() + cnt1,
1154  tnp0);
1155  }
1156  break;
1158  {
1159  DNekMatSharedPtr I1 = m_interpTraceI1[dir][i];
1160 
1161  for (int j = 0; j < m_interpNfaces[dir][i]; ++j)
1162  {
1163  Blas::Dgemm('N',
1164  'T',
1165  fnp0,
1166  tnp1,
1167  fnp1,
1168  1.0,
1169  locfaces.get() + cnt + j * fnp0 * fnp1,
1170  fnp0,
1171  I1->GetPtr().get(),
1172  tnp1,
1173  0.0,
1174  tmp.get() + cnt1 + j * tnp0 * tnp1,
1175  tnp0);
1176  }
1177 
1179  Blas::Dgemv('T',
1180  fnp0,
1181  tnp1 * m_interpNfaces[dir][i],
1182  1.0,
1183  tmp.get() + cnt1,
1184  tnp0,
1185  I0.get(),
1186  1,
1187  0.0,
1188  tmp.get() + cnt1 + tnp0 - 1,
1189  tnp0);
1190  }
1191  break;
1192  }
1193  cnt += nfromfacepts;
1194  cnt1 += m_interpNfaces[dir][i] * tnp0 * tnp1;
1195  }
1196  }
1197 
1198  Vmath::Scatr(m_LocTraceToTraceMap[dir].num_elements(),
1199  tmp.get(),
1200  m_LocTraceToTraceMap[dir].get(),
1201  faces.get());
1202 }
1203 
1204 /**
1205  * @brief Add contributions from trace coefficients to the elemental field
1206  * storage.
1207  *
1208  * @param trace Array of global trace coefficients.
1209  * @param field Array containing field coefficients storage.
1210  */
1213 {
1214  int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
1215  for (int i = 0; i < nvals; ++i)
1216  {
1217  field[m_traceCoeffsToElmtMap[0][i]] +=
1218  m_traceCoeffsToElmtSign[0][i] *
1219  trace[m_traceCoeffsToElmtTrace[0][i]];
1220  }
1221 }
1222 
1223 /**
1224  * @brief Add contributions from backwards or forwards oriented trace
1225  * coefficients to the elemental field storage.
1226  *
1227  * @param dir Selects forwards (0) or backwards (1) direction
1228  * @param trace Array of global trace coefficients.
1229  * @param field Array containing field coefficients storage.
1230  */
1232  const int dir,
1233  const Array<OneD, const NekDouble> &trace,
1234  Array<OneD, NekDouble> &field)
1235 {
1236  int nvals = m_nTraceCoeffs[dir];
1237  for (int i = 0; i < nvals; ++i)
1238  {
1239  field[m_traceCoeffsToElmtMap[dir][i]] +=
1240  m_traceCoeffsToElmtSign[dir][i] *
1241  trace[m_traceCoeffsToElmtTrace[dir][i]];
1242  }
1243 }
1244 
1245 }
1246 }
int m_nTracePts
The number of global trace points.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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:1926
boost::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 ...
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:630
Array< OneD, Array< OneD, DNekMatSharedPtr > > m_interpTraceI0
Interpolation matrices for either 2D edges or first coordinate of 3D face.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
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:1934
void FwdLocTracesFromField(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > faces)
Gather the forwards-oriented local traces in physical space from field using m_fieldToLocTraceMap.
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:1917
Array< OneD, Array< OneD, int > > m_traceCoeffsToElmtTrace
Mapping from forwards/backwards trace coefficients to the position of the trace element in global sto...
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
void InterpLocEdgesToTrace(const int dir, const Array< OneD, const NekDouble > &locfaces, Array< OneD, NekDouble > edges)
Interpolate local trace edges to global trace edge point distributions where required.
LocTraceToTraceMap(const ExpList &locExp, const ExpListSharedPtr &trace, const Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > &elmtToTrace, const vector< bool > &LeftAdjacents)
Set up trace to trace mapping components.
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_interpEndPtI0
Mapping to hold first coordinate direction endpoint interpolation, which can be more optimal if using...
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
Array< OneD, Array< OneD, int > > m_traceCoeffsToElmtSign
Sign array for mapping from forwards/backwards trace coefficients to local trace storage.
int m_nLocTracePts
The number of local trace points.
int m_nTraceCoeffs[2]
Number of forwards/backwards trace coefficients.
boost::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:48
int m_nFwdLocTracePts
The number of forward trace points. A local trace element is `forward' if it is the side selected for...
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
void Setup2D(const ExpList &locExp, const ExpListSharedPtr &trace, const Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > &elmtToTrace, const vector< bool > &LeftAdjacents)
Set up member variables for a two-dimensional problem.
Array< OneD, Array< OneD, int > > m_traceCoeffsToElmtMap
Mapping from forwards/backwards trace coefficients to elemental coefficient storage.
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:659
Defines a specification for a set of points.
Definition: Points.h:58
void LocTracesFromField(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > faces)
Gather the local traces in physical space from field using m_fieldToLocTraceMap.
Array< OneD, Array< OneD, DNekMatSharedPtr > > m_interpTraceI1
Interpolation matrices for the second coordinate of 3D face, not used in 2D.
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:434
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Array< OneD, Array< OneD, int > > m_interpNfaces
Number of edges/faces on a 2D/3D element that require interpolation.
unsigned int GetNumPoints() const
Definition: Points.h:106
Array< OneD, Array< OneD, int > > m_LocTraceToTraceMap
A mapping from local trace points to the global trace. Dimension 0 holds forward traces, dimension 1 backward.
void Setup3D(const ExpList &locExp, const ExpListSharedPtr &trace, const Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > &elmtToTrace, const vector< bool > &LeftAdjacents)
Set up member variables for a three-dimensional problem.
Array< OneD, Array< OneD, InterpLocTraceToTrace > > m_interpTrace
A mapping holding the type of interpolation needed for each local trace. Dimension 0 holds forward tr...
void AddTraceCoeffsToFieldCoeffs(const Array< OneD, const NekDouble > &trace, Array< OneD, NekDouble > &field)
Add contributions from trace coefficients to the elemental field storage.
Array< OneD, Array< OneD, TraceInterpPoints > > m_interpPoints
Interpolation points key distributions for each of the local to global mappings.
void InterpLocFacesToTrace(const int dir, const Array< OneD, const NekDouble > &locfaces, Array< OneD, NekDouble > faces)
Interpolate local faces to trace face point distributions where required.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Array< OneD, int > m_fieldToLocTraceMap
A mapping from the local trace points, arranged as all forwards traces followed by backwards traces...
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
PointsType GetPointsType() const
Definition: Points.h:111
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_interpEndPtI1
Mapping to hold second coordinate direction endpoint interpolation, which can be more optimal if usin...
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49