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 {
47  namespace StdRegions
48  {
49  StdExpansion3D::StdExpansion3D()
50  {
51  }
52 
53  StdExpansion3D::StdExpansion3D(int numcoeffs,
54  const LibUtilities::BasisKey &Ba,
55  const LibUtilities::BasisKey &Bb,
56  const LibUtilities::BasisKey &Bc) :
57  StdExpansion(numcoeffs,3,Ba,Bb,Bc)
58  {
59  }
60 
61  StdExpansion3D::StdExpansion3D(const StdExpansion3D &T):
62  StdExpansion(T)
63  {
64  }
65 
66  StdExpansion3D::~StdExpansion3D()
67  {
68  }
69  void StdExpansion3D::PhysTensorDeriv(
70  const Array<OneD, const NekDouble> &inarray,
71  Array<OneD, NekDouble> &out_dx,
72  Array<OneD, NekDouble> &out_dy,
73  Array<OneD, NekDouble> &out_dz)
74  {
75  const int nquad0 = m_base[0]->GetNumPoints();
76  const int nquad1 = m_base[1]->GetNumPoints();
77  const int nquad2 = m_base[2]->GetNumPoints();
78 
79  Array<OneD, NekDouble> wsp(nquad0*nquad1*nquad2);
80 
81  // copy inarray to wsp in case inarray is used as outarray
82  Vmath::Vcopy(nquad0*nquad1*nquad2, &inarray[0], 1, &wsp[0], 1);
83 
84  if (out_dx.num_elements() > 0)
85  {
86  NekDouble *D0 = &((m_base[0]->GetD())->GetPtr())[0];
87 
88  Blas::Dgemm('N','N', nquad0,nquad1*nquad2,nquad0,1.0,
89  D0,nquad0,&wsp[0],nquad0,0.0,&out_dx[0],nquad0);
90  }
91 
92  if (out_dy.num_elements() > 0)
93  {
94  NekDouble *D1 = &((m_base[1]->GetD())->GetPtr())[0];
95  for (int j = 0; j < nquad2; ++j)
96  {
97  Blas::Dgemm('N', 'T', nquad0, nquad1, nquad1,
98  1.0, &wsp[j*nquad0*nquad1], nquad0,
99  D1, nquad1,
100  0.0, &out_dy[j*nquad0*nquad1], nquad0);
101  }
102  }
103 
104  if (out_dz.num_elements() > 0)
105  {
106  NekDouble *D2 = &((m_base[2]->GetD())->GetPtr())[0];
107 
108  Blas::Dgemm('N','T',nquad0*nquad1,nquad2,nquad2,1.0,
109  &wsp[0],nquad0*nquad1,D2,nquad2,0.0,&out_dz[0],
110  nquad0*nquad1);
111  }
112  }
113 
114  void StdExpansion3D::BwdTrans_SumFacKernel(
115  const Array<OneD, const NekDouble>& base0,
116  const Array<OneD, const NekDouble>& base1,
117  const Array<OneD, const NekDouble>& base2,
118  const Array<OneD, const NekDouble>& inarray,
119  Array<OneD, NekDouble>& outarray,
120  Array<OneD, NekDouble>& wsp,
121  bool doCheckCollDir0,
122  bool doCheckCollDir1,
123  bool doCheckCollDir2)
124  {
125  v_BwdTrans_SumFacKernel(base0, base1, base2, inarray, outarray, wsp, doCheckCollDir0, doCheckCollDir1, doCheckCollDir2);
126  }
127 
128  void StdExpansion3D::IProductWRTBase_SumFacKernel(
129  const Array<OneD, const NekDouble>& base0,
130  const Array<OneD, const NekDouble>& base1,
131  const Array<OneD, const NekDouble>& base2,
132  const Array<OneD, const NekDouble>& inarray,
133  Array<OneD, NekDouble> &outarray,
134  Array<OneD, NekDouble> &wsp,
135  bool doCheckCollDir0,
136  bool doCheckCollDir1,
137  bool doCheckCollDir2)
138  {
139  v_IProductWRTBase_SumFacKernel(base0, base1, base2, inarray, outarray, wsp, doCheckCollDir0, doCheckCollDir1, doCheckCollDir2);
140  }
141 
142  NekDouble StdExpansion3D::v_PhysEvaluate(
143  const Array<OneD, const NekDouble> &coords,
144  const Array<OneD, const NekDouble> &physvals)
145  {
146  Array<OneD, NekDouble> eta = Array<OneD, NekDouble>(3);
147  Array<OneD, DNekMatSharedPtr> I(3);
148 
149  WARNINGL2(coords[0] >= -1,"coord[0] < -1");
150  WARNINGL2(coords[0] <= 1,"coord[0] > 1");
151  WARNINGL2(coords[1] >= -1,"coord[1] < -1");
152  WARNINGL2(coords[1] <= 1,"coord[1] > 1");
153  WARNINGL2(coords[2] >= -1,"coord[2] < -1");
154  WARNINGL2(coords[2] <= 1,"coord[2] > 1");
155 
156  // Obtain local collapsed corodinate from
157  // cartesian coordinate.
158  LocCoordToLocCollapsed(coords,eta);
159 
160  // Get Lagrange interpolants.
161  I[0] = m_base[0]->GetI(eta);
162  I[1] = m_base[1]->GetI(eta+1);
163  I[2] = m_base[2]->GetI(eta+2);
164 
165  return v_PhysEvaluate(I,physvals);
166  }
167 
168  NekDouble StdExpansion3D::v_PhysEvaluate(
169  const Array<OneD, DNekMatSharedPtr > &I,
170  const Array<OneD, const NekDouble> &physvals)
171  {
172  NekDouble value;
173 
174  int Qx = m_base[0]->GetNumPoints();
175  int Qy = m_base[1]->GetNumPoints();
176  int Qz = m_base[2]->GetNumPoints();
177 
178  Array<OneD, NekDouble> sumFactorization_qr = Array<OneD, NekDouble>(Qy*Qz);
179  Array<OneD, NekDouble> sumFactorization_r = Array<OneD, NekDouble>(Qz);
180 
181  // Lagrangian interpolation matrix
182  NekDouble *interpolatingNodes = 0;
183 
184  // Interpolate first coordinate direction
185  interpolatingNodes = &I[0]->GetPtr()[0];
186 
187  Blas::Dgemv('T',Qx,Qy*Qz,1.0,&physvals[0],Qx,&interpolatingNodes[0], 1, 0.0, &sumFactorization_qr[0], 1);
188 
189  // Interpolate in second coordinate direction
190  interpolatingNodes = &I[1]->GetPtr()[0];
191 
192  Blas::Dgemv('T',Qy,Qz,1.0,&sumFactorization_qr[0],Qy,&interpolatingNodes[0],1,0.0,&sumFactorization_r[0], 1);
193 
194  // Interpolate in third coordinate direction
195  interpolatingNodes = &I[2]->GetPtr()[0];
196  value = Blas::Ddot(Qz, interpolatingNodes, 1, &sumFactorization_r[0], 1);
197 
198  return value;
199  }
200 
201 
202  /**
203  * @param inarray Input coefficients.
204  * @param output Output coefficients.
205  * @param mkey Matrix key
206  */
207  void StdExpansion3D::v_LaplacianMatrixOp_MatFree(
208  const Array<OneD, const NekDouble> &inarray,
209  Array<OneD,NekDouble> &outarray,
210  const StdRegions::StdMatrixKey &mkey)
211  {
212  if ( mkey.GetNVarCoeff() == 0 &&
213  !mkey.ConstFactorExists(eFactorSVVCutoffRatio))
214  {
215  // This implementation is only valid when there are no
216  // coefficients associated to the Laplacian operator
217  int nqtot = GetTotPoints();
218 
219  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
220  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
221  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata();
222 
223  // Allocate temporary storage
224  Array<OneD,NekDouble> wsp0(7*nqtot);
225  Array<OneD,NekDouble> wsp1(wsp0+nqtot);
226 
227  if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
228  m_base[2]->Collocation()))
229  {
230  // LAPLACIAN MATRIX OPERATION
231  // wsp0 = u = B * u_hat
232  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
233  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
234  BwdTrans_SumFacKernel(base0,base1,base2,inarray,wsp0,wsp1,true,true,true);
235  LaplacianMatrixOp_MatFree_Kernel(wsp0,outarray,wsp1);
236  }
237  else
238  {
239  LaplacianMatrixOp_MatFree_Kernel(inarray,outarray,wsp1);
240  }
241  }
242  else
243  {
244  StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
245  }
246  }
247 
248 
249  void StdExpansion3D::v_HelmholtzMatrixOp_MatFree(
250  const Array<OneD, const NekDouble> &inarray,
251  Array<OneD,NekDouble> &outarray,
252  const StdRegions::StdMatrixKey &mkey)
253  {
254  if(mkey.GetNVarCoeff() == 0)
255  {
256  int nquad0 = m_base[0]->GetNumPoints();
257  int nquad1 = m_base[1]->GetNumPoints();
258  int nquad2 = m_base[2]->GetNumPoints();
259  int nmodes0 = m_base[0]->GetNumModes();
260  int nmodes1 = m_base[1]->GetNumModes();
261  int nmodes2 = m_base[2]->GetNumModes();
262  int wspsize = max(nquad0*nmodes2*(nmodes1+nquad1),
263  nquad0*nquad1*(nquad2+nmodes0)+
264  nmodes0*nmodes1*nquad2);
265 
266  NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
267 
268  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata ();
269  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata ();
270  const Array<OneD, const NekDouble>& base2 = m_base[2]->GetBdata ();
271  Array<OneD,NekDouble> wsp0(8*wspsize);
272  Array<OneD,NekDouble> wsp1(wsp0+1*wspsize);
273  Array<OneD,NekDouble> wsp2(wsp0+2*wspsize);
274 
275  if(!(m_base[0]->Collocation() && m_base[1]->Collocation() &&
276  m_base[2]->Collocation()))
277  {
278  // MASS MATRIX OPERATION
279  // The following is being calculated:
280  // wsp0 = B * u_hat = u
281  // wsp1 = W * wsp0
282  // outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
283  BwdTrans_SumFacKernel (base0,base1,base2,inarray,
284  wsp0,wsp2,true,true,true);
285  MultiplyByQuadratureMetric (wsp0,wsp1);
286  IProductWRTBase_SumFacKernel (base0,base1,base2,wsp1,
287  outarray,wsp2,true,true,true);
288  LaplacianMatrixOp_MatFree_Kernel(wsp0,wsp1,wsp2);
289  }
290  else
291  {
292  // specialised implementation for the classical spectral
293  // element method
294  MultiplyByQuadratureMetric (inarray,outarray);
295  LaplacianMatrixOp_MatFree_Kernel(inarray,wsp1,wsp2);
296  }
297 
298  // outarray = lambda * outarray + wsp1
299  // = (lambda * M + L ) * u_hat
300  Vmath::Svtvp(m_ncoeffs,lambda,&outarray[0],1,&wsp1[0],1,
301  &outarray[0],1);
302  }
303  else
304  {
305  StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(inarray,outarray,mkey);
306  }
307 
308  }
309 
310  const NormalVector & StdExpansion3D::v_GetSurfaceNormal(
311  const int id) const
312  {
313  return v_GetFaceNormal(id);
314  }
315 
316  const NormalVector & StdExpansion3D::v_GetFaceNormal(const int face) const
317  {
318  std::map<int, NormalVector>::const_iterator x;
319  x = m_faceNormals.find(face);
320  ASSERTL0 (x != m_faceNormals.end(),
321  "face normal not computed.");
322  return x->second;
323  }
324 
325  NekDouble StdExpansion3D::v_Integral(
326  const Array<OneD, const NekDouble>& inarray)
327  {
328  const int nqtot = GetTotPoints();
329  Array<OneD, NekDouble> tmp(GetTotPoints());
330  MultiplyByStdQuadratureMetric(inarray, tmp);
331  return Vmath::Vsum(nqtot, tmp, 1);
332  }
333 
334 
335  void StdExpansion3D::v_NegateFaceNormal(const int face)
336  {
337  m_negatedNormals[face] = true;
338  for (int i = 0; i < GetCoordim(); ++i)
339  {
340  Vmath::Neg(m_faceNormals[face][i].num_elements(),
341  m_faceNormals[face][i], 1);
342  }
343  }
344  }//end namespace
345 }//end namespace
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2