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 
38 #include <StdRegions/StdExpansion3D.h>
39 
40 #ifdef max
41 #undef max
42 #endif
43 
44 namespace Nektar
45 {
46  namespace StdRegions
47  {
48  StdExpansion3D::StdExpansion3D()
49  {
50  }
51 
52  StdExpansion3D::StdExpansion3D(int numcoeffs,
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 
60  StdExpansion3D::StdExpansion3D(const StdExpansion3D &T):
61  StdExpansion(T)
62  {
63  }
64 
65  StdExpansion3D::~StdExpansion3D()
66  {
67  }
68  void StdExpansion3D::PhysTensorDeriv(
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 
113  void StdExpansion3D::BwdTrans_SumFacKernel(
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,
119  Array<OneD, NekDouble>& wsp,
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 
127  void StdExpansion3D::IProductWRTBase_SumFacKernel(
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,
133  Array<OneD, NekDouble> &wsp,
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 
141  NekDouble StdExpansion3D::v_PhysEvaluate(
142  const Array<OneD, const NekDouble> &coords,
143  const Array<OneD, const NekDouble> &physvals)
144  {
145  Array<OneD, NekDouble> eta = Array<OneD, NekDouble>(3);
146  Array<OneD, DNekMatSharedPtr> I(3);
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 
167  NekDouble StdExpansion3D::v_PhysEvaluate(
168  const Array<OneD, DNekMatSharedPtr > &I,
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  */
206  void StdExpansion3D::v_LaplacianMatrixOp_MatFree(
207  const Array<OneD, const NekDouble> &inarray,
208  Array<OneD,NekDouble> &outarray,
209  const StdRegions::StdMatrixKey &mkey)
210  {
211  if ( mkey.GetNVarCoeff() == 0 &&
212  !mkey.ConstFactorExists(eFactorSVVCutoffRatio))
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  {
243  StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
244  }
245  }
246 
247 
248  void StdExpansion3D::v_HelmholtzMatrixOp_MatFree(
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 
265  NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
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  {
304  StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
305  }
306 
307  }
308 
309  const NormalVector & StdExpansion3D::v_GetSurfaceNormal(
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 
324  NekDouble StdExpansion3D::v_Integral(
325  const Array<OneD, const NekDouble>& inarray)
326  {
327  const int nqtot = GetTotPoints();
328  Array<OneD, NekDouble> tmp(GetTotPoints());
329  MultiplyByStdQuadratureMetric(inarray, tmp);
330  return Vmath::Vsum(nqtot, tmp, 1);
331  }
332 
333 
334  void StdExpansion3D::v_NegateFaceNormal(const int face)
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 
344  LibUtilities::BasisKey EvaluateQuadFaceBasisKey(
345  const int facedir,
346  const LibUtilities::BasisType faceDirBasisType,
347  const int numpoints,
348  const int nummodes)
349  {
350 
351  switch(faceDirBasisType)
352  {
353  case LibUtilities::eModified_A:
354  {
355  const LibUtilities::PointsKey pkey(
356  numpoints, LibUtilities::eGaussLobattoLegendre);
357  return LibUtilities::BasisKey(
358  LibUtilities::eModified_A, nummodes, pkey);
359  }
360  case LibUtilities::eModified_B:
361  case LibUtilities::eModified_C:
362  {
363  const LibUtilities::PointsKey pkey(
364  numpoints+1, LibUtilities::eGaussLobattoLegendre);
365  return LibUtilities::BasisKey(
366  LibUtilities::eModified_A, nummodes, pkey);
367  }
368  case LibUtilities::eGLL_Lagrange:
369  {
370  const LibUtilities::PointsKey pkey(
371  numpoints, LibUtilities::eGaussLobattoLegendre);
372  return LibUtilities::BasisKey(
373  LibUtilities::eGLL_Lagrange, nummodes, pkey);
374  }
375  case LibUtilities::eOrtho_A:
376  {
377  const LibUtilities::PointsKey pkey(
378  numpoints, LibUtilities::eGaussLobattoLegendre);
379  return LibUtilities::BasisKey(
380  LibUtilities::eOrtho_A, nummodes, pkey);
381  }
382  case LibUtilities::eOrtho_B:
383  case LibUtilities::eOrtho_C:
384  {
385  const LibUtilities::PointsKey pkey(
386  numpoints+1, LibUtilities::eGaussLobattoLegendre);
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.
398  return LibUtilities::NullBasisKey;
399  }
400 
401  LibUtilities::BasisKey EvaluateTriFaceBasisKey(
402  const int facedir,
403  const LibUtilities::BasisType faceDirBasisType,
404  const int numpoints,
405  const int nummodes)
406  {
407  switch(faceDirBasisType)
408  {
409  case LibUtilities::eModified_A:
410  {
411  const LibUtilities::PointsKey pkey(
412  numpoints, LibUtilities::eGaussLobattoLegendre);
413  return LibUtilities::BasisKey(
414  LibUtilities::eModified_A, nummodes, pkey);
415  }
416  case LibUtilities::eModified_B:
417  case LibUtilities::eModified_C:
418  {
419  switch (facedir)
420  {
421  case 0:
422  {
423  const LibUtilities::PointsKey pkey(
424  numpoints+1,
425  LibUtilities::eGaussLobattoLegendre);
426  return LibUtilities::BasisKey(
427  LibUtilities::eModified_A, nummodes, pkey);
428  }
429  case 1:
430  {
431  const LibUtilities::PointsKey pkey(
432  numpoints+1,
433  LibUtilities::eGaussLobattoLegendre);
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 
445  case LibUtilities::eGLL_Lagrange:
446  {
447  switch (facedir)
448  {
449  case 0:
450  {
451  const LibUtilities::PointsKey pkey(
452  numpoints,
453  LibUtilities::eGaussLobattoLegendre);
454  return LibUtilities::BasisKey(
455  LibUtilities::eOrtho_A, nummodes, pkey);
456  }
457  case 1:
458  {
459  const LibUtilities::PointsKey pkey(
460  numpoints,
461  LibUtilities::eGaussRadauMAlpha1Beta0);
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 
473  case LibUtilities::eOrtho_A:
474  case LibUtilities::eOrtho_B:
475  case LibUtilities::eOrtho_C:
476  {
477  switch (facedir)
478  {
479  case 0:
480  {
481  const LibUtilities::PointsKey pkey(
482  numpoints,
483  LibUtilities::eGaussLobattoLegendre);
484  return LibUtilities::BasisKey(
485  LibUtilities::eOrtho_A, nummodes, pkey);
486  }
487  case 1:
488  {
489  const LibUtilities::PointsKey pkey(
490  numpoints,
491  LibUtilities::eGaussRadauMAlpha1Beta0);
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.
510  return LibUtilities::NullBasisKey;
511  }
512  }//end namespace
513 }//end namespace
Generated on Wed Mar 4 2015 16:33:17 for Nektar++ by   doxygen 1.8.1.2