51 namespace MultiRegions
64 LocTraceToTraceMap::LocTraceToTraceMap(
69 const vector<bool> &LeftAdjacents)
74 int m_expdim = locExpVector[0]->GetShapeDimension();
82 Setup2D(locExp, trace, elmtToTrace, LeftAdjacents);
85 Setup3D(locExp, trace, elmtToTrace, LeftAdjacents);
88 ASSERTL0(
false,
"Number of dimensions greater than 3")
93 LocTraceToTraceMap::~LocTraceToTraceMap()
105 void LocTraceToTraceMap::Setup2D(
110 const vector<bool> &LeftAdjacents)
124 const boost::shared_ptr<LocalRegions::ExpansionVector> exp =
127 int cnt, n, e, phys_offset;
129 int nexp = exp->size();
130 m_nTracePts = trace->GetTotPoints();
137 m_nFwdLocTracePts = 0;
140 for (cnt = n = 0; n < nexp; ++n)
144 for (
int i = 0; i < exp2d->GetNedges(); ++i, ++cnt)
146 int nLocPts = exp2d->GetEdgeNumPoints(i);
147 m_nLocTracePts += nLocPts;
149 if (LeftAdjacents[cnt])
151 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
152 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
153 m_nFwdLocTracePts += nLocPts;
157 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
158 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
168 m_nTraceCoeffs[0] = nFwdCoeffs;
169 m_nTraceCoeffs[1] = nBwdCoeffs;
172 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
174 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
176 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
179 map<TraceInterpPoints, vector<pair<int, int> >,
cmpop> TraceInterpMap;
182 vector<vector<int> > TraceOrder;
183 TraceOrder.resize(nexp);
190 for (cnt = n = 0; n < nexp; ++n)
193 nedge = exp2d->GetNedges();
194 TraceOrder[n].resize(nedge);
197 for (e = 0; e < nedge; ++e, ++cnt)
209 basisDir = e == 0 ? 0 : 1;
216 fromPointsKey0 = exp2d->GetBasis(basisDir)->GetPointsKey();
219 toPointsKey0 = edge->GetBasis(0)->GetPointsKey();
225 fromPointsKey0, fromPointsKey1, toPointsKey0, toPointsKey1);
227 pair<int, int> epf(n, e);
228 TraceInterpMap[fpoint].push_back(epf);
229 TraceOrder[n][e] = cnt;
236 exp2d->GetEdgeToElementMap(
237 e, orient, map, sign, edge->GetBasisNumModes(0));
239 int order_f = edge->GetNcoeffs();
240 int foffset = trace->GetCoeff_Offset(edge->GetElmtId());
242 double fac = (*exp)[n]->EdgeNormalNegated(e) ? -1.0 : 1.0;
247 if (exp2d->GetEdgeExp(e)->GetRightAdjacentElementExp())
249 if (locExp1d->GetRightAdjacentElementExp()
251 ->GetGlobalID() == exp2d->GetGeom()->GetGlobalID())
257 if (LeftAdjacents[cnt])
259 for (
int i = 0; i < order_f; ++i)
261 m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
262 m_traceCoeffsToElmtTrace[0][fwdcnt] = foffset + i;
263 m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
268 for (
int i = 0; i < order_f; ++i)
270 m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
271 m_traceCoeffsToElmtTrace[1][bwdcnt] = foffset + i;
272 m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
278 int nInterpType = TraceInterpMap.size();
280 for (
int i = 0; i < 2; ++i)
289 int nedgepts, nedgepts1;
301 for (it = TraceInterpMap.begin(); it != TraceInterpMap.end(); ++it, ++cnt1)
309 for (
int f = 0; f < it->second.size(); ++f, ++cnt2)
311 n = it->second[f].first;
312 e = it->second[f].second;
321 exp2d->GetEdgePhysMap(e, edgeids);
323 nedgepts = exp2d->GetEdgeNumPoints(e);
324 nedgepts1 = edge->GetTotPoints();
329 exp2d->ReOrientEdgePhysMap(elmtToTrace[n][e]->GetNverts(),
334 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
336 if (LeftAdjacents[TraceOrder[n][e]])
338 for (
int i = 0; i < nedgepts; ++i)
340 m_fieldToLocTraceMap[cntFwd + i] = phys_offset + edgeids[i];
343 for (
int i = 0; i < nedgepts1; ++i)
345 m_LocTraceToTraceMap[0][cntFwd1 + i] =
346 offset + locTraceToTraceMap[i];
350 cntFwd1 += nedgepts1;
355 for (
int i = 0; i < nedgepts; ++i)
357 m_fieldToLocTraceMap[m_nFwdLocTracePts + cntBwd + i] =
358 phys_offset + edgeids[i];
361 for (
int i = 0; i < nedgepts1; ++i)
363 m_LocTraceToTraceMap[1][cntBwd1 + i] =
364 offset + locTraceToTraceMap[i];
368 cntBwd1 += nedgepts1;
372 m_interpNfaces[set][cnt1] += 1;
374 if ((fwdSet ==
false && set == 0) || (bwdSet ==
false && set == 1))
376 m_interpPoints[set][cnt1] = it->first;
378 if (fromPointsKey0 == toPointsKey0)
385 m_interpTraceI0[set][cnt1] =
405 m_interpEndPtI0[set][cnt1] =
410 m_interpTraceI0[set][cnt1]->GetPtr().
get() +
413 &m_interpEndPtI0[set][cnt1][0],
440 void LocTraceToTraceMap::Setup3D(
445 const vector<bool> &LeftAdjacents)
462 const boost::shared_ptr<LocalRegions::ExpansionVector> exp =
465 int cnt, n, e, phys_offset;
467 int nexp = exp->size();
468 m_nTracePts = trace->GetTotPoints();
475 m_nFwdLocTracePts = 0;
478 for (cnt = n = 0; n < nexp; ++n)
482 for (
int i = 0; i < exp3d->GetNfaces(); ++i, ++cnt)
484 int nLocPts = exp3d->GetFaceNumPoints(i);
485 m_nLocTracePts += nLocPts;
487 if (LeftAdjacents[cnt])
489 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
490 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
491 m_nFwdLocTracePts += nLocPts;
495 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
496 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
506 m_nTraceCoeffs[0] = nFwdCoeffs;
507 m_nTraceCoeffs[1] = nBwdCoeffs;
510 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
512 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
514 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
517 map<TraceInterpPoints, vector<pair<int, int> >,
cmpop> TraceInterpMap;
520 vector<vector<int> > TraceOrder;
521 TraceOrder.resize(nexp);
528 for (cnt = n = 0; n < nexp; ++n)
531 nface = exp3d->GetNfaces();
532 TraceOrder[n].resize(nface);
535 for (e = 0; e < nface; ++e, ++cnt)
544 int dir0 = exp3d->GetGeom3D()->GetDir(e, 0);
545 int dir1 = exp3d->GetGeom3D()->GetDir(e, 1);
547 fromPointsKey0 = exp3d->GetBasis(dir0)->GetPointsKey();
548 fromPointsKey1 = exp3d->GetBasis(dir1)->GetPointsKey();
552 toPointsKey0 = face->GetBasis(0)->GetPointsKey();
553 toPointsKey1 = face->GetBasis(1)->GetPointsKey();
557 toPointsKey0 = face->GetBasis(1)->GetPointsKey();
558 toPointsKey1 = face->GetBasis(0)->GetPointsKey();
562 fromPointsKey0, fromPointsKey1, toPointsKey0, toPointsKey1);
564 pair<int, int> epf(n, e);
565 TraceInterpMap[fpoint].push_back(epf);
566 TraceOrder[n][e] = cnt;
572 exp3d->GetFaceToElementMap(e,
576 face->GetBasisNumModes(0),
577 face->GetBasisNumModes(1));
579 int order_f = face->GetNcoeffs();
580 int foffset = trace->GetCoeff_Offset(face->GetElmtId());
582 int fac = (*exp)[n]->FaceNormalNegated(e) ? -1.0 : 1.0;
584 if (exp3d->GetFaceExp(e)->GetRightAdjacentElementExp())
586 if (exp3d->GetFaceExp(e)
587 ->GetRightAdjacentElementExp()
589 ->GetGlobalID() == exp3d->GetGeom3D()->GetGlobalID())
595 if (LeftAdjacents[cnt])
597 for (
int i = 0; i < order_f; ++i)
599 m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
600 m_traceCoeffsToElmtTrace[0][fwdcnt] = foffset + i;
601 m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
606 for (
int i = 0; i < order_f; ++i)
608 m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
609 m_traceCoeffsToElmtTrace[1][bwdcnt] = foffset + i;
610 m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
616 int nInterpType = TraceInterpMap.size();
617 for (
int i = 0; i < 2; ++i)
628 int nfacepts, nfacepts1;
639 for (it = TraceInterpMap.begin(); it != TraceInterpMap.end(); ++it, ++cnt1)
649 for (
int f = 0; f < it->second.size(); ++f, ++cnt2)
651 n = it->second[f].first;
652 e = it->second[f].second;
662 exp3d->GetFacePhysMap(e, faceids);
663 nfacepts = exp3d->GetFaceNumPoints(e);
664 nfacepts1 = face->GetTotPoints();
668 exp3d->ReOrientFacePhysMap(elmtToTrace[n][e]->GetNverts(),
674 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
676 if (LeftAdjacents[TraceOrder[n][e]])
678 for (
int i = 0; i < nfacepts; ++i)
680 m_fieldToLocTraceMap[cntFwd + i] = phys_offset + faceids[i];
683 for (
int i = 0; i < nfacepts1; ++i)
685 m_LocTraceToTraceMap[0][cntFwd1 + i] =
686 offset + locTraceToTraceMap[i];
690 cntFwd1 += nfacepts1;
695 for (
int i = 0; i < nfacepts; ++i)
697 m_fieldToLocTraceMap[m_nFwdLocTracePts + cntBwd + i] =
698 phys_offset + faceids[i];
701 for (
int i = 0; i < nfacepts1; ++i)
703 m_LocTraceToTraceMap[1][cntBwd1 + i] =
704 offset + locTraceToTraceMap[i];
708 cntBwd1 += nfacepts1;
712 m_interpNfaces[set][cnt1] += 1;
714 if (((fwdSet ==
false) && (set == 0)) ||
715 ((bwdSet ==
false) && (set == 1)))
717 m_interpPoints[set][cnt1] = it->first;
719 if (fromPointsKey0 == toPointsKey0)
721 if (fromPointsKey1 == toPointsKey1)
728 m_interpTraceI1[set][cnt1] =
745 m_interpEndPtI1[set][cnt1] =
749 m_interpTraceI1[set][cnt1]->GetPtr().
get() +
752 &m_interpEndPtI1[set][cnt1][0],
760 if (fromPointsKey1 == toPointsKey1)
763 m_interpTraceI0[set][cnt1] =
780 m_interpEndPtI0[set][cnt1] =
784 m_interpTraceI0[set][cnt1]->GetPtr().
get() +
787 &m_interpEndPtI0[set][cnt1][0],
795 m_interpTraceI0[set][cnt1] =
798 m_interpTraceI1[set][cnt1] =
812 m_interpTrace[set][cnt1] =
816 m_interpEndPtI0[set][cnt1] =
820 m_interpTraceI0[set][cnt1]->GetPtr().
get() +
823 &m_interpEndPtI0[set][cnt1][0],
850 void LocTraceToTraceMap::LocTracesFromField(
855 m_fieldToLocTraceMap,
866 void LocTraceToTraceMap::FwdLocTracesFromField(
869 Vmath::Gathr(m_nFwdLocTracePts, field, m_fieldToLocTraceMap, faces);
880 void LocTraceToTraceMap::InterpLocEdgesToTrace(
886 "option dir out of range, "
887 " dir=0 is fwd, dir=1 is bwd");
895 for (
int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
898 if (m_interpNfaces[dir][i])
902 m_interpPoints[dir][i].get<0>();
904 m_interpPoints[dir][i].get<2>();
908 int nedges = m_interpNfaces[dir][i];
911 switch (m_interpTrace[dir][i])
916 locedges.get() + cnt,
933 locedges.get() + cnt,
944 for (
int k = 0; k < nedges; ++k)
947 &locedges[cnt + k * fnp],
949 &tmp[cnt1 + k * tnp],
953 fnp, locedges.get() + cnt + k * fnp, 1, &I0[0], 1);
959 "Invalid interpolation type for 2D elements");
964 cnt1 += nedges * tnp;
970 m_LocTraceToTraceMap[dir].get(),
982 void LocTraceToTraceMap::InterpLocFacesToTrace(
988 "option dir out of range, "
989 " dir=0 is fwd, dir=1 is bwd");
997 for (
int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
1000 if (m_interpNfaces[dir][i])
1004 m_interpPoints[dir][i].get<0>();
1006 m_interpPoints[dir][i].get<1>();
1008 m_interpPoints[dir][i].get<2>();
1010 m_interpPoints[dir][i].get<3>();
1016 int nfromfacepts = m_interpNfaces[dir][i] * fnp0 * fnp1;
1020 switch (m_interpTrace[dir][i])
1025 locfaces.get() + cnt,
1042 locfaces.get() + cnt,
1051 int nfaces = m_interpNfaces[dir][i];
1052 for (
int k = 0; k < fnp0; ++k)
1055 locfaces.get() + cnt + k,
1057 tmp.get() + cnt1 + k,
1063 tnp1 * m_interpNfaces[dir][i],
1070 tmp.get() + cnt1 + tnp0 - 1,
1077 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1085 locfaces.get() + cnt + j * fnp0 * fnp1,
1090 tmp.get() + cnt1 + j * tnp0 * tnp1,
1098 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1102 locfaces.get() + cnt + j * fnp0 * fnp1,
1104 tmp.get() + cnt1 + j * tnp0 * tnp1,
1108 for (
int k = 0; k < tnp0; ++k)
1110 tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1112 locfaces.get() + cnt +
1113 j * fnp0 * fnp1 + k,
1128 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1136 locfaces.get() + cnt + j * fnp0 * fnp1,
1141 wsp.get() + j * fnp0 * tnp1,
1147 tnp1 * m_interpNfaces[dir][i],
1163 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1171 locfaces.get() + cnt + j * fnp0 * fnp1,
1176 tmp.get() + cnt1 + j * tnp0 * tnp1,
1183 tnp1 * m_interpNfaces[dir][i],
1190 tmp.get() + cnt1 + tnp0 - 1,
1195 cnt += nfromfacepts;
1196 cnt1 += m_interpNfaces[dir][i] * tnp0 * tnp1;
1202 m_LocTraceToTraceMap[dir].get(),
1213 void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1216 int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
1217 for (
int i = 0; i < nvals; ++i)
1219 field[m_traceCoeffsToElmtMap[0][i]] +=
1220 m_traceCoeffsToElmtSign[0][i] *
1221 trace[m_traceCoeffsToElmtTrace[0][i]];
1233 void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1238 int nvals = m_nTraceCoeffs[dir];
1239 for (
int i = 0; i < nvals; ++i)
1241 field[m_traceCoeffsToElmtMap[dir][i]] +=
1242 m_traceCoeffsToElmtSign[dir][i] *
1243 trace[m_traceCoeffsToElmtTrace[dir][i]];
#define ASSERTL0(condition, msg)
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
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]].
#define sign(a, b)
return the sign(b)*a
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
std::vector< ExpansionSharedPtr > ExpansionVector
Gauss Radau pinned at x=-1, .
Base class for all multi-elemental spectral/hp expansions.
boost::shared_ptr< Expansion3D > Expansion3DSharedPtr
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].
Defines a specification for a set of points.
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)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
unsigned int GetNumPoints() const
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
1D Gauss-Lobatto-Legendre quadrature points
PointsType GetPointsType() const
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr