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