Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  using std::max;
256 
257  int nquad0 = m_base[0]->GetNumPoints();
258  int nquad1 = m_base[1]->GetNumPoints();
259  int nquad2 = m_base[2]->GetNumPoints();
260  int nmodes0 = m_base[0]->GetNumModes();
261  int nmodes1 = m_base[1]->GetNumModes();
262  int nmodes2 = m_base[2]->GetNumModes();
263  int wspsize = max(nquad0*nmodes2*(nmodes1+nquad1),
264  nquad0*nquad1*(nquad2+nmodes0)+
265  nmodes0*nmodes1*nquad2);
266 
268 
269  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
270  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
271  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
272  Array<OneD,NekDouble> wsp0(8*wspsize);
273  Array<OneD,NekDouble> wsp1(wsp0+1*wspsize);
274  Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
275 
276  if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
277  m_base[2]->Collocation()))
278  {
279  // MASS MATRIX OPERATION
280  // The following is being calculated:
281  // wsp0 = B * u_hat = u
282  // wsp1 = W * wsp0
283  // outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
284  BwdTrans_SumFacKernel (base0,base1,base2,inarray,
285  wsp0,wsp2,true,true,true);
286  MultiplyByQuadratureMetric (wsp0,wsp1);
287  IProductWRTBase_SumFacKernel (base0,base1,base2,wsp1,
288  outarray,wsp2,true,true,true);
289  LaplacianMatrixOp_MatFree_Kernel(wsp0,wsp1,wsp2);
290  }
291  else
292  {
293  // specialised implementation for the classical spectral
294  // element method
295  MultiplyByQuadratureMetric (inarray,outarray);
296  LaplacianMatrixOp_MatFree_Kernel(inarray,wsp1,wsp2);
297  }
298 
299  // outarray = lambda * outarray + wsp1
300  // = (lambda * M + L ) * u_hat
301  Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
302  &outarray[0],1);
303  }
304  else
305  {
307  }
308 
309  }
310 
312  const int id) const
313  {
314  return v_GetFaceNormal(id);
315  }
316 
317  const NormalVector & StdExpansion3D::v_GetFaceNormal(const int face) const
318  {
319  std::map<int, NormalVector>::const_iterator x;
320  x = m_faceNormals.find(face);
321  ASSERTL0 (x != m_faceNormals.end(),
322  "face normal not computed.");
323  return x->second;
324  }
325 
327  const Array<OneD, const NekDouble>& inarray)
328  {
329  const int nqtot = GetTotPoints();
331  MultiplyByStdQuadratureMetric(inarray, tmp);
332  return Vmath::Vsum(nqtot, tmp, 1);
333  }
334 
335 
337  {
338  m_negatedNormals[face] = true;
339  for (int i = 0; i < GetCoordim(); ++i)
340  {
341  Vmath::Neg(m_faceNormals[face][i].num_elements(),
342  m_faceNormals[face][i], 1);
343  }
344  }
345 
347  {
348  return m_negatedNormals[face];
349  }
350 
352  const int facedir,
353  const LibUtilities::BasisType faceDirBasisType,
354  const int numpoints,
355  const int nummodes)
356  {
357 
358  switch(faceDirBasisType)
359  {
361  {
362  const LibUtilities::PointsKey pkey(
364  return LibUtilities::BasisKey(
365  LibUtilities::eModified_A, nummodes, pkey);
366  }
369  {
370  const LibUtilities::PointsKey pkey(
372  return LibUtilities::BasisKey(
373  LibUtilities::eModified_A, nummodes, pkey);
374  }
376  {
377  const LibUtilities::PointsKey pkey(
379  return LibUtilities::BasisKey(
380  LibUtilities::eGLL_Lagrange, nummodes, pkey);
381  }
383  {
384  const LibUtilities::PointsKey pkey(
386  return LibUtilities::BasisKey(
387  LibUtilities::eOrtho_A, nummodes, pkey);
388  }
391  {
392  const LibUtilities::PointsKey pkey(
394  return LibUtilities::BasisKey(
395  LibUtilities::eOrtho_A, nummodes, pkey);
396  }
397  default:
398  {
399  ASSERTL0(false, "expansion type unknown");
400  break;
401  }
402  }
403 
404  // Keep things happy by returning a value.
406  }
407 
409  const int facedir,
410  const LibUtilities::BasisType faceDirBasisType,
411  const int numpoints,
412  const int nummodes)
413  {
414  switch(faceDirBasisType)
415  {
417  {
418  const LibUtilities::PointsKey pkey(
420  return LibUtilities::BasisKey(
421  LibUtilities::eModified_A, nummodes, pkey);
422  }
425  {
426  switch (facedir)
427  {
428  case 0:
429  {
430  const LibUtilities::PointsKey pkey(
431  numpoints+1,
433  return LibUtilities::BasisKey(
434  LibUtilities::eModified_A, nummodes, pkey);
435  }
436  case 1:
437  {
438 // const LibUtilities::PointsKey pkey(
439 // numpoints+1,
440 // LibUtilities::eGaussLobattoLegendre);
441  const LibUtilities::PointsKey pkey(
442  numpoints,
444  return LibUtilities::BasisKey(
445  LibUtilities::eModified_B, nummodes, pkey);
446  }
447  default:
448  {
449 
450  ASSERTL0(false,"invalid value to flag");
451  break;
452  }
453  }
454  }
455 
457  {
458  switch (facedir)
459  {
460  case 0:
461  {
462  const LibUtilities::PointsKey pkey(
463  numpoints,
465  return LibUtilities::BasisKey(
466  LibUtilities::eOrtho_A, nummodes, pkey);
467  }
468  case 1:
469  {
470  const LibUtilities::PointsKey pkey(
471  numpoints,
473  return LibUtilities::BasisKey(
474  LibUtilities::eOrtho_B, nummodes, pkey);
475  }
476  default:
477  {
478  ASSERTL0(false,"invalid value to flag");
479  break;
480  }
481  }
482  }
483 
487  {
488  switch (facedir)
489  {
490  case 0:
491  {
492  const LibUtilities::PointsKey pkey(
493  numpoints,
495  return LibUtilities::BasisKey(
496  LibUtilities::eOrtho_A, nummodes, pkey);
497  }
498  case 1:
499  {
500  const LibUtilities::PointsKey pkey(
501  numpoints,
503  return LibUtilities::BasisKey(
504  LibUtilities::eOrtho_B, nummodes, pkey);
505  }
506  default:
507  {
508  ASSERTL0(false,"invalid value to flag");
509  break;
510  }
511  }
512  }
513  default:
514  {
515  ASSERTL0(false,"expansion type unknown");
516  break;
517  }
518  }
519 
520  // Keep things happy by returning a value.
522  }
523  }//end namespace
524 }//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:942
const NormalVector & v_GetSurfaceNormal(const int id) const
void MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:949
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 bool v_FaceNormalNegated(const int face)
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:723
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:1047
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