49 namespace MultiRegions
67 const vector<bool> &LeftAdjacents)
72 int m_expdim = locExpVector[0]->GetShapeDimension();
80 Setup2D(locExp, trace, elmtToTrace, LeftAdjacents);
83 Setup3D(locExp, trace, elmtToTrace, LeftAdjacents);
86 ASSERTL0(
false,
"Number of dimensions greater than 3")
108 const vector<bool> &LeftAdjacents)
122 const boost::shared_ptr<LocalRegions::ExpansionVector> exp =
125 int cnt, n, e, phys_offset;
127 int nexp = exp->size();
138 for (cnt = n = 0; n < nexp; ++n)
142 for (
int i = 0; i < exp2d->GetNedges(); ++i, ++cnt)
144 int nLocPts = exp2d->GetEdgeNumPoints(i);
147 if (LeftAdjacents[cnt])
149 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
150 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
155 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
156 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
167 m_nTraceCoeffs[1] = nBwdCoeffs;
170 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
172 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
174 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
177 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 (it = TraceInterpMap.begin(); it != TraceInterpMap.end(); ++it, ++cnt1)
307 for (
int f = 0; f < it->second.size(); ++f, ++cnt2)
309 n = it->second[f].first;
310 e = it->second[f].second;
319 exp2d->GetEdgePhysMap(e, edgeids);
321 nedgepts = exp2d->GetEdgeNumPoints(e);
322 nedgepts1 = edge->GetTotPoints();
327 exp2d->ReOrientEdgePhysMap(elmtToTrace[n][e]->GetNverts(),
332 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
334 if (LeftAdjacents[TraceOrder[n][e]])
336 for (
int i = 0; i < nedgepts; ++i)
341 for (
int i = 0; i < nedgepts1; ++i)
343 m_LocTraceToTraceMap[0][cntFwd1 + i] =
344 offset + locTraceToTraceMap[i];
348 cntFwd1 += nedgepts1;
353 for (
int i = 0; i < nedgepts; ++i)
356 phys_offset + edgeids[i];
359 for (
int i = 0; i < nedgepts1; ++i)
361 m_LocTraceToTraceMap[1][cntBwd1 + i] =
362 offset + locTraceToTraceMap[i];
366 cntBwd1 += nedgepts1;
370 m_interpNfaces[set][cnt1] += 1;
372 if ((fwdSet ==
false && set == 0) || (bwdSet ==
false && set == 1))
374 m_interpPoints[set][cnt1] = it->first;
376 if (fromPointsKey0 == toPointsKey0)
383 m_interpTraceI0[set][cnt1] =
403 m_interpEndPtI0[set][cnt1] =
408 m_interpTraceI0[set][cnt1]->GetPtr().
get() +
411 &m_interpEndPtI0[set][cnt1][0],
443 const vector<bool> &LeftAdjacents)
460 const boost::shared_ptr<LocalRegions::ExpansionVector> exp =
463 int cnt, n, e, phys_offset;
465 int nexp = exp->size();
476 for (cnt = n = 0; n < nexp; ++n)
480 for (
int i = 0; i < exp3d->GetNfaces(); ++i, ++cnt)
482 int nLocPts = exp3d->GetFaceNumPoints(i);
485 if (LeftAdjacents[cnt])
487 nFwdPts += elmtToTrace[n][i]->GetTotPoints();
488 nFwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
493 nBwdPts += elmtToTrace[n][i]->GetTotPoints();
494 nBwdCoeffs += elmtToTrace[n][i]->GetNcoeffs();
505 m_nTraceCoeffs[1] = nBwdCoeffs;
508 m_traceCoeffsToElmtMap[1] = m_traceCoeffsToElmtMap[0] + nFwdCoeffs;
510 m_traceCoeffsToElmtTrace[1] = m_traceCoeffsToElmtTrace[0] + nFwdCoeffs;
512 m_traceCoeffsToElmtSign[1] = m_traceCoeffsToElmtSign[0] + nFwdCoeffs;
515 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 (it = TraceInterpMap.begin(); it != TraceInterpMap.end(); ++it, ++cnt1)
647 for (
int f = 0; f < it->second.size(); ++f, ++cnt2)
649 n = it->second[f].first;
650 e = it->second[f].second;
660 exp3d->GetFacePhysMap(e, faceids);
661 nfacepts = exp3d->GetFaceNumPoints(e);
662 nfacepts1 = face->GetTotPoints();
666 exp3d->ReOrientFacePhysMap(elmtToTrace[n][e]->GetNverts(),
672 int offset = trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
674 if (LeftAdjacents[TraceOrder[n][e]])
676 for (
int i = 0; i < nfacepts; ++i)
681 for (
int i = 0; i < nfacepts1; ++i)
683 m_LocTraceToTraceMap[0][cntFwd1 + i] =
684 offset + locTraceToTraceMap[i];
688 cntFwd1 += nfacepts1;
693 for (
int i = 0; i < nfacepts; ++i)
696 phys_offset + faceids[i];
699 for (
int i = 0; i < nfacepts1; ++i)
701 m_LocTraceToTraceMap[1][cntBwd1 + i] =
702 offset + locTraceToTraceMap[i];
706 cntBwd1 += nfacepts1;
710 m_interpNfaces[set][cnt1] += 1;
712 if (((fwdSet ==
false) && (set == 0)) ||
713 ((bwdSet ==
false) && (set == 1)))
715 m_interpPoints[set][cnt1] = it->first;
717 if (fromPointsKey0 == toPointsKey0)
719 if (fromPointsKey1 == toPointsKey1)
726 m_interpTraceI1[set][cnt1] =
743 m_interpEndPtI1[set][cnt1] =
747 m_interpTraceI1[set][cnt1]->GetPtr().
get() +
750 &m_interpEndPtI1[set][cnt1][0],
758 if (fromPointsKey1 == toPointsKey1)
761 m_interpTraceI0[set][cnt1] =
778 m_interpEndPtI0[set][cnt1] =
782 m_interpTraceI0[set][cnt1]->GetPtr().
get() +
785 &m_interpEndPtI0[set][cnt1][0],
793 m_interpTraceI0[set][cnt1] =
796 m_interpTraceI1[set][cnt1] =
810 m_interpTrace[set][cnt1] =
814 m_interpEndPtI0[set][cnt1] =
818 m_interpTraceI0[set][cnt1]->GetPtr().
get() +
821 &m_interpEndPtI0[set][cnt1][0],
884 "option dir out of range, "
885 " dir=0 is fwd, dir=1 is bwd");
914 locedges.get() + cnt,
931 locedges.get() + cnt,
942 for (
int k = 0; k < nedges; ++k)
945 &locedges[cnt + k * fnp],
947 &tmp[cnt1 + k * tnp],
951 fnp, locedges.get() + cnt + k * fnp, 1, &I0[0], 1);
957 "Invalid interpolation type for 2D elements");
962 cnt1 += nedges * tnp;
986 "option dir out of range, "
987 " dir=0 is fwd, dir=1 is bwd");
1023 locfaces.get() + cnt,
1040 locfaces.get() + cnt,
1050 for (
int k = 0; k < fnp0; ++k)
1053 locfaces.get() + cnt + k,
1055 tmp.get() + cnt1 + k,
1068 tmp.get() + cnt1 + tnp0 - 1,
1083 locfaces.get() + cnt + j * fnp0 * fnp1,
1088 tmp.get() + cnt1 + j * tnp0 * tnp1,
1100 locfaces.get() + cnt + j * fnp0 * fnp1,
1102 tmp.get() + cnt1 + j * tnp0 * tnp1,
1106 for (
int k = 0; k < tnp0; ++k)
1108 tmp[cnt1 + k + (j + 1) * tnp0 * tnp1 - tnp0] =
1110 locfaces.get() + cnt +
1111 j * fnp0 * fnp1 + k,
1134 locfaces.get() + cnt + j * fnp0 * fnp1,
1139 wsp.get() + j * fnp0 * tnp1,
1145 tnp1 * m_interpNfaces[dir][i],
1169 locfaces.get() + cnt + j * fnp0 * fnp1,
1174 tmp.get() + cnt1 + j * tnp0 * tnp1,
1181 tnp1 * m_interpNfaces[dir][i],
1188 tmp.get() + cnt1 + tnp0 - 1,
1193 cnt += nfromfacepts;
1215 for (
int i = 0; i < nvals; ++i)
1237 for (
int i = 0; i < nvals; ++i)
virtual ~LocTraceToTraceMap()
int m_nTracePts
The number of global trace points.
#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]].
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
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
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.
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
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
Gauss Radau pinned at x=-1, .
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.
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
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].
Defines a specification for a set of points.
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)
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
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...
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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
PointsType GetPointsType() const
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