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