Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdExpansion2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdExpansion2D.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 2D expansion. Typically this inolves physiocal
34 // space operations.
35 //
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 #include <StdRegions/StdExpansion2D.h>
39 
40 #ifdef max
41 #undef max
42 #endif
43 
44 namespace Nektar
45 {
46  namespace StdRegions
47  {
48 
49  StdExpansion2D::StdExpansion2D()
50  {
51  }
52 
53  StdExpansion2D::StdExpansion2D(int numcoeffs,
54  const LibUtilities::BasisKey &Ba,
55  const LibUtilities::BasisKey &Bb):
56  StdExpansion(numcoeffs,2, Ba, Bb)
57  {
58  }
59 
60  StdExpansion2D::StdExpansion2D(const StdExpansion2D &T):
61  StdExpansion(T)
62  {
63  }
64 
65  StdExpansion2D::~StdExpansion2D()
66  {
67  }
68 
69  //----------------------------
70  // Differentiation Methods
71  //----------------------------
72  void StdExpansion2D::PhysTensorDeriv(const Array<OneD, const NekDouble>& inarray,
73  Array<OneD, NekDouble> &outarray_d0,
74  Array<OneD, NekDouble> &outarray_d1)
75  {
76  int nquad0 = m_base[0]->GetNumPoints();
77  int nquad1 = m_base[1]->GetNumPoints();
78 
79  if (outarray_d0.num_elements() > 0) // calculate du/dx_0
80  {
81  DNekMatSharedPtr D0 = m_base[0]->GetD();
82  if(inarray.data() == outarray_d0.data())
83  {
84  Array<OneD, NekDouble> wsp(nquad0 * nquad1);
85  Vmath::Vcopy(nquad0 * nquad1,inarray.get(),1,wsp.get(),1);
86  Blas::Dgemm('N', 'N', nquad0, nquad1, nquad0, 1.0,
87  &(D0->GetPtr())[0], nquad0, &wsp[0], nquad0, 0.0,
88  &outarray_d0[0], nquad0);
89  }
90  else
91  {
92  Blas::Dgemm('N', 'N', nquad0, nquad1, nquad0, 1.0,
93  &(D0->GetPtr())[0], nquad0, &inarray[0], nquad0, 0.0,
94  &outarray_d0[0], nquad0);
95  }
96  }
97 
98  if (outarray_d1.num_elements() > 0) // calculate du/dx_1
99  {
100  DNekMatSharedPtr D1 = m_base[1]->GetD();
101  if(inarray.data() == outarray_d1.data())
102  {
103  Array<OneD, NekDouble> wsp(nquad0 * nquad1);
104  Vmath::Vcopy(nquad0 * nquad1,inarray.get(),1,wsp.get(),1);
105  Blas:: Dgemm('N', 'T', nquad0, nquad1, nquad1, 1.0, &wsp[0], nquad0,
106  &(D1->GetPtr())[0], nquad1, 0.0, &outarray_d1[0], nquad0);
107  }
108  else
109  {
110  Blas:: Dgemm('N', 'T', nquad0, nquad1, nquad1, 1.0, &inarray[0], nquad0,
111  &(D1->GetPtr())[0], nquad1, 0.0, &outarray_d1[0], nquad0);
112  }
113  }
114  }
115 
116  NekDouble StdExpansion2D::v_PhysEvaluate(const Array<OneD, const NekDouble>& coords, const Array<OneD, const NekDouble> & physvals)
117  {
118  Array<OneD, NekDouble> coll(2);
119  Array<OneD, DNekMatSharedPtr> I(2);
120 
121  ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
122  ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] > 1");
123  ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
124  ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] > 1");
125 
126  LocCoordToLocCollapsed(coords,coll);
127 
128  I[0] = m_base[0]->GetI(coll);
129  I[1] = m_base[1]->GetI(coll+1);
130 
131  return v_PhysEvaluate(I,physvals);
132  }
133 
134  NekDouble StdExpansion2D::v_PhysEvaluate(
135  const Array<OneD, DNekMatSharedPtr > &I,
136  const Array<OneD, const NekDouble> &physvals)
137  {
138  NekDouble val;
139  int i;
140  int nq0 = m_base[0]->GetNumPoints();
141  int nq1 = m_base[1]->GetNumPoints();
142  Array<OneD, NekDouble> wsp1(nq1);
143 
144  // interpolate first coordinate direction
145  for (i = 0; i < nq1;++i)
146  {
147  wsp1[i] = Blas::Ddot(nq0, &(I[0]->GetPtr())[0], 1,
148  &physvals[i * nq0], 1);
149  }
150 
151  // interpolate in second coordinate direction
152  val = Blas::Ddot(nq1, I[1]->GetPtr(), 1, wsp1, 1);
153 
154  return val;
155  }
156 
157  //////////////////////////////
158  // Integration Methods
159  //////////////////////////////
160 
161  NekDouble StdExpansion2D::Integral(const Array<OneD, const NekDouble>& inarray,
162  const Array<OneD, const NekDouble>& w0,
163  const Array<OneD, const NekDouble>& w1)
164  {
165  int i;
166  NekDouble Int = 0.0;
167  int nquad0 = m_base[0]->GetNumPoints();
168  int nquad1 = m_base[1]->GetNumPoints();
169  Array<OneD, NekDouble> tmp(nquad0 * nquad1);
170 
171  // multiply by integration constants
172  for (i = 0; i < nquad1; ++i)
173  {
174  Vmath::Vmul(nquad0, &inarray[0] + i*nquad0, 1, w0.get(),
175  1, &tmp[0] + i*nquad0, 1);
176  }
177 
178  for (i = 0; i < nquad0; ++i)
179  {
180  Vmath::Vmul(nquad1, &tmp[0]+ i, nquad0, w1.get(), 1,
181  &tmp[0] + i, nquad0);
182  }
183  Int = Vmath::Vsum(nquad0 * nquad1, tmp, 1);
184 
185  return Int;
186  }
187 
188  void StdExpansion2D::BwdTrans_SumFacKernel(
189  const Array<OneD, const NekDouble>& base0,
190  const Array<OneD, const NekDouble>& base1,
191  const Array<OneD, const NekDouble>& inarray,
192  Array<OneD, NekDouble> &outarray,
193  Array<OneD, NekDouble> &wsp,
194  bool doCheckCollDir0,
195  bool doCheckCollDir1)
196  {
197  v_BwdTrans_SumFacKernel(base0, base1, inarray, outarray, wsp, doCheckCollDir0, doCheckCollDir1);
198  }
199 
200  void StdExpansion2D::IProductWRTBase_SumFacKernel(
201  const Array<OneD, const NekDouble>& base0,
202  const Array<OneD, const NekDouble>& base1,
203  const Array<OneD, const NekDouble>& inarray,
204  Array<OneD, NekDouble> &outarray,
205  Array<OneD, NekDouble> &wsp,
206  bool doCheckCollDir0,
207  bool doCheckCollDir1)
208  {
209  v_IProductWRTBase_SumFacKernel(base0, base1, inarray, outarray, wsp, doCheckCollDir0, doCheckCollDir1);
210  }
211 
212  void StdExpansion2D::v_LaplacianMatrixOp_MatFree(
213  const Array<OneD, const NekDouble> &inarray,
214  Array<OneD,NekDouble> &outarray,
215  const StdRegions::StdMatrixKey &mkey)
216  {
217  if (mkey.GetNVarCoeff() == 0
218  &&!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
219  {
220  // This implementation is only valid when there are no
221  // coefficients associated to the Laplacian operator
222  int nquad0 = m_base[0]->GetNumPoints();
223  int nquad1 = m_base[1]->GetNumPoints();
224  int nqtot = nquad0*nquad1;
225  int nmodes0 = m_base[0]->GetNumModes();
226  int nmodes1 = m_base[1]->GetNumModes();
227  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
228 
229  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
230  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
231 
232  // Allocate temporary storage
233  Array<OneD,NekDouble> wsp0(4*wspsize); // size wspsize
234  Array<OneD,NekDouble> wsp1(wsp0+wspsize); // size 3*wspsize
235 
236  if(!(m_base[0]->Collocation() && m_base[1]->Collocation()))
237  {
238  // LAPLACIAN MATRIX OPERATION
239  // wsp0 = u = B * u_hat
240  // wsp1 = du_dxi1 = D_xi1 * wsp0 = D_xi1 * u
241  // wsp2 = du_dxi2 = D_xi2 * wsp0 = D_xi2 * u
242  BwdTrans_SumFacKernel(base0,base1,inarray,wsp0,wsp1,true,true);
243  LaplacianMatrixOp_MatFree_Kernel(wsp0, outarray, wsp1);
244  }
245  else
246  {
247  LaplacianMatrixOp_MatFree_Kernel(inarray, outarray, wsp1);
248  }
249  }
250  else
251  {
252  StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(
253  inarray,outarray,mkey);
254  }
255  }
256 
257 
258 
259  void StdExpansion2D::v_HelmholtzMatrixOp_MatFree(
260  const Array<OneD, const NekDouble> &inarray,
261  Array<OneD,NekDouble> &outarray,
262  const StdRegions::StdMatrixKey &mkey)
263  {
264  if (mkey.GetNVarCoeff() == 0
265  &&!mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio))
266  {
267  int nquad0 = m_base[0]->GetNumPoints();
268  int nquad1 = m_base[1]->GetNumPoints();
269  int nqtot = nquad0*nquad1;
270  int nmodes0 = m_base[0]->GetNumModes();
271  int nmodes1 = m_base[1]->GetNumModes();
272  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),
273  nquad0*nmodes1);
274  NekDouble lambda =
275  mkey.GetConstFactor(StdRegions::eFactorLambda);
276 
277  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
278  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
279 
280  // Allocate temporary storage
281  Array<OneD,NekDouble> wsp0(5*wspsize); // size wspsize
282  Array<OneD,NekDouble> wsp1(wsp0 + wspsize); // size wspsize
283  Array<OneD,NekDouble> wsp2(wsp0 + 2*wspsize);// size 3*wspsize
284 
285  if (!(m_base[0]->Collocation() && m_base[1]->Collocation()))
286  {
287  // MASS MATRIX OPERATION
288  // The following is being calculated:
289  // wsp0 = B * u_hat = u
290  // wsp1 = W * wsp0
291  // outarray = B^T * wsp1 = B^T * W * B * u_hat = M * u_hat
292  BwdTrans_SumFacKernel (base0, base1, inarray,
293  wsp0, wsp2,true,true);
294  MultiplyByQuadratureMetric (wsp0, wsp1);
295  IProductWRTBase_SumFacKernel(base0, base1, wsp1, outarray,
296  wsp2, true, true);
297 
298  LaplacianMatrixOp_MatFree_Kernel(wsp0, wsp1, wsp2);
299  }
300  else
301  {
302  MultiplyByQuadratureMetric(inarray,outarray);
303  LaplacianMatrixOp_MatFree_Kernel(inarray, wsp1, wsp2);
304  }
305 
306  // outarray = lambda * outarray + wsp1
307  // = (lambda * M + L ) * u_hat
308  Vmath::Svtvp(m_ncoeffs, lambda, &outarray[0], 1,
309  &wsp1[0], 1, &outarray[0], 1);
310  }
311  else
312  {
313  StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(
314  inarray,outarray,mkey);
315  }
316  }
317 
318  } //end namespace
319 } //end namespace
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2