Nektar++
StdExpansion3D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdExpansion3D.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: Daughter of StdExpansion. This class contains routine
33 // which are common to 3D expansion. Typically this inolves physiocal
34 // space operations.
35 //
36 ///////////////////////////////////////////////////////////////////////////////
37 
39 
40 #ifdef max
41 #undef max
42 #endif
43 
44 namespace Nektar
45 {
46  namespace StdRegions
47  {
49  {
50  }
51 
53  const LibUtilities::BasisKey &Ba,
54  const LibUtilities::BasisKey &Bb,
55  const LibUtilities::BasisKey &Bc) :
56  StdExpansion(numcoeffs,3,Ba,Bb,Bc)
57  {
58  }
59 
61  StdExpansion(T)
62  {
63  }
64 
66  {
67  }
69  const Array<OneD, const NekDouble> &inarray,
70  Array<OneD, NekDouble> &out_dx,
71  Array<OneD, NekDouble> &out_dy,
72  Array<OneD, NekDouble> &out_dz)
73  {
74  const int nquad0 = m_base[0]->GetNumPoints();
75  const int nquad1 = m_base[1]->GetNumPoints();
76  const int nquad2 = m_base[2]->GetNumPoints();
77 
78  Array<OneD, NekDouble> wsp(nquad0*nquad1*nquad2);
79 
80  // copy inarray to wsp in case inarray is used as outarray
81  Vmath::Vcopy(nquad0*nquad1*nquad2, &inarray[0], 1, &wsp[0], 1);
82 
83  if (out_dx.num_elements() > 0)
84  {
85  NekDouble *D0 = &((m_base[0]->GetD())->GetPtr())[0];
86 
87  Blas::Dgemm('N','N', nquad0,nquad1*nquad2,nquad0,1.0,
88  D0,nquad0,&wsp[0],nquad0,0.0,&out_dx[0],nquad0);
89  }
90 
91  if (out_dy.num_elements() > 0)
92  {
93  NekDouble *D1 = &((m_base[1]->GetD())->GetPtr())[0];
94  for (int j = 0; j < nquad2; ++j)
95  {
96  Blas::Dgemm('N', 'T', nquad0, nquad1, nquad1,
97  1.0, &wsp[j*nquad0*nquad1], nquad0,
98  D1, nquad1,
99  0.0, &out_dy[j*nquad0*nquad1], nquad0);
100  }
101  }
102 
103  if (out_dz.num_elements() > 0)
104  {
105  NekDouble *D2 = &((m_base[2]->GetD())->GetPtr())[0];
106 
107  Blas::Dgemm('N','T',nquad0*nquad1,nquad2,nquad2,1.0,
108  &wsp[0],nquad0*nquad1,D2,nquad2,0.0,&out_dz[0],
109  nquad0*nquad1);
110  }
111  }
112 
114  const Array<OneD, const NekDouble>& base0,
115  const Array<OneD, const NekDouble>& base1,
116  const Array<OneD, const NekDouble>& base2,
117  const Array<OneD, const NekDouble>& inarray,
118  Array<OneD, NekDouble>& outarray,
120  bool doCheckCollDir0,
121  bool doCheckCollDir1,
122  bool doCheckCollDir2)
123  {
124  v_BwdTrans_SumFacKernel(base0, base1, base2, inarray, outarray, wsp, doCheckCollDir0, doCheckCollDir1, doCheckCollDir2);
125  }
126 
128  const Array<OneD, const NekDouble>& base0,
129  const Array<OneD, const NekDouble>& base1,
130  const Array<OneD, const NekDouble>& base2,
131  const Array<OneD, const NekDouble>& inarray,
132  Array<OneD, NekDouble> &outarray,
134  bool doCheckCollDir0,
135  bool doCheckCollDir1,
136  bool doCheckCollDir2)
137  {
138  v_IProductWRTBase_SumFacKernel(base0, base1, base2, inarray, outarray, wsp, doCheckCollDir0, doCheckCollDir1, doCheckCollDir2);
139  }
140 
142  const Array<OneD, const NekDouble> &coords,
143  const Array<OneD, const NekDouble> &physvals)
144  {
147 
148  WARNINGL2(coords[0] >= -1,"coord[0] < -1");
149  WARNINGL2(coords[0] <= 1,"coord[0] > 1");
150  WARNINGL2(coords[1] >= -1,"coord[1] < -1");
151  WARNINGL2(coords[1] <= 1,"coord[1] > 1");
152  WARNINGL2(coords[2] >= -1,"coord[2] < -1");
153  WARNINGL2(coords[2] <= 1,"coord[2] > 1");
154 
155  // Obtain local collapsed corodinate from
156  // cartesian coordinate.
157  LocCoordToLocCollapsed(coords,eta);
158 
159  // Get Lagrange interpolants.
160  I[0] = m_base[0]->GetI(eta);
161  I[1] = m_base[1]->GetI(eta+1);
162  I[2] = m_base[2]->GetI(eta+2);
163 
164  return v_PhysEvaluate(I,physvals);
165  }
166 
169  const Array<OneD, const NekDouble> &physvals)
170  {
171  NekDouble value;
172 
173  int Qx = m_base[0]->GetNumPoints();
174  int Qy = m_base[1]->GetNumPoints();
175  int Qz = m_base[2]->GetNumPoints();
176 
177  Array<OneD, NekDouble> sumFactorization_qr = Array<OneD, NekDouble>(Qy*Qz);
178  Array<OneD, NekDouble> sumFactorization_r = Array<OneD, NekDouble>(Qz);
179 
180  // Lagrangian interpolation matrix
181  NekDouble *interpolatingNodes = 0;
182 
183  // Interpolate first coordinate direction
184  interpolatingNodes = &I[0]->GetPtr()[0];
185 
186  Blas::Dgemv('T',Qx,Qy*Qz,1.0,&physvals[0],Qx,&interpolatingNodes[0], 1, 0.0, &sumFactorization_qr[0], 1);
187 
188  // Interpolate in second coordinate direction
189  interpolatingNodes = &I[1]->GetPtr()[0];
190 
191  Blas::Dgemv('T',Qy,Qz,1.0,&sumFactorization_qr[0],Qy,&interpolatingNodes[0],1,0.0,&sumFactorization_r[0], 1);
192 
193  // Interpolate in third coordinate direction
194  interpolatingNodes = &I[2]->GetPtr()[0];
195  value = Blas::Ddot(Qz, interpolatingNodes, 1, &sumFactorization_r[0], 1);
196 
197  return value;
198  }
199 
200 
201  /**
202  * @param inarray Input coefficients.
203  * @param output Output coefficients.
204  * @param mkey Matrix key
205  */
207  const Array<OneD, const NekDouble> &inarray,
208  Array<OneD,NekDouble> &outarray,
209  const StdRegions::StdMatrixKey &mkey)
210  {
211  if ( mkey.GetNVarCoeff() == 0 &&
213  {
214  // This implementation is only valid when there are no
215  // coefficients associated to the Laplacian operator
216  int nqtot = GetTotPoints();
217 
218  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
219  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
220  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
221 
222  // Allocate temporary storage
223  Array<OneD,NekDouble> wsp0(7*nqtot);
224  Array<OneD,NekDouble> wsp1(wsp0+nqtot);
225 
226  if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
227  m_base[2]->Collocation()))
228  {
229  // LAPLACIAN MATRIX OPERATION
230  // wsp0 = u = B * u_hat
231  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
232  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
233  BwdTrans_SumFacKernel(base0,base1,base2,inarray,wsp0,wsp1,true,true,true);
234  LaplacianMatrixOp_MatFree_Kernel(wsp0,outarray,wsp1);
235  }
236  else
237  {
238  LaplacianMatrixOp_MatFree_Kernel(inarray,outarray,wsp1);
239  }
240  }
241  else
242  {
244  }
245  }
246 
247 
249  const Array<OneD, const NekDouble> &inarray,
250  Array<OneD,NekDouble> &outarray,
251  const StdRegions::StdMatrixKey &mkey)
252  {
253  if(mkey.GetNVarCoeff() == 0)
254  {
255  int nquad0 = m_base[0]->GetNumPoints();
256  int nquad1 = m_base[1]->GetNumPoints();
257  int nquad2 = m_base[2]->GetNumPoints();
258  int nmodes0 = m_base[0]->GetNumModes();
259  int nmodes1 = m_base[1]->GetNumModes();
260  int nmodes2 = m_base[2]->GetNumModes();
261  int wspsize = max(nquad0*nmodes2*(nmodes1+nquad1),
262  nquad0*nquad1*(nquad2+nmodes0)+
263  nmodes0*nmodes1*nquad2);
264 
266 
267  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
268  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
269  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
270  Array<OneD,NekDouble> wsp0(8*wspsize);
271  Array<OneD,NekDouble> wsp1(wsp0+1*wspsize);
272  Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
273 
274  if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
275  m_base[2]->Collocation()))
276  {
277  // MASS MATRIX OPERATION
278  // The following is being calculated:
279  // wsp0 = B * u_hat = u
280  // wsp1 = W * wsp0
281  // outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
282  BwdTrans_SumFacKernel (base0,base1,base2,inarray,
283  wsp0,wsp2,true,true,true);
284  MultiplyByQuadratureMetric (wsp0,wsp1);
285  IProductWRTBase_SumFacKernel (base0,base1,base2,wsp1,
286  outarray,wsp2,true,true,true);
287  LaplacianMatrixOp_MatFree_Kernel(wsp0,wsp1,wsp2);
288  }
289  else
290  {
291  // specialised implementation for the classical spectral
292  // element method
293  MultiplyByQuadratureMetric (inarray,outarray);
294  LaplacianMatrixOp_MatFree_Kernel(inarray,wsp1,wsp2);
295  }
296 
297  // outarray = lambda * outarray + wsp1
298  // = (lambda * M + L ) * u_hat
299  Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
300  &outarray[0],1);
301  }
302  else
303  {
305  }
306 
307  }
308 
310  const int id) const
311  {
312  return v_GetFaceNormal(id);
313  }
314 
315  const NormalVector & StdExpansion3D::v_GetFaceNormal(const int face) const
316  {
317  std::map<int, NormalVector>::const_iterator x;
318  x = m_faceNormals.find(face);
319  ASSERTL0 (x != m_faceNormals.end(),
320  "face normal not computed.");
321  return x->second;
322  }
323 
325  const Array<OneD, const NekDouble>& inarray)
326  {
327  const int nqtot = GetTotPoints();
329  MultiplyByStdQuadratureMetric(inarray, tmp);
330  return Vmath::Vsum(nqtot, tmp, 1);
331  }
332 
333 
335  {
336  m_negatedNormals[face] = true;
337  for (int i = 0; i < GetCoordim(); ++i)
338  {
339  Vmath::Neg(m_faceNormals[face][i].num_elements(),
340  m_faceNormals[face][i], 1);
341  }
342  }
343 
345  const int facedir,
346  const LibUtilities::BasisType faceDirBasisType,
347  const int numpoints,
348  const int nummodes)
349  {
350 
351  switch(faceDirBasisType)
352  {
354  {
355  const LibUtilities::PointsKey pkey(
357  return LibUtilities::BasisKey(
358  LibUtilities::eModified_A, nummodes, pkey);
359  }
362  {
363  const LibUtilities::PointsKey pkey(
365  return LibUtilities::BasisKey(
366  LibUtilities::eModified_A, nummodes, pkey);
367  }
369  {
370  const LibUtilities::PointsKey pkey(
372  return LibUtilities::BasisKey(
373  LibUtilities::eGLL_Lagrange, nummodes, pkey);
374  }
376  {
377  const LibUtilities::PointsKey pkey(
379  return LibUtilities::BasisKey(
380  LibUtilities::eOrtho_A, nummodes, pkey);
381  }
384  {
385  const LibUtilities::PointsKey pkey(
387  return LibUtilities::BasisKey(
388  LibUtilities::eOrtho_A, nummodes, pkey);
389  }
390  default:
391  {
392  ASSERTL0(false, "expansion type unknown");
393  break;
394  }
395  }
396 
397  // Keep things happy by returning a value.
399  }
400 
402  const int facedir,
403  const LibUtilities::BasisType faceDirBasisType,
404  const int numpoints,
405  const int nummodes)
406  {
407  switch(faceDirBasisType)
408  {
410  {
411  const LibUtilities::PointsKey pkey(
413  return LibUtilities::BasisKey(
414  LibUtilities::eModified_A, nummodes, pkey);
415  }
418  {
419  switch (facedir)
420  {
421  case 0:
422  {
423  const LibUtilities::PointsKey pkey(
424  numpoints+1,
426  return LibUtilities::BasisKey(
427  LibUtilities::eModified_A, nummodes, pkey);
428  }
429  case 1:
430  {
431  const LibUtilities::PointsKey pkey(
432  numpoints+1,
434  return LibUtilities::BasisKey(
435  LibUtilities::eModified_B, nummodes, pkey);
436  }
437  default:
438  {
439  ASSERTL0(false,"invalid value to flag");
440  break;
441  }
442  }
443  }
444 
446  {
447  switch (facedir)
448  {
449  case 0:
450  {
451  const LibUtilities::PointsKey pkey(
452  numpoints,
454  return LibUtilities::BasisKey(
455  LibUtilities::eOrtho_A, nummodes, pkey);
456  }
457  case 1:
458  {
459  const LibUtilities::PointsKey pkey(
460  numpoints,
462  return LibUtilities::BasisKey(
463  LibUtilities::eOrtho_B, nummodes, pkey);
464  }
465  default:
466  {
467  ASSERTL0(false,"invalid value to flag");
468  break;
469  }
470  }
471  }
472 
476  {
477  switch (facedir)
478  {
479  case 0:
480  {
481  const LibUtilities::PointsKey pkey(
482  numpoints,
484  return LibUtilities::BasisKey(
485  LibUtilities::eOrtho_A, nummodes, pkey);
486  }
487  case 1:
488  {
489  const LibUtilities::PointsKey pkey(
490  numpoints,
492  return LibUtilities::BasisKey(
493  LibUtilities::eOrtho_B, nummodes, pkey);
494  }
495  default:
496  {
497  ASSERTL0(false,"invalid value to flag");
498  break;
499  }
500  }
501  }
502  default:
503  {
504  ASSERTL0(false,"expansion type unknown");
505  break;
506  }
507  }
508 
509  // Keep things happy by returning a value.
511  }
512  }//end namespace
513 }//end namespace
void HelmholtzMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void LaplacianMatrixOp_MatFree_GenericImpl(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Principle Modified Functions .
Definition: BasisType.h:51
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:902
const NormalVector & v_GetSurfaceNormal(const int id) const
void MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:909
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
LibUtilities::BasisKey EvaluateTriFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Principle Modified Functions .
Definition: BasisType.h:49
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)=0
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Principle Orthogonal Functions .
Definition: BasisType.h:47
std::map< int, bool > m_negatedNormals
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:131
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points...
The base class for all shapes.
Definition: StdExpansion.h:69
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Principle Modified Functions .
Definition: BasisType.h:50
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Principle Orthogonal Functions .
Definition: BasisType.h:48
Principle Orthogonal Functions .
Definition: BasisType.h:46
Defines a specification for a set of points.
Definition: Points.h:58
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
void LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
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)
Definition: VmathArray.hpp:434
std::map< int, NormalVector > m_faceNormals
const NormalVector & v_GetFaceNormal(const int face) const
virtual void v_NegateFaceNormal(const int face)
#define WARNINGL2(condition, msg)
Definition: ErrorUtil.hpp:214
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Lagrange for SEM basis .
Definition: BasisType.h:53
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:714
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Describes the specification for a Basis.
Definition: Basis.h:50
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)=0