50 namespace MultiRegions
63 LocTraceToTraceMap::LocTraceToTraceMap(
68 const vector<bool> &LeftAdjacents)
73 int m_expdim = locExpVector[0]->GetShapeDimension();
81 Setup2D(locExp, trace, elmtToTrace, LeftAdjacents);
84 Setup3D(locExp, trace, elmtToTrace, LeftAdjacents);
87 ASSERTL0(
false,
"Number of dimensions greater than 3")
92 LocTraceToTraceMap::~LocTraceToTraceMap()
104 void LocTraceToTraceMap::Setup2D(
109 const vector<bool> &LeftAdjacents)
123 const std::shared_ptr<LocalRegions::ExpansionVector> exp =
126 int cnt, n, e, phys_offset;
128 int nexp = exp->size();
129 m_nTracePts = trace->GetTotPoints();
136 m_nFwdLocTracePts = 0;
139 for (cnt = n = 0; n < nexp; ++n)
143 for (
int i = 0; i < exp2d->GetNedges(); ++i, ++cnt)
145 int nLocPts = exp2d->GetEdgeNumPoints(i);
146 m_nLocTracePts += nLocPts;
148 if (LeftAdjacents[cnt])
150 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
151 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
152 m_nFwdLocTracePts += nLocPts;
156 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
157 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
167 m_nTraceCoeffs[0] = nFwdCoeffs;
168 m_nTraceCoeffs[1] = nBwdCoeffs;
171 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
173 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
175 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
178 map<TraceInterpPoints, vector<pair<int, int> >,
cmpop> TraceInterpMap;
180 vector<vector<int> > TraceOrder;
181 TraceOrder.resize(nexp);
188 for (cnt = n = 0; n < nexp; ++n)
191 nedge = exp2d->GetNedges();
192 TraceOrder[n].resize(nedge);
195 for (e = 0; e < nedge; ++e, ++cnt)
207 basisDir = e == 0 ? 0 : 1;
214 fromPointsKey0 = exp2d->GetBasis(basisDir)->GetPointsKey();
217 toPointsKey0 = edge->GetBasis(0)->GetPointsKey();
223 fromPointsKey0, fromPointsKey1, toPointsKey0, toPointsKey1);
225 pair<int, int> epf(n, e);
226 TraceInterpMap[fpoint].push_back(epf);
227 TraceOrder[n][e] = cnt;
234 exp2d->GetEdgeToElementMap(
235 e, orient, map, sign, edge->GetBasisNumModes(0));
237 int order_f = edge->GetNcoeffs();
238 int foffset = trace->GetCoeff_Offset(edge->GetElmtId());
240 double fac = (*exp)[n]->EdgeNormalNegated(e) ? -1.0 : 1.0;
245 if (exp2d->GetEdgeExp(e)->GetRightAdjacentElementExp())
247 if (locExp1d->GetRightAdjacentElementExp()
249 ->GetGlobalID() == exp2d->GetGeom()->GetGlobalID())
255 if (LeftAdjacents[cnt])
257 for (
int i = 0; i < order_f; ++i)
259 m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
260 m_traceCoeffsToElmtTrace[0][fwdcnt] = foffset + i;
261 m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
266 for (
int i = 0; i < order_f; ++i)
268 m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
269 m_traceCoeffsToElmtTrace[1][bwdcnt] = foffset + i;
270 m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
276 int nInterpType = TraceInterpMap.size();
278 for (
int i = 0; i < 2; ++i)
287 int nedgepts, nedgepts1;
299 for (
auto it = TraceInterpMap.begin(); it != TraceInterpMap.end();
308 for (
int f = 0; f < it->second.size(); ++f, ++cnt2)
310 n = it->second[f].first;
311 e = it->second[f].second;
320 exp2d->GetEdgePhysMap(e, edgeids);
322 nedgepts = exp2d->GetEdgeNumPoints(e);
323 nedgepts1 = edge->GetTotPoints();
328 exp2d->ReOrientEdgePhysMap(elmtToTrace[n][e]->GetNverts(),
330 toPointsKey0.GetNumPoints(),
333 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
335 if (LeftAdjacents[TraceOrder[n][e]])
337 for (
int i = 0; i < nedgepts; ++i)
339 m_fieldToLocTraceMap[cntFwd + i] = phys_offset + edgeids[i];
342 for (
int i = 0; i < nedgepts1; ++i)
344 m_LocTraceToTraceMap[0][cntFwd1 + i] =
345 offset + locTraceToTraceMap[i];
349 cntFwd1 += nedgepts1;
354 for (
int i = 0; i < nedgepts; ++i)
356 m_fieldToLocTraceMap[m_nFwdLocTracePts + cntBwd + i] =
357 phys_offset + edgeids[i];
360 for (
int i = 0; i < nedgepts1; ++i)
362 m_LocTraceToTraceMap[1][cntBwd1 + i] =
363 offset + locTraceToTraceMap[i];
367 cntBwd1 += nedgepts1;
371 m_interpNfaces[
set][cnt1] += 1;
373 if ((fwdSet ==
false &&
set == 0) || (bwdSet ==
false &&
set == 1))
375 m_interpPoints[
set][cnt1] = it->first;
377 if (fromPointsKey0 == toPointsKey0)
384 m_interpTraceI0[
set][cnt1] =
392 (toPointsKey0.GetPointsType() ==
396 toPointsKey0.GetNumPoints())
402 int tnp0 = toPointsKey0.GetNumPoints();
404 m_interpEndPtI0[
set][cnt1] =
409 m_interpTraceI0[
set][cnt1]->GetPtr().
get() +
412 &m_interpEndPtI0[
set][cnt1][0],
439 void LocTraceToTraceMap::Setup3D(
444 const vector<bool> &LeftAdjacents)
461 const std::shared_ptr<LocalRegions::ExpansionVector> exp =
464 int cnt, n, e, phys_offset;
466 int nexp = exp->size();
467 m_nTracePts = trace->GetTotPoints();
474 m_nFwdLocTracePts = 0;
477 for (cnt = n = 0; n < nexp; ++n)
481 for (
int i = 0; i < exp3d->GetNfaces(); ++i, ++cnt)
483 int nLocPts = exp3d->GetFaceNumPoints(i);
484 m_nLocTracePts += nLocPts;
486 if (LeftAdjacents[cnt])
488 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
489 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
490 m_nFwdLocTracePts += nLocPts;
494 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
495 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
505 m_nTraceCoeffs[0] = nFwdCoeffs;
506 m_nTraceCoeffs[1] = nBwdCoeffs;
509 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
511 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
513 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
516 map<TraceInterpPoints, vector<pair<int, int> >,
cmpop> TraceInterpMap;
518 vector<vector<int> > TraceOrder;
519 TraceOrder.resize(nexp);
526 for (cnt = n = 0; n < nexp; ++n)
529 nface = exp3d->GetNfaces();
530 TraceOrder[n].resize(nface);
533 for (e = 0; e < nface; ++e, ++cnt)
542 int dir0 = exp3d->GetGeom3D()->GetDir(e, 0);
543 int dir1 = exp3d->GetGeom3D()->GetDir(e, 1);
545 fromPointsKey0 = exp3d->GetBasis(dir0)->GetPointsKey();
546 fromPointsKey1 = exp3d->GetBasis(dir1)->GetPointsKey();
550 toPointsKey0 = face->GetBasis(0)->GetPointsKey();
551 toPointsKey1 = face->GetBasis(1)->GetPointsKey();
555 toPointsKey0 = face->GetBasis(1)->GetPointsKey();
556 toPointsKey1 = face->GetBasis(0)->GetPointsKey();
560 fromPointsKey0, fromPointsKey1, toPointsKey0, toPointsKey1);
562 pair<int, int> epf(n, e);
563 TraceInterpMap[fpoint].push_back(epf);
564 TraceOrder[n][e] = cnt;
570 exp3d->GetFaceToElementMap(e,
574 face->GetBasisNumModes(0),
575 face->GetBasisNumModes(1));
577 int order_f = face->GetNcoeffs();
578 int foffset = trace->GetCoeff_Offset(face->GetElmtId());
580 int fac = (*exp)[n]->FaceNormalNegated(e) ? -1.0 : 1.0;
582 if (exp3d->GetFaceExp(e)->GetRightAdjacentElementExp())
584 if (exp3d->GetFaceExp(e)
585 ->GetRightAdjacentElementExp()
587 ->GetGlobalID() == exp3d->GetGeom3D()->GetGlobalID())
593 if (LeftAdjacents[cnt])
595 for (
int i = 0; i < order_f; ++i)
597 m_traceCoeffsToElmtMap[0][fwdcnt] = coeffoffset + map[i];
598 m_traceCoeffsToElmtTrace[0][fwdcnt] = foffset + i;
599 m_traceCoeffsToElmtSign[0][fwdcnt++] = fac * sign[i];
604 for (
int i = 0; i < order_f; ++i)
606 m_traceCoeffsToElmtMap[1][bwdcnt] = coeffoffset + map[i];
607 m_traceCoeffsToElmtTrace[1][bwdcnt] = foffset + i;
608 m_traceCoeffsToElmtSign[1][bwdcnt++] = fac * sign[i];
614 int nInterpType = TraceInterpMap.size();
615 for (
int i = 0; i < 2; ++i)
626 int nfacepts, nfacepts1;
637 for (
auto it = TraceInterpMap.begin(); it != TraceInterpMap.end();
648 for (
int f = 0; f < it->second.size(); ++f, ++cnt2)
650 n = it->second[f].first;
651 e = it->second[f].second;
661 exp3d->GetFacePhysMap(e, faceids);
662 nfacepts = exp3d->GetFaceNumPoints(e);
663 nfacepts1 = face->GetTotPoints();
667 exp3d->ReOrientFacePhysMap(elmtToTrace[n][e]->GetNverts(),
669 toPointsKey0.GetNumPoints(),
670 toPointsKey1.GetNumPoints(),
673 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
675 if (LeftAdjacents[TraceOrder[n][e]])
677 for (
int i = 0; i < nfacepts; ++i)
679 m_fieldToLocTraceMap[cntFwd + i] = phys_offset + faceids[i];
682 for (
int i = 0; i < nfacepts1; ++i)
684 m_LocTraceToTraceMap[0][cntFwd1 + i] =
685 offset + locTraceToTraceMap[i];
689 cntFwd1 += nfacepts1;
694 for (
int i = 0; i < nfacepts; ++i)
696 m_fieldToLocTraceMap[m_nFwdLocTracePts + cntBwd + i] =
697 phys_offset + faceids[i];
700 for (
int i = 0; i < nfacepts1; ++i)
702 m_LocTraceToTraceMap[1][cntBwd1 + i] =
703 offset + locTraceToTraceMap[i];
707 cntBwd1 += nfacepts1;
711 m_interpNfaces[
set][cnt1] += 1;
713 if (((fwdSet ==
false) && (
set == 0)) ||
714 ((bwdSet ==
false) && (
set == 1)))
716 m_interpPoints[
set][cnt1] = it->first;
718 if (fromPointsKey0 == toPointsKey0)
720 if (fromPointsKey1 == toPointsKey1)
727 m_interpTraceI1[
set][cnt1] =
733 if ((fromPointsKey1.GetPointsType() ==
735 (toPointsKey1.GetPointsType() ==
738 if (fromPointsKey1.GetNumPoints() + 1 ==
739 toPointsKey1.GetNumPoints())
742 int fnp1 = fromPointsKey1.GetNumPoints();
743 int tnp1 = toPointsKey1.GetNumPoints();
744 m_interpEndPtI1[
set][cnt1] =
748 m_interpTraceI1[
set][cnt1]->GetPtr().
get() +
751 &m_interpEndPtI1[
set][cnt1][0],
759 if (fromPointsKey1 == toPointsKey1)
762 m_interpTraceI0[
set][cnt1] =
770 (toPointsKey0.GetPointsType() ==
774 toPointsKey0.GetNumPoints())
778 int tnp0 = toPointsKey0.GetNumPoints();
779 m_interpEndPtI0[
set][cnt1] =
783 m_interpTraceI0[
set][cnt1]->GetPtr().
get() +
786 &m_interpEndPtI0[
set][cnt1][0],
794 m_interpTraceI0[
set][cnt1] =
797 m_interpTraceI1[
set][cnt1] =
805 (toPointsKey0.GetPointsType() ==
809 toPointsKey0.GetNumPoints())
811 m_interpTrace[
set][cnt1] =
814 int tnp0 = toPointsKey0.GetNumPoints();
815 m_interpEndPtI0[
set][cnt1] =
819 m_interpTraceI0[
set][cnt1]->GetPtr().
get() +
822 &m_interpEndPtI0[
set][cnt1][0],
849 void LocTraceToTraceMap::LocTracesFromField(
854 m_fieldToLocTraceMap,
865 void LocTraceToTraceMap::FwdLocTracesFromField(
868 Vmath::Gathr(m_nFwdLocTracePts, field, m_fieldToLocTraceMap, faces);
879 void LocTraceToTraceMap::InterpLocEdgesToTrace(
885 "option dir out of range, " 886 " dir=0 is fwd, dir=1 is bwd");
894 for (
int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
897 if (m_interpNfaces[dir][i])
901 std::get<0>(m_interpPoints[dir][i]);
903 std::get<2>(m_interpPoints[dir][i]);
906 int tnp = toPointsKey0.GetNumPoints();
907 int nedges = m_interpNfaces[dir][i];
910 switch (m_interpTrace[dir][i])
915 locedges.get() + cnt,
932 locedges.get() + cnt,
943 for (
int k = 0; k < nedges; ++k)
946 &locedges[cnt + k * fnp],
948 &tmp[cnt1 + k * tnp],
952 fnp, locedges.get() + cnt + k * fnp, 1, &I0[0], 1);
958 "Invalid interpolation type for 2D elements");
963 cnt1 += nedges * tnp;
969 m_LocTraceToTraceMap[dir].get(),
981 void LocTraceToTraceMap::InterpLocFacesToTrace(
987 "option dir out of range, " 988 " dir=0 is fwd, dir=1 is bwd");
996 for (
int i = 0; i < m_interpTrace[dir].num_elements(); ++i)
999 if (m_interpNfaces[dir][i])
1003 std::get<0>(m_interpPoints[dir][i]);
1005 std::get<1>(m_interpPoints[dir][i]);
1007 std::get<2>(m_interpPoints[dir][i]);
1009 std::get<3>(m_interpPoints[dir][i]);
1012 int fnp1 = fromPointsKey1.GetNumPoints();
1013 int tnp0 = toPointsKey0.GetNumPoints();
1014 int tnp1 = toPointsKey1.GetNumPoints();
1015 int nfromfacepts = m_interpNfaces[dir][i] * fnp0 * fnp1;
1019 switch (m_interpTrace[dir][i])
1024 locfaces.get() + cnt,
1041 locfaces.get() + cnt,
1050 int nfaces = m_interpNfaces[dir][i];
1051 for (
int k = 0; k < fnp0; ++k)
1054 locfaces.get() + cnt + k,
1056 tmp.get() + cnt1 + k,
1062 tnp1 * m_interpNfaces[dir][i],
1069 tmp.get() + cnt1 + tnp0 - 1,
1076 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1084 locfaces.get() + cnt + j * fnp0 * fnp1,
1089 tmp.get() + cnt1 + j * tnp0 * tnp1,
1097 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1101 locfaces.get() + cnt + j * fnp0 * fnp1,
1103 tmp.get() + cnt1 + j * tnp0 * tnp1,
1107 for (
int k = 0; k < tnp0; ++k)
1109 tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1111 locfaces.get() + cnt +
1112 j * fnp0 * fnp1 + k,
1127 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1135 locfaces.get() + cnt + j * fnp0 * fnp1,
1140 wsp.get() + j * fnp0 * tnp1,
1146 tnp1 * m_interpNfaces[dir][i],
1162 for (
int j = 0; j < m_interpNfaces[dir][i]; ++j)
1170 locfaces.get() + cnt + j * fnp0 * fnp1,
1175 tmp.get() + cnt1 + j * tnp0 * tnp1,
1182 tnp1 * m_interpNfaces[dir][i],
1189 tmp.get() + cnt1 + tnp0 - 1,
1194 cnt += nfromfacepts;
1195 cnt1 += m_interpNfaces[dir][i] * tnp0 * tnp1;
1201 m_LocTraceToTraceMap[dir].get(),
1212 void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1215 int nvals = m_nTraceCoeffs[0] + m_nTraceCoeffs[1];
1216 for (
int i = 0; i < nvals; ++i)
1218 field[m_traceCoeffsToElmtMap[0][i]] +=
1219 m_traceCoeffsToElmtSign[0][i] *
1220 trace[m_traceCoeffsToElmtTrace[0][i]];
1232 void LocTraceToTraceMap::AddTraceCoeffsToFieldCoeffs(
1237 int nvals = m_nTraceCoeffs[dir];
1238 for (
int i = 0; i < nvals; ++i)
1240 field[m_traceCoeffsToElmtMap[dir][i]] +=
1241 m_traceCoeffsToElmtSign[dir][i] *
1242 trace[m_traceCoeffsToElmtTrace[dir][i]];
#define ASSERTL0(condition, msg)
PointsType GetPointsType() const
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
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< DNekMat > DNekMatSharedPtr
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].
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
std::vector< ExpansionSharedPtr > ExpansionVector
Gauss Radau pinned at x=-1, .
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Base class for all multi-elemental spectral/hp expansions.
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.
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].
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
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
unsigned int GetNumPoints() const
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
1D Gauss-Lobatto-Legendre quadrature points