Nektar++
Expansion1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Expansion1D.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: File for Expansion1D routines
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 namespace LocalRegions
45 {
46 const NormalVector &Expansion1D::v_GetTraceNormal(const int vert) const
47 {
48  std::map<int, NormalVector>::const_iterator x;
49  x = m_traceNormals.find(vert);
50  ASSERTL1(x != m_traceNormals.end(), "Vertex normal not computed.");
51  return x->second;
52 }
53 
54 DNekMatSharedPtr Expansion1D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
55 {
56  DNekMatSharedPtr returnval;
57 
58  switch (mkey.GetMatrixType())
59  {
61  {
62  ASSERTL1(IsBoundaryInteriorExpansion(),
63  "HybridDGHelmholtz matrix not set up "
64  "for non boundary-interior expansions");
65  int i;
66  NekDouble lambdaval =
69  int ncoeffs = GetNcoeffs();
70 
71  int coordim = GetCoordim();
72 
73  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
77  DNekMat LocMat(ncoeffs, ncoeffs);
78 
79  returnval =
81  DNekMat &Mat = *returnval;
82 
83  Vmath::Zero(ncoeffs * ncoeffs, Mat.GetPtr(), 1);
84 
85  for (i = 0; i < coordim; ++i)
86  {
87  DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
88 
89  Mat = Mat + Dmat * invMass * Transpose(Dmat);
90  }
91 
92  // Add end Mass Matrix Contribution
93  DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
94  Mat = Mat + lambdaval * Mass;
95 
97  GetBoundaryMap(bmap);
98 
99  // Add tau*F_e using elemental mass matrices
100  for (i = 0; i < 2; ++i)
101  {
102  Mat(bmap[i], bmap[i]) = Mat(bmap[i], bmap[i]) + tau;
103  }
104  }
105  break;
107  {
108  int j, k;
109  int nbndry = NumDGBndryCoeffs();
110  int ncoeffs = GetNcoeffs();
112  factors[StdRegions::eFactorLambda] =
114  factors[StdRegions::eFactorTau] =
116 
117  Array<OneD, NekDouble> lambda(nbndry);
118  DNekVec Lambda(nbndry, lambda, eWrapper);
119  Array<OneD, NekDouble> ulam(ncoeffs);
120  DNekVec Ulam(ncoeffs, ulam, eWrapper);
121  Array<OneD, NekDouble> f(ncoeffs);
122  DNekVec F(ncoeffs, f, eWrapper);
123 
124  // declare matrix space
125  returnval =
127  DNekMat &Umat = *returnval;
128 
129  // Helmholtz matrix
130  DNekScalMat &invHmat =
131  *GetLocMatrix(StdRegions::eInvHybridDGHelmholtz, factors);
132 
133  // for each degree of freedom of the lambda space
134  // calculate Umat entry
135  // Generate Lambda to U_lambda matrix
136  for (j = 0; j < nbndry; ++j)
137  {
138  Vmath::Zero(nbndry, &lambda[0], 1);
139  Vmath::Zero(ncoeffs, &f[0], 1);
140  lambda[j] = 1.0;
141 
142  AddHDGHelmholtzTraceTerms(factors[StdRegions::eFactorTau],
143  lambda, f);
144 
145  Ulam = invHmat * F; // generate Ulam from lambda
146 
147  // fill column of matrix
148  for (k = 0; k < ncoeffs; ++k)
149  {
150  Umat(k, j) = Ulam[k];
151  }
152  }
153  }
154  break;
158  {
159  int j = 0;
160  int k = 0;
161  int dir = 0;
162  int nbndry = NumDGBndryCoeffs();
163  int ncoeffs = GetNcoeffs();
164 
165  Array<OneD, NekDouble> lambda(nbndry);
166  DNekVec Lambda(nbndry, lambda, eWrapper);
167 
168  Array<OneD, NekDouble> ulam(ncoeffs);
169  DNekVec Ulam(ncoeffs, ulam, eWrapper);
170  Array<OneD, NekDouble> f(ncoeffs);
171  DNekVec F(ncoeffs, f, eWrapper);
173  factors[StdRegions::eFactorLambda] =
175  factors[StdRegions::eFactorTau] =
177 
178  // declare matrix space
179  returnval =
181  DNekMat &Qmat = *returnval;
182 
183  // Lambda to U matrix
184  DNekScalMat &lamToU =
185  *GetLocMatrix(StdRegions::eHybridDGLamToU, factors);
186 
187  // Inverse mass matrix
188  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
189 
190  // Weak Derivative matrix
192  switch (mkey.GetMatrixType())
193  {
195  dir = 0;
196  Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
197  break;
199  dir = 1;
200  Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
201  break;
203  dir = 2;
204  Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
205  break;
206  default:
207  ASSERTL0(false, "Direction not known");
208  break;
209  }
210 
211  // for each degree of freedom of the lambda space
212  // calculate Qmat entry
213  // Generate Lambda to Q_lambda matrix
214  for (j = 0; j < nbndry; ++j)
215  {
216  Vmath::Zero(nbndry, &lambda[0], 1);
217  lambda[j] = 1.0;
218 
219  // for lambda[j] = 1 this is the solution to ulam
220  for (k = 0; k < ncoeffs; ++k)
221  {
222  Ulam[k] = lamToU(k, j);
223  }
224 
225  // -D^T ulam
226  Vmath::Neg(ncoeffs, &ulam[0], 1);
227  F = Transpose(*Dmat) * Ulam;
228 
229  // + \tilde{G} \lambda
230  AddNormTraceInt(dir, lambda, f);
231 
232  // multiply by inverse mass matrix
233  Ulam = invMass * F;
234 
235  // fill column of matrix (Qmat is in column major format)
236  Vmath::Vcopy(ncoeffs, &ulam[0], 1,
237  &(Qmat.GetPtr())[0] + j * ncoeffs, 1);
238  }
239  }
240  break;
242  {
243  int j;
244  int nbndry = NumBndryCoeffs();
245 
247  factors[StdRegions::eFactorLambda] =
249  factors[StdRegions::eFactorTau] =
251 
253  Array<OneD, NekDouble> lam(2);
254  GetBoundaryMap(bmap);
255 
256  // declare matrix space
257  returnval =
259  DNekMat &BndMat = *returnval;
260 
261  // Matrix to map Lambda to U
262  DNekScalMat &LamToU =
263  *GetLocMatrix(StdRegions::eHybridDGLamToU, factors);
264 
265  // Matrix to map Lambda to Q
266  DNekScalMat &LamToQ =
267  *GetLocMatrix(StdRegions::eHybridDGLamToQ0, factors);
268 
269  lam[0] = 1.0;
270  lam[1] = 0.0;
271  for (j = 0; j < nbndry; ++j)
272  {
273  BndMat(0, j) =
274  -LamToQ(bmap[0], j) - factors[StdRegions::eFactorTau] *
275  (LamToU(bmap[0], j) - lam[j]);
276  }
277 
278  lam[0] = 0.0;
279  lam[1] = 1.0;
280  for (j = 0; j < nbndry; ++j)
281  {
282  BndMat(1, j) =
283  LamToQ(bmap[1], j) - factors[StdRegions::eFactorTau] *
284  (LamToU(bmap[1], j) - lam[j]);
285  }
286  }
287  break;
288  default:
289  ASSERTL0(false,
290  "This matrix type cannot be generated from this class");
291  break;
292  }
293 
294  return returnval;
295 }
296 
297 void Expansion1D::AddNormTraceInt(const int dir,
299  Array<OneD, NekDouble> &outarray)
300 {
301  boost::ignore_unused(dir);
302 
303  int k;
304  int nbndry = NumBndryCoeffs();
305  int nquad = GetNumPoints(0);
306  const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
308 
309  GetBoundaryMap(vmap);
310 
311  // add G \lambda term (can assume G is diagonal since one
312  // of the basis is zero at boundary otherwise)
313  for (k = 0; k < nbndry; ++k)
314  {
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]) *
318  inarray[vmap[k]];
319  }
320 }
321 
322 void Expansion1D::AddHDGHelmholtzTraceTerms(
323  const NekDouble tau, const Array<OneD, const NekDouble> &inarray,
324  Array<OneD, NekDouble> &outarray)
325 {
326  int i, n;
327  int nbndry = NumBndryCoeffs();
328  int nquad = GetNumPoints(0);
329  int ncoeffs = GetNcoeffs();
330  int coordim = GetCoordim();
332 
333  ASSERTL0(&inarray[0] != &outarray[0],
334  "Input and output arrays use the same memory");
335 
336  const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
337  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
338 
339  GetBoundaryMap(vmap);
340 
341  // Add F = \tau <phi_i,phi_j> (note phi_i is zero if phi_j is non-zero)
342  for (i = 0; i < nbndry; ++i)
343  {
344  outarray[vmap[i]] += tau * Basis[(vmap[i] + 1) * nquad - 1] *
345  Basis[(vmap[i] + 1) * nquad - 1] *
346  inarray[vmap[i]];
347  outarray[vmap[i]] += tau * Basis[vmap[i] * nquad] *
348  Basis[vmap[i] * nquad] * inarray[vmap[i]];
349  }
350 
351  //===============================================================
352  // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
353  // \sum_i D_i M^{-1} G_i term
354 
358  Array<OneD, NekDouble> tmpcoeff(ncoeffs, 0.0);
359  DNekVec Coeffs(ncoeffs, outarray, eWrapper);
360  DNekVec Tmpcoeff(ncoeffs, tmpcoeff, eWrapper);
361 
362  for (n = 0; n < coordim; ++n)
363  {
364  // evaluate M^{-1} G
365  for (i = 0; i < ncoeffs; ++i)
366  {
367  // lower boundary (negative normal)
368  tmpcoeff[i] -= invMass(i, vmap[0]) * Basis[vmap[0] * nquad] *
369  Basis[vmap[0] * nquad] * inarray[vmap[0]];
370 
371  // upper boundary (positive normal)
372  tmpcoeff[i] += invMass(i, vmap[1]) *
373  Basis[(vmap[1] + 1) * nquad - 1] *
374  Basis[(vmap[1] + 1) * nquad - 1] * inarray[vmap[1]];
375  }
376 
377  DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
378  Coeffs = Coeffs + Dmat * Tmpcoeff;
379  }
380 }
381 
382 void Expansion1D::v_AddRobinMassMatrix(
383  const int vert, const Array<OneD, const NekDouble> &primCoeffs,
384  DNekMatSharedPtr &inoutmat)
385 {
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");
391 
392  // Get local Element mapping for vertex point
393  int map = GetVertexMap(vert);
394 
395  // Now need to identify a map which takes the local edge
396  // mass matrix to the matrix stored in inoutmat;
397  // This can currently be deduced from the size of the matrix
398  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
399  // matrix system
400  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
401  // boundary CG system
402 
403  int rows = inoutmat->GetRows();
404 
405  if (rows == GetNcoeffs())
406  {
407  // no need to do anything
408  }
409  else if (rows == NumBndryCoeffs()) // same as NumDGBndryCoeffs()
410  {
411  int i;
413  GetBoundaryMap(bmap);
414 
415  for (i = 0; i < 2; ++i)
416  {
417  if (map == bmap[i])
418  {
419  map = i;
420  break;
421  }
422  }
423  ASSERTL1(i != 2, "Did not find number in map");
424  }
425 
426  // assumes end points have unit magnitude
427  (*inoutmat)(map, map) += primCoeffs[0];
428 }
429 
430 /**
431  * Given an edge and vector of element coefficients:
432  * - maps those elemental coefficients corresponding to the trace into
433  * an vector.
434  * - update the element coefficients
435  * - multiplies the edge vector by the edge mass matrix
436  * - maps the edge coefficients back onto the elemental coefficients
437  */
438 void Expansion1D::v_AddRobinEdgeContribution(
439  const int vert, const Array<OneD, const NekDouble> &primCoeffs,
440  const Array<OneD, NekDouble> &incoeffs, Array<OneD, NekDouble> &coeffs)
441 {
442  ASSERTL1(IsBoundaryInteriorExpansion(),
443  "Not set up for non boundary-interior expansions");
444 
445  int map = GetVertexMap(vert);
446  coeffs[map] += primCoeffs[0] * incoeffs[map];
447 }
448 
449 NekDouble Expansion1D::v_VectorFlux(
450  const Array<OneD, Array<OneD, NekDouble>> &vec)
451 {
453  GetLeftAdjacentElementExp()->GetTraceNormal(
454  GetLeftAdjacentElementTrace());
455 
456  int nq = m_base[0]->GetNumPoints();
457  Array<OneD, NekDouble> Fn(nq);
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);
460 
461  return Integral(Fn);
462 }
463 
464 /** @brief: This method gets all of the factors which are
465  required as part of the Gradient Jump Penalty
466  stabilisation and involves the product of the normal and
467  geometric factors along the element trace.
468 */
469 void Expansion1D::v_NormalTraceDerivFactors(
470  Array<OneD, Array<OneD, NekDouble>> &factors,
471  Array<OneD, Array<OneD, NekDouble>> &d0factors,
472  Array<OneD, Array<OneD, NekDouble>> &d1factors)
473 {
474  boost::ignore_unused(d0factors, d1factors); // for 2D&3D shapes
475  int nquad = GetNumPoints(0);
477  m_metricinfo->GetDerivFactors(GetPointsKeys());
478 
479  if (factors.size() <= 2)
480  {
482  factors[0] = Array<OneD, NekDouble>(1);
483  factors[1] = Array<OneD, NekDouble>(1);
484  }
485 
486  // Outwards normal
488  GetTraceNormal(0);
490  GetTraceNormal(1);
491 
492  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
493  {
494  factors[0][0] = gmat[0][nquad - 1] * normal_0[0][0];
495  factors[1][0] = gmat[0][0] * normal_1[0][0];
496 
497  for (int n = 1; n < normal_0.size(); ++n)
498  {
499  factors[0][0] += gmat[n][0] * normal_0[n][0];
500  factors[1][0] += gmat[n][nquad - 1] * normal_1[n][0];
501  }
502  }
503  else
504  {
505  factors[0][0] = gmat[0][0] * normal_0[0][0];
506  factors[1][0] = gmat[0][0] * normal_1[0][0];
507 
508  for (int n = 1; n < normal_0.size(); ++n)
509  {
510  factors[0][0] += gmat[n][0] * normal_0[n][0];
511  factors[1][0] += gmat[n][0] * normal_1[n][0];
512  }
513  }
514 }
515 
516 void Expansion1D::v_ReOrientTracePhysMap(const StdRegions::Orientation orient,
517  Array<OneD, int> &idmap, const int nq0,
518  const int nq1)
519 {
520  boost::ignore_unused(orient, nq0, nq1);
521 
522  if (idmap.size() != 1)
523  {
524  idmap = Array<OneD, int>(1);
525  }
526 
527  idmap[0] = 0;
528 }
529 
530 void Expansion1D::v_TraceNormLen(const int traceid, NekDouble &h, NekDouble &p)
531 {
532  boost::ignore_unused(traceid);
533  h = GetGeom()->GetVertex(1)->dist(*GetGeom()->GetVertex(0));
534  p = m_ncoeffs - 1;
535 }
536 
537 } // namespace LocalRegions
538 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Represents a basis of a given type.
Definition: Basis.h:213
const Array< OneD, const NekDouble > & GetBdata() const
Return basis definition array m_bdata.
Definition: Basis.h:316
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:126
@ eDeformed
Geometry is curved or has non-constant factors.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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.
Definition: Vmath.cpp:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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
Definition: Vmath.cpp:574
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255