35 #include <boost/core/ignore_unused.hpp>
44 namespace LocalRegions
46 const NormalVector &Expansion1D::v_GetTraceNormal(
const int vert)
const
48 std::map<int, NormalVector>::const_iterator x;
49 x = m_traceNormals.find(vert);
50 ASSERTL1(x != m_traceNormals.end(),
"Vertex normal not computed.");
62 ASSERTL1(IsBoundaryInteriorExpansion(),
63 "HybridDGHelmholtz matrix not set up "
64 "for non boundary-interior expansions");
69 int ncoeffs = GetNcoeffs();
71 int coordim = GetCoordim();
77 DNekMat LocMat(ncoeffs, ncoeffs);
85 for (i = 0; i < coordim; ++i)
89 Mat = Mat + Dmat * invMass *
Transpose(Dmat);
94 Mat = Mat + lambdaval * Mass;
100 for (i = 0; i < 2; ++i)
102 Mat(bmap[i], bmap[i]) = Mat(bmap[i], bmap[i]) + tau;
109 int nbndry = NumDGBndryCoeffs();
110 int ncoeffs = GetNcoeffs();
136 for (j = 0; j < nbndry; ++j)
148 for (k = 0; k < ncoeffs; ++k)
150 Umat(k, j) = Ulam[k];
162 int nbndry = NumDGBndryCoeffs();
163 int ncoeffs = GetNcoeffs();
207 ASSERTL0(
false,
"Direction not known");
214 for (j = 0; j < nbndry; ++j)
220 for (k = 0; k < ncoeffs; ++k)
222 Ulam[k] = lamToU(k, j);
230 AddNormTraceInt(dir, lambda, f);
237 &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
244 int nbndry = NumBndryCoeffs();
254 GetBoundaryMap(bmap);
271 for (j = 0; j < nbndry; ++j)
275 (LamToU(bmap[0], j) - lam[j]);
280 for (j = 0; j < nbndry; ++j)
284 (LamToU(bmap[1], j) - lam[j]);
290 "This matrix type cannot be generated from this class");
297 void Expansion1D::AddNormTraceInt(
const int dir,
301 boost::ignore_unused(dir);
304 int nbndry = NumBndryCoeffs();
305 int nquad = GetNumPoints(0);
309 GetBoundaryMap(vmap);
313 for (k = 0; k < nbndry; ++k)
315 outarray[vmap[k]] += (
Basis[(vmap[k] + 1) * nquad - 1] *
316 Basis[(vmap[k] + 1) * nquad - 1] -
317 Basis[vmap[k] * nquad] *
Basis[vmap[k] * nquad]) *
322 void Expansion1D::AddHDGHelmholtzTraceTerms(
327 int nbndry = NumBndryCoeffs();
328 int nquad = GetNumPoints(0);
329 int ncoeffs = GetNcoeffs();
330 int coordim = GetCoordim();
333 ASSERTL0(&inarray[0] != &outarray[0],
334 "Input and output arrays use the same memory");
339 GetBoundaryMap(vmap);
342 for (i = 0; i < nbndry; ++i)
344 outarray[vmap[i]] += tau *
Basis[(vmap[i] + 1) * nquad - 1] *
345 Basis[(vmap[i] + 1) * nquad - 1] *
347 outarray[vmap[i]] += tau *
Basis[vmap[i] * nquad] *
348 Basis[vmap[i] * nquad] * inarray[vmap[i]];
362 for (n = 0; n < coordim; ++n)
365 for (i = 0; i < ncoeffs; ++i)
368 tmpcoeff[i] -= invMass(i, vmap[0]) *
Basis[vmap[0] * nquad] *
369 Basis[vmap[0] * nquad] * inarray[vmap[0]];
372 tmpcoeff[i] += invMass(i, vmap[1]) *
373 Basis[(vmap[1] + 1) * nquad - 1] *
374 Basis[(vmap[1] + 1) * nquad - 1] * inarray[vmap[1]];
378 Coeffs = Coeffs + Dmat * Tmpcoeff;
382 void Expansion1D::v_AddRobinMassMatrix(
386 ASSERTL0(IsBoundaryInteriorExpansion(),
387 "Robin boundary conditions are only implemented for "
388 "boundary-interior expanisons");
389 ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
390 "Assuming that input matrix was square");
393 int map = GetVertexMap(vert);
403 int rows = inoutmat->GetRows();
405 if (rows == GetNcoeffs())
409 else if (rows == NumBndryCoeffs())
413 GetBoundaryMap(bmap);
415 for (i = 0; i < 2; ++i)
423 ASSERTL1(i != 2,
"Did not find number in map");
427 (*inoutmat)(map, map) += primCoeffs[0];
438 void Expansion1D::v_AddRobinEdgeContribution(
442 ASSERTL1(IsBoundaryInteriorExpansion(),
443 "Not set up for non boundary-interior expansions");
445 int map = GetVertexMap(vert);
446 coeffs[map] += primCoeffs[0] * incoeffs[map];
453 GetLeftAdjacentElementExp()->GetTraceNormal(
454 GetLeftAdjacentElementTrace());
456 int nq = m_base[0]->GetNumPoints();
458 Vmath::Vmul(nq, &vec[0][0], 1, &normals[0][0], 1, &Fn[0], 1);
459 Vmath::Vvtvp(nq, &vec[1][0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
469 void Expansion1D::v_NormalTraceDerivFactors(
474 boost::ignore_unused(d0factors, d1factors);
475 int nquad = GetNumPoints(0);
477 m_metricinfo->GetDerivFactors(GetPointsKeys());
479 if (factors.size() <= 2)
494 factors[0][0] = gmat[0][nquad - 1] * normal_0[0][0];
495 factors[1][0] = gmat[0][0] * normal_1[0][0];
497 for (
int n = 1; n < normal_0.size(); ++n)
499 factors[0][0] += gmat[n][0] * normal_0[n][0];
500 factors[1][0] += gmat[n][nquad - 1] * normal_1[n][0];
505 factors[0][0] = gmat[0][0] * normal_0[0][0];
506 factors[1][0] = gmat[0][0] * normal_1[0][0];
508 for (
int n = 1; n < normal_0.size(); ++n)
510 factors[0][0] += gmat[n][0] * normal_0[n][0];
511 factors[1][0] += gmat[n][0] * normal_1[n][0];
520 boost::ignore_unused(orient, nq0, nq1);
522 if (idmap.size() != 1)
532 boost::ignore_unused(traceid);
533 h = GetGeom()->GetVertex(1)->dist(*GetGeom()->GetVertex(0));
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Represents a basis of a given type.
const Array< OneD, const NekDouble > & GetBdata() const
Return basis definition array m_bdata.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
MatrixType GetMatrixType() const
NekDouble GetConstFactor(const ConstFactorType &factor) const
@ eDeformed
Geometry is curved or has non-constant factors.
std::map< ConstFactorType, NekDouble > ConstFactorMap
The above copyright notice and this permission notice shall be included.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)