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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: File for Expansion1D routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40  namespace LocalRegions
41  {
42 
44  {
45  DNekMatSharedPtr returnval;
46 
47  switch(mkey.GetMatrixType())
48  {
50  {
52  "HybridDGHelmholtz matrix not set up "
53  "for non boundary-interior expansions");
54  int i;
57  int ncoeffs = GetNcoeffs();
58 
59  int coordim = GetCoordim();
60 
65  DNekMat LocMat(ncoeffs,ncoeffs);
66 
67  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,ncoeffs);
68  DNekMat &Mat = *returnval;
69 
70  Vmath::Zero(ncoeffs*ncoeffs,Mat.GetPtr(),1);
71 
72  for(i=0; i < coordim; ++i)
73  {
74  DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
75 
76  Mat = Mat + Dmat*invMass*Transpose(Dmat);
77  }
78 
79  // Add end Mass Matrix Contribution
81  Mat = Mat + lambdaval*Mass;
82 
84  GetBoundaryMap(bmap);
85 
86  // Add tau*F_e using elemental mass matrices
87  for(i = 0; i < 2; ++i)
88  {
89  Mat(bmap[i],bmap[i]) = Mat(bmap[i],bmap[i]) + tau;
90  }
91  }
92  break;
94  {
95  int j,k;
96  int nbndry = NumDGBndryCoeffs();
97  int ncoeffs = GetNcoeffs();
101 
102  Array<OneD,NekDouble> lambda(nbndry);
103  DNekVec Lambda(nbndry,lambda,eWrapper);
104  Array<OneD,NekDouble> ulam(ncoeffs);
105  DNekVec Ulam(ncoeffs,ulam,eWrapper);
106  Array<OneD,NekDouble> f(ncoeffs);
107  DNekVec F(ncoeffs,f,eWrapper);
108 
109  // declare matrix space
110  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
111  DNekMat &Umat = *returnval;
112 
113  // Helmholtz matrix
115 
116  // for each degree of freedom of the lambda space
117  // calculate Umat entry
118  // Generate Lambda to U_lambda matrix
119  for(j = 0; j < nbndry; ++j)
120  {
121  Vmath::Zero(nbndry,&lambda[0],1);
122  Vmath::Zero(ncoeffs,&f[0],1);
123  lambda[j] = 1.0;
124 
126 
127  Ulam = invHmat*F; // generate Ulam from lambda
128 
129  // fill column of matrix
130  for(k = 0; k < ncoeffs; ++k)
131  {
132  Umat(k,j) = Ulam[k];
133  }
134  }
135  }
136  break;
140  {
141  int j,k,dir;
142  int nbndry = NumDGBndryCoeffs();
143  int ncoeffs = GetNcoeffs();
144 
145  Array<OneD,NekDouble> lambda(nbndry);
146  DNekVec Lambda(nbndry,lambda,eWrapper);
147 
148  Array<OneD,NekDouble> ulam(ncoeffs);
149  DNekVec Ulam(ncoeffs,ulam,eWrapper);
150  Array<OneD,NekDouble> f(ncoeffs);
151  DNekVec F(ncoeffs,f,eWrapper);
155 
156  // declare matrix space
157  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
158  DNekMat &Qmat = *returnval;
159 
160  // Lambda to U matrix
162 
163  // Inverse mass matrix
165 
166  //Weak Derivative matrix
168  switch(mkey.GetMatrixType())
169  {
171  dir = 0;
173  break;
175  dir = 1;
177  break;
179  dir = 2;
181  break;
182  default:
183  ASSERTL0(false,"Direction not known");
184  break;
185  }
186 
187  // for each degree of freedom of the lambda space
188  // calculate Qmat entry
189  // Generate Lambda to Q_lambda matrix
190  for(j = 0; j < nbndry; ++j)
191  {
192  Vmath::Zero(nbndry,&lambda[0],1);
193  lambda[j] = 1.0;
194 
195  // for lambda[j] = 1 this is the solution to ulam
196  for(k = 0; k < ncoeffs; ++k)
197  {
198  Ulam[k] = lamToU(k,j);
199  }
200 
201 
202  // -D^T ulam
203  Vmath::Neg(ncoeffs,&ulam[0],1);
204  F = Transpose(*Dmat)*Ulam;
205 
206  // + \tilde{G} \lambda
207  AddNormTraceInt(dir,lambda,f);
208 
209  // multiply by inverse mass matrix
210  Ulam = invMass*F;
211 
212  // fill column of matrix (Qmat is in column major format)
213  Vmath::Vcopy(ncoeffs,&ulam[0],1,&(Qmat.GetPtr())[0]+j*ncoeffs,1);
214  }
215  }
216  break;
218  {
219  int j;
220  int nbndry = NumBndryCoeffs();
221 
225 
227  Array<OneD, NekDouble> lam(2);
228  GetBoundaryMap(bmap);
229 
230  // declare matrix space
231  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
232  DNekMat &BndMat = *returnval;
233 
234  // Matrix to map Lambda to U
236 
237  // Matrix to map Lambda to Q
239 
240  lam[0] = 1.0; lam[1] = 0.0;
241  for(j = 0; j < nbndry; ++j)
242  {
243  BndMat(0,j) = -LamToQ(bmap[0],j) - factors[StdRegions::eFactorTau]*(LamToU(bmap[0],j) - lam[j]);
244  }
245 
246  lam[0] = 0.0; lam[1] = 1.0;
247  for(j = 0; j < nbndry; ++j)
248  {
249  BndMat(1,j) = LamToQ(bmap[1],j) - factors[StdRegions::eFactorTau]*(LamToU(bmap[1],j) - lam[j]);
250  }
251  }
252  break;
253  default:
254  ASSERTL0(false,"This matrix type cannot be generated from this class");
255  break;
256  }
257 
258  return returnval;
259  }
260 
261  void Expansion1D::AddNormTraceInt(const int dir,
263  Array<OneD,NekDouble> &outarray)
264  {
265 
266  int k;
267  int nbndry = NumBndryCoeffs();
268  int nquad = GetNumPoints(0);
269  const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
271 
272  GetBoundaryMap(vmap);
273 
274  // add G \lambda term (can assume G is diagonal since one
275  // of the basis is zero at boundary otherwise)
276  for(k = 0; k < nbndry; ++k)
277  {
278  outarray[vmap[k]] += (Basis[(vmap[k]+1)*nquad-1]*Basis[(vmap[k]+1)*nquad-1] - Basis[vmap[k]*nquad]*Basis[vmap[k]*nquad])*inarray[vmap[k]];
279  }
280  }
281 
283  const Array<OneD,const NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
284  {
285  int i,n;
286  int nbndry = NumBndryCoeffs();
287  int nquad = GetNumPoints(0);
288  int ncoeffs = GetNcoeffs();
289  int coordim = GetCoordim();
291 
292  ASSERTL0(&inarray[0] != &outarray[0],"Input and output arrays use the same memory");
293 
294 
295  const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
297 
298  GetBoundaryMap(vmap);
299 
300  // Add F = \tau <phi_i,phi_j> (note phi_i is zero if phi_j is non-zero)
301  for(i = 0; i < nbndry; ++i)
302  {
303  outarray[vmap[i]] += tau*Basis[(vmap[i]+1)*nquad-1]*Basis[(vmap[i]+1)*nquad-1]*inarray[vmap[i]];
304  outarray[vmap[i]] += tau*Basis[vmap[i]*nquad]*Basis[vmap[i]*nquad]*inarray[vmap[i]];
305  }
306 
307 
308  //===============================================================
309  // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
310  // \sum_i D_i M^{-1} G_i term
311 
315  Array<OneD, NekDouble> tmpcoeff(ncoeffs,0.0);
316  DNekVec Coeffs (ncoeffs,outarray,eWrapper);
317  DNekVec Tmpcoeff(ncoeffs,tmpcoeff,eWrapper);
318 
319  for(n = 0; n < coordim; ++n)
320  {
321  // evaluate M^{-1} G
322  for(i = 0; i < ncoeffs; ++i)
323  {
324  // lower boundary (negative normal)
325  tmpcoeff[i] -= invMass(i,vmap[0])*Basis[vmap[0]*nquad]*Basis[vmap[0]*nquad]*inarray[vmap[0]];
326 
327  // upper boundary (positive normal)
328  tmpcoeff[i] += invMass(i,vmap[1])*Basis[(vmap[1]+1)*nquad-1]*Basis[(vmap[1]+1)*nquad-1]*inarray[vmap[1]];
329 
330  }
331 
332  DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
333  Coeffs = Coeffs + Dmat*Tmpcoeff;
334  }
335  }
336 
337  void Expansion1D::v_AddRobinMassMatrix(const int vert, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
338  {
339  ASSERTL0(IsBoundaryInteriorExpansion(),"Robin boundary conditions are only implemented for boundary-interior expanisons");
340  ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
341  "Assuming that input matrix was square");
342 
343  // Get local Element mapping for vertex point
344  int map = GetVertexMap(vert);
345 
346  // Now need to identify a map which takes the local edge
347  // mass matrix to the matrix stored in inoutmat;
348  // This can currently be deduced from the size of the matrix
349  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
350  // matrix system
351  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
352  // boundary CG system
353 
354  int rows = inoutmat->GetRows();
355 
356  if (rows == GetNcoeffs())
357  {
358  // no need to do anything
359  }
360  else if(rows == NumBndryCoeffs()) // same as NumDGBndryCoeffs()
361  {
362  int i;
364  GetBoundaryMap(bmap);
365 
366  for(i = 0; i < 2; ++i)
367  {
368  if(map == bmap[i])
369  {
370  map = i;
371  break;
372  }
373  }
374  ASSERTL1(i != 2,"Did not find number in map");
375  }
376 
377  // assumes end points have unit magnitude
378  (*inoutmat)(map,map) += primCoeffs[0];
379 
380  }
381 
382  /**
383  * Given an edge and vector of element coefficients:
384  * - maps those elemental coefficients corresponding to the edge into
385  * an edge-vector.
386  * - resets the element coefficients
387  * - multiplies the edge vector by the edge mass matrix
388  * - maps the edge coefficients back onto the elemental coefficients
389  */
391  {
393  "Not set up for non boundary-interior expansions");
394 
395  int map = GetVertexMap(vert);
396  Vmath::Zero(GetNcoeffs(), coeffs, 1);
397  coeffs[map] = primCoeffs[0];
398  }
399  } //end of namespace
400 } //end of namespace
401 
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
const LibUtilities::BasisSharedPtr & GetBasis(int dir) const
This function gets the shared point to basis in the dir direction.
Definition: StdExpansion.h:118
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
virtual void v_AddRobinMassMatrix(const int vert, const Array< OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion1D.cpp:43
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:229
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:800
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
void AddNormTraceInt(const int dir, Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:83
virtual void v_AddRobinEdgeContribution(const int vert, const Array< OneD, const NekDouble > &primCoeffs, Array< OneD, NekDouble > &coeffs)
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:131
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void AddHDGHelmholtzTraceTerms(const NekDouble tau, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:790