Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
39 
40 namespace Nektar
41 {
42  namespace LocalRegions
43  {
44 
45  DNekMatSharedPtr Expansion1D::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
46  {
47  DNekMatSharedPtr returnval;
48 
49  switch(mkey.GetMatrixType())
50  {
52  {
53  ASSERTL1(IsBoundaryInteriorExpansion(),
54  "HybridDGHelmholtz matrix not set up "
55  "for non boundary-interior expansions");
56  int i;
59  int ncoeffs = GetNcoeffs();
60 
61  int coordim = GetCoordim();
62 
63  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
67  DNekMat LocMat(ncoeffs,ncoeffs);
68 
69  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,ncoeffs);
70  DNekMat &Mat = *returnval;
71 
72  Vmath::Zero(ncoeffs*ncoeffs,Mat.GetPtr(),1);
73 
74  for(i=0; i < coordim; ++i)
75  {
76  DNekScalMat &Dmat = *GetLocMatrix(DerivType[i]);
77 
78  Mat = Mat + Dmat*invMass*Transpose(Dmat);
79  }
80 
81  // Add end Mass Matrix Contribution
82  DNekScalMat &Mass = *GetLocMatrix(StdRegions::eMass);
83  Mat = Mat + lambdaval*Mass;
84 
86  GetBoundaryMap(bmap);
87 
88  // Add tau*F_e using elemental mass matrices
89  for(i = 0; i < 2; ++i)
90  {
91  Mat(bmap[i],bmap[i]) = Mat(bmap[i],bmap[i]) + tau;
92  }
93  }
94  break;
96  {
97  int j,k;
98  int nbndry = NumDGBndryCoeffs();
99  int ncoeffs = GetNcoeffs();
103 
104  Array<OneD,NekDouble> lambda(nbndry);
105  DNekVec Lambda(nbndry,lambda,eWrapper);
106  Array<OneD,NekDouble> ulam(ncoeffs);
107  DNekVec Ulam(ncoeffs,ulam,eWrapper);
108  Array<OneD,NekDouble> f(ncoeffs);
109  DNekVec F(ncoeffs,f,eWrapper);
110 
111  // declare matrix space
112  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
113  DNekMat &Umat = *returnval;
114 
115  // Helmholtz matrix
116  DNekScalMat &invHmat = *GetLocMatrix(StdRegions::eInvHybridDGHelmholtz, factors);
117 
118  // for each degree of freedom of the lambda space
119  // calculate Umat entry
120  // Generate Lambda to U_lambda matrix
121  for(j = 0; j < nbndry; ++j)
122  {
123  Vmath::Zero(nbndry,&lambda[0],1);
124  Vmath::Zero(ncoeffs,&f[0],1);
125  lambda[j] = 1.0;
126 
127  AddHDGHelmholtzTraceTerms(factors[StdRegions::eFactorTau],lambda,f);
128 
129  Ulam = invHmat*F; // generate Ulam from lambda
130 
131  // fill column of matrix
132  for(k = 0; k < ncoeffs; ++k)
133  {
134  Umat(k,j) = Ulam[k];
135  }
136  }
137  }
138  break;
142  {
143  int j,k,dir;
144  int nbndry = NumDGBndryCoeffs();
145  int ncoeffs = GetNcoeffs();
146 
147  Array<OneD,NekDouble> lambda(nbndry);
148  DNekVec Lambda(nbndry,lambda,eWrapper);
149 
150  Array<OneD,NekDouble> ulam(ncoeffs);
151  DNekVec Ulam(ncoeffs,ulam,eWrapper);
152  Array<OneD,NekDouble> f(ncoeffs);
153  DNekVec F(ncoeffs,f,eWrapper);
157 
158  // declare matrix space
159  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(ncoeffs,nbndry);
160  DNekMat &Qmat = *returnval;
161 
162  // Lambda to U matrix
163  DNekScalMat &lamToU = *GetLocMatrix(StdRegions::eHybridDGLamToU, factors);
164 
165  // Inverse mass matrix
166  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
167 
168  //Weak Derivative matrix
170  switch(mkey.GetMatrixType())
171  {
173  dir = 0;
174  Dmat = GetLocMatrix(StdRegions::eWeakDeriv0);
175  break;
177  dir = 1;
178  Dmat = GetLocMatrix(StdRegions::eWeakDeriv1);
179  break;
181  dir = 2;
182  Dmat = GetLocMatrix(StdRegions::eWeakDeriv2);
183  break;
184  default:
185  ASSERTL0(false,"Direction not known");
186  break;
187  }
188 
189  // for each degree of freedom of the lambda space
190  // calculate Qmat entry
191  // Generate Lambda to Q_lambda matrix
192  for(j = 0; j < nbndry; ++j)
193  {
194  Vmath::Zero(nbndry,&lambda[0],1);
195  lambda[j] = 1.0;
196 
197  // for lambda[j] = 1 this is the solution to ulam
198  for(k = 0; k < ncoeffs; ++k)
199  {
200  Ulam[k] = lamToU(k,j);
201  }
202 
203 
204  // -D^T ulam
205  Vmath::Neg(ncoeffs,&ulam[0],1);
206  F = Transpose(*Dmat)*Ulam;
207 
208  // + \tilde{G} \lambda
209  AddNormTraceInt(dir,lambda,f);
210 
211  // multiply by inverse mass matrix
212  Ulam = invMass*F;
213 
214  // fill column of matrix (Qmat is in column major format)
215  Vmath::Vcopy(ncoeffs,&ulam[0],1,&(Qmat.GetPtr())[0]+j*ncoeffs,1);
216  }
217  }
218  break;
220  {
221  int j;
222  int nbndry = NumBndryCoeffs();
223 
227 
229  Array<OneD, NekDouble> lam(2);
230  GetBoundaryMap(bmap);
231 
232  // declare matrix space
233  returnval = MemoryManager<DNekMat>::AllocateSharedPtr(nbndry, nbndry);
234  DNekMat &BndMat = *returnval;
235 
236  // Matrix to map Lambda to U
237  DNekScalMat &LamToU = *GetLocMatrix(StdRegions::eHybridDGLamToU, factors);
238 
239  // Matrix to map Lambda to Q
240  DNekScalMat &LamToQ = *GetLocMatrix(StdRegions::eHybridDGLamToQ0, factors);
241 
242  lam[0] = 1.0; lam[1] = 0.0;
243  for(j = 0; j < nbndry; ++j)
244  {
245  BndMat(0,j) = -LamToQ(bmap[0],j) - factors[StdRegions::eFactorTau]*(LamToU(bmap[0],j) - lam[j]);
246  }
247 
248  lam[0] = 0.0; lam[1] = 1.0;
249  for(j = 0; j < nbndry; ++j)
250  {
251  BndMat(1,j) = LamToQ(bmap[1],j) - factors[StdRegions::eFactorTau]*(LamToU(bmap[1],j) - lam[j]);
252  }
253  }
254  break;
255  default:
256  ASSERTL0(false,"This matrix type cannot be generated from this class");
257  break;
258  }
259 
260  return returnval;
261  }
262 
263  void Expansion1D::AddNormTraceInt(const int dir,
265  Array<OneD,NekDouble> &outarray)
266  {
267 
268  int k;
269  int nbndry = NumBndryCoeffs();
270  int nquad = GetNumPoints(0);
271  const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
273 
274  GetBoundaryMap(vmap);
275 
276  // add G \lambda term (can assume G is diagonal since one
277  // of the basis is zero at boundary otherwise)
278  for(k = 0; k < nbndry; ++k)
279  {
280  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]];
281  }
282  }
283 
284  void Expansion1D::AddHDGHelmholtzTraceTerms(const NekDouble tau,
285  const Array<OneD,const NekDouble> &inarray, Array<OneD,NekDouble> &outarray)
286  {
287  int i,n;
288  int nbndry = NumBndryCoeffs();
289  int nquad = GetNumPoints(0);
290  int ncoeffs = GetNcoeffs();
291  int coordim = GetCoordim();
293 
294  ASSERTL0(&inarray[0] != &outarray[0],"Input and output arrays use the same memory");
295 
296 
297  const Array<OneD, const NekDouble> &Basis = GetBasis(0)->GetBdata();
298  DNekScalMat &invMass = *GetLocMatrix(StdRegions::eInvMass);
299 
300  GetBoundaryMap(vmap);
301 
302  // Add F = \tau <phi_i,phi_j> (note phi_i is zero if phi_j is non-zero)
303  for(i = 0; i < nbndry; ++i)
304  {
305  outarray[vmap[i]] += tau*Basis[(vmap[i]+1)*nquad-1]*Basis[(vmap[i]+1)*nquad-1]*inarray[vmap[i]];
306  outarray[vmap[i]] += tau*Basis[vmap[i]*nquad]*Basis[vmap[i]*nquad]*inarray[vmap[i]];
307  }
308 
309 
310  //===============================================================
311  // Add -\sum_i D_i^T M^{-1} G_i + E_i M^{-1} G_i =
312  // \sum_i D_i M^{-1} G_i term
313 
317  Array<OneD, NekDouble> tmpcoeff(ncoeffs,0.0);
318  DNekVec Coeffs (ncoeffs,outarray,eWrapper);
319  DNekVec Tmpcoeff(ncoeffs,tmpcoeff,eWrapper);
320 
321  for(n = 0; n < coordim; ++n)
322  {
323  // evaluate M^{-1} G
324  for(i = 0; i < ncoeffs; ++i)
325  {
326  // lower boundary (negative normal)
327  tmpcoeff[i] -= invMass(i,vmap[0])*Basis[vmap[0]*nquad]*Basis[vmap[0]*nquad]*inarray[vmap[0]];
328 
329  // upper boundary (positive normal)
330  tmpcoeff[i] += invMass(i,vmap[1])*Basis[(vmap[1]+1)*nquad-1]*Basis[(vmap[1]+1)*nquad-1]*inarray[vmap[1]];
331 
332  }
333 
334  DNekScalMat &Dmat = *GetLocMatrix(DerivType[n]);
335  Coeffs = Coeffs + Dmat*Tmpcoeff;
336  }
337  }
338 
339  void Expansion1D::v_AddRobinMassMatrix(const int vert, const Array<OneD, const NekDouble > &primCoeffs, DNekMatSharedPtr &inoutmat)
340  {
341  ASSERTL0(IsBoundaryInteriorExpansion(),"Robin boundary conditions are only implemented for boundary-interior expanisons");
342  ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
343  "Assuming that input matrix was square");
344 
345  // Get local Element mapping for vertex point
346  int map = GetVertexMap(vert);
347 
348  // Now need to identify a map which takes the local edge
349  // mass matrix to the matrix stored in inoutmat;
350  // This can currently be deduced from the size of the matrix
351  // - if inoutmat.m_rows() == v_NCoeffs() it is a full
352  // matrix system
353  // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
354  // boundary CG system
355 
356  int rows = inoutmat->GetRows();
357 
358  if (rows == GetNcoeffs())
359  {
360  // no need to do anything
361  }
362  else if(rows == NumBndryCoeffs()) // same as NumDGBndryCoeffs()
363  {
364  int i;
366  GetBoundaryMap(bmap);
367 
368  for(i = 0; i < 2; ++i)
369  {
370  if(map == bmap[i])
371  {
372  map = i;
373  break;
374  }
375  }
376  ASSERTL1(i != 2,"Did not find number in map");
377  }
378 
379  // assumes end points have unit magnitude
380  (*inoutmat)(map,map) += primCoeffs[0];
381 
382  }
383 
384  /**
385  * Given an edge and vector of element coefficients:
386  * - maps those elemental coefficients corresponding to the edge into
387  * an edge-vector.
388  * - resets the element coefficients
389  * - multiplies the edge vector by the edge mass matrix
390  * - maps the edge coefficients back onto the elemental coefficients
391  */
392  void Expansion1D::v_AddRobinEdgeContribution(const int vert, const Array<OneD, const NekDouble > &primCoeffs, Array<OneD, NekDouble> &coeffs)
393  {
394  ASSERTL1(IsBoundaryInteriorExpansion(),
395  "Not set up for non boundary-interior expansions");
396 
397  int map = GetVertexMap(vert);
398  Vmath::Zero(GetNcoeffs(), coeffs, 1);
399  coeffs[map] = primCoeffs[0];
400  }
401  } //end of namespace
402 } //end of namespace
403 
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:251
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047