Nektar++
Interp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Interp.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: Definition of Interpolation methods
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
39 
44 
45 namespace Nektar
46 {
47  namespace LibUtilities
48  {
49 
50  // Physical Space Interpolation methods
51 
52 
53  // 1D Interpolation
54  void Interp1D(const BasisKey &fbasis0,
55  const Array<OneD, const NekDouble>& from,
56  const BasisKey &tbasis0,
58  {
59  Interp1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
60  }
61 
62  void Interp1D(const PointsKey &fpoints0,
63  const Array<OneD, const NekDouble>& from,
64  const PointsKey &tpoints0,
66  {
67  if(fpoints0 == tpoints0) //check to see if the same
68  {
69  Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
70  }
71  else // interpolate
72  {
74 
75  I0 = PointsManager()[fpoints0]->GetI(tpoints0);
76 
77  NekVector<NekDouble> in(fpoints0.GetNumPoints(),from,eWrapper);
78  NekVector<NekDouble> out(tpoints0.GetNumPoints(),to,eWrapper);
79 
80  out = (*I0)*in;
81  }
82  }
83 
84  void Interp1D(const BasisKey &fbasis0,
85  const NekDouble *from,
86  const BasisKey &tbasis0,
87  NekDouble *to)
88  {
89  Interp1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
90  }
91 
92  void Interp1D(const PointsKey &fpoints0,
93  const NekDouble *from,
94  const PointsKey &tpoints0,
95  NekDouble *to)
96  {
97  if(fpoints0 == tpoints0) //check to see if the same
98  {
99  Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
100  }
101  else // interpolate
102  {
103 
104  DNekMatSharedPtr I0;
105 
106  I0 = PointsManager()[fpoints0]
107  ->GetI(tpoints0);
108 
109  Blas::Dgemv('N', tpoints0.GetNumPoints(), fpoints0.GetNumPoints(),
110  1.0, I0->GetPtr().get(), tpoints0.GetNumPoints(),
111  from, 1, 0.0, to, 1);
112  }
113  }
114 
115  // 2D Interpolation
116  void Interp2D(const BasisKey &fbasis0,
117  const BasisKey &fbasis1,
118  const Array<OneD, const NekDouble>& from,
119  const BasisKey &tbasis0,
120  const BasisKey &tbasis1,
122  {
123  Interp2D(fbasis0.GetPointsKey(),fbasis1.GetPointsKey(),from.data(),
124  tbasis0.GetPointsKey(),tbasis1.GetPointsKey(),to.data());
125  }
126 
127  void Interp2D(const PointsKey &fpoints0,
128  const PointsKey &fpoints1,
129  const Array<OneD, const NekDouble>& from,
130  const PointsKey &tpoints0,
131  const PointsKey &tpoints1,
133  {
134  Interp2D(fpoints0,fpoints1,from.data(),tpoints0,tpoints1,to.data());
135  }
136 
137  void Interp2D(const PointsKey &fpoints0,
138  const PointsKey &fpoints1,
139  const NekDouble *from,
140  const PointsKey &tpoints0,
141  const PointsKey &tpoints1,
142  NekDouble *to)
143  {
144  DNekMatSharedPtr I0,I1;
145  Array<OneD, NekDouble> wsp(tpoints1.GetNumPoints()*fpoints0.GetNumPoints()); // fnp0*tnp1
146 
147  int fnp0 = fpoints0.GetNumPoints();
148  int fnp1 = fpoints1.GetNumPoints();
149  int tnp0 = tpoints0.GetNumPoints();
150  int tnp1 = tpoints1.GetNumPoints();
151 
152  if(fpoints1 == tpoints1)
153  {
154  Vmath::Vcopy(fnp0*tnp1,from,1,wsp.get(),1);
155  }
156  else
157  {
158  I1 = PointsManager()[fpoints1]->GetI(tpoints1);
159  Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0, from, fnp0,
160  I1->GetPtr().get(), tnp1, 0.0, wsp.get(), fnp0);
161  }
162 
163  if(fpoints0 == tpoints0)
164  {
165  Vmath::Vcopy(tnp0*tnp1,wsp.get(),1,to,1);
166  }
167  else
168  {
169  I0 = PointsManager()[fpoints0]->GetI(tpoints0);
170  Blas::Dgemm('N', 'N', tnp0, tnp1, fnp0, 1.0, I0->GetPtr().get(),
171  tnp0, wsp.get(), fnp0, 0.0, to, tnp0);
172  }
173  }
174 
175 
176 
177  // 3D interpolation
178  void Interp3D(const BasisKey &fbasis0,
179  const BasisKey &fbasis1,
180  const BasisKey &fbasis2,
181  const Array<OneD, const NekDouble>& from,
182  const BasisKey &tbasis0,
183  const BasisKey &tbasis1,
184  const BasisKey &tbasis2,
186  {
187  Interp3D(fbasis0.GetPointsKey(),fbasis1.GetPointsKey(),
188  fbasis2.GetPointsKey(),from.data(),
189  tbasis0.GetPointsKey(),tbasis1.GetPointsKey(),
190  tbasis2.GetPointsKey(),to.data());
191  }
192 
193 
194  void Interp3D(const PointsKey &fpoints0,
195  const PointsKey &fpoints1,
196  const PointsKey &fpoints2,
197  const Array<OneD, const NekDouble>& from,
198  const PointsKey &tpoints0,
199  const PointsKey &tpoints1,
200  const PointsKey &tpoints2,
202  {
203  Interp3D(fpoints0,fpoints1,fpoints2,from.data(),
204  tpoints0,tpoints1,tpoints2,to.data());
205  }
206 
207  void Interp3D(const PointsKey &fpoints0,
208  const PointsKey &fpoints1,
209  const PointsKey &fpoints2,
210  const NekDouble *from,
211  const PointsKey &tpoints0,
212  const PointsKey &tpoints1,
213  const PointsKey &tpoints2,
214  NekDouble *to)
215  {
216  int i;
217  DNekMatSharedPtr I0, I1, I2;
218 
219  int fnp0 = fpoints0.GetNumPoints();
220  int fnp1 = fpoints1.GetNumPoints();
221  int fnp2 = fpoints2.GetNumPoints();
222  int tnp0 = tpoints0.GetNumPoints();
223  int tnp1 = tpoints1.GetNumPoints();
224  int tnp2 = tpoints2.GetNumPoints();
225 
226  Array<OneD, NekDouble> wsp1(tnp0*tnp1*fnp2);
227  Array<OneD, NekDouble> wsp2(tnp0*fnp1*fnp2);
228 
229  I0 = PointsManager()[fpoints0]->GetI(tpoints0);
230  I1 = PointsManager()[fpoints1]->GetI(tpoints1);
231  I2 = PointsManager()[fpoints2]->GetI(tpoints2);
232 
233  Blas::Dgemm('N', 'N', tnp0, fnp1*fnp2, fnp0, 1.0, I0->GetPtr().get(),
234  tnp0, from, fnp0, 0.0, wsp2.get(), tnp0);
235 
236  for(i = 0; i < fnp2; i++)
237  {
238  Blas::Dgemm('N', 'T', tnp0, tnp1, fnp1, 1.0, wsp2.get()+i*tnp0*fnp1,
239  tnp0, I1->GetPtr().get(), tnp1, 0.0, wsp1.get()+i*tnp0*tnp1, tnp0);
240  }
241 
242  Blas::Dgemm('N', 'T', tnp0*tnp1, tnp2, fnp2, 1.0, wsp1.get(),
243  tnp0*tnp1, I2->GetPtr().get(), tnp2, 0.0, to, tnp0*tnp1);
244  }
245 
246 
247 
248  } // end of namespace
249 } // end of namespace
250 
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis...
Definition: Interp.cpp:116
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis...
Definition: Interp.cpp:178
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:54
unsigned int GetNumPoints() const
Definition: Points.h:106
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