Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PhysGalerkinProject.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PhysGalerkinProject.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 Physical Space Galerkin Projection methods
33 //
34 ///////////////////////////////////////////////////////////////////////////////
38 
43 
44 namespace Nektar
45 {
46  namespace LibUtilities
47  {
48 
49  // Physical Space Interpolation methods
50 
51 
52  // 1D Interpolation
53  void PhysGalerkinProject1D(const BasisKey &fbasis0,
54  const Array<OneD, const NekDouble>& from,
55  const BasisKey &tbasis0,
56  Array<OneD, NekDouble> &to)
57  {
58  PhysGalerkinProject1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
59  }
60 
61  void PhysGalerkinProject1D(const PointsKey &fpoints0,
62  const Array<OneD, const NekDouble>& from,
63  const PointsKey &tpoints0,
64  Array<OneD, NekDouble> &to)
65  {
66  if(fpoints0 == tpoints0) //check to see if the same
67  {
68  Vmath::Vcopy(fpoints0.GetNumPoints(),from,1,to,1);
69  }
70  else // interpolate
71  {
72  DNekMatSharedPtr GP0;
73 
74  GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
75 
76  NekVector<NekDouble> in(fpoints0.GetNumPoints(),from,eWrapper);
77  NekVector<NekDouble> out(tpoints0.GetNumPoints(),to,eWrapper);
78 
79  GP0->Transpose();
80  out = (*GP0)*in;
81  }
82  }
83 
84  void PhysGalerkinProject1D(const BasisKey &fbasis0,
85  const NekDouble *from,
86  const BasisKey &tbasis0,
87  NekDouble *to)
88  {
89  PhysGalerkinProject1D(fbasis0.GetPointsKey(),from,tbasis0.GetPointsKey(),to);
90  }
91 
92  void PhysGalerkinProject1D(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 GP0;
105 
106  GP0 = PointsManager()[tpoints0]
107  ->GetGalerkinProjection(fpoints0);
108 
109  Blas::Dgemv('T', tpoints0.GetNumPoints(), fpoints0.GetNumPoints(),
110  1.0, GP0->GetPtr().get(), tpoints0.GetNumPoints(),
111  from, 1, 0.0, to, 1);
112  }
113  }
114 
115  // 2D Galerkin Projection
116  void PhysGalerkinProject2D(const BasisKey &fbasis0,
117  const BasisKey &fbasis1,
118  const Array<OneD, const NekDouble>& from,
119  const BasisKey &tbasis0,
120  const BasisKey &tbasis1,
121  Array<OneD, NekDouble> &to)
122  {
123  PhysGalerkinProject2D(fbasis0.GetPointsKey(),fbasis1.GetPointsKey(),from.data(),
124  tbasis0.GetPointsKey(),tbasis1.GetPointsKey(),to.data());
125  }
126 
127  void PhysGalerkinProject2D(const PointsKey &fpoints0,
128  const PointsKey &fpoints1,
129  const Array<OneD, const NekDouble>& from,
130  const PointsKey &tpoints0,
131  const PointsKey &tpoints1,
132  Array<OneD, NekDouble> &to)
133  {
134  PhysGalerkinProject2D(fpoints0,fpoints1,from.data(),tpoints0,tpoints1,to.data());
135  }
136 
137  void PhysGalerkinProject2D(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 GP0,GP1;
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  GP1 = PointsManager()[tpoints1]->GetGalerkinProjection(fpoints1);
159  Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0, from, fnp0,
160  GP1->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  GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
170  Blas::Dgemm('N', 'N', tnp0, tnp1, fnp0, 1.0,
171  GP0->GetPtr().get(),
172  tnp0, wsp.get(), fnp0, 0.0, to, tnp0);
173  }
174  }
175 
176 
177  // 3D Galerkin Projection
178  void PhysGalerkinProject3D(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,
185  Array<OneD, NekDouble> &to)
186  {
188  fbasis1.GetPointsKey(),
189  fbasis2.GetPointsKey(),
190  from.data(),
191  tbasis0.GetPointsKey(),
192  tbasis1.GetPointsKey(),
193  tbasis2.GetPointsKey(),
194  to.data());
195  }
196 
197  void PhysGalerkinProject3D(const PointsKey &fpoints0,
198  const PointsKey &fpoints1,
199  const PointsKey &fpoints2,
200  const Array<OneD, const NekDouble>& from,
201  const PointsKey &tpoints0,
202  const PointsKey &tpoints1,
203  const PointsKey &tpoints2,
204  Array<OneD, NekDouble> &to)
205  {
206  PhysGalerkinProject3D(fpoints0,fpoints1,fpoints2,from.data(),
207  tpoints0,tpoints1,tpoints2,to.data());
208  }
209 
210  void PhysGalerkinProject3D(const PointsKey &fpoints0,
211  const PointsKey &fpoints1,
212  const PointsKey &fpoints2,
213  const NekDouble *from,
214  const PointsKey &tpoints0,
215  const PointsKey &tpoints1,
216  const PointsKey &tpoints2,
217  NekDouble *to)
218  {
219  DNekMatSharedPtr GP0,GP1,GP2;
220 
221  int fnp0 = fpoints0.GetNumPoints();
222  int fnp1 = fpoints1.GetNumPoints();
223  int fnp2 = fpoints2.GetNumPoints();
224  int tnp0 = tpoints0.GetNumPoints();
225  int tnp1 = tpoints1.GetNumPoints();
226  int tnp2 = tpoints2.GetNumPoints();
227 
228  Array<OneD, NekDouble> wsp1(fnp0*tnp1*tnp2);
229  Array<OneD, NekDouble> wsp2(fnp0*fnp1*tnp2);
230 
231  GP2 = PointsManager()[tpoints2]->GetGalerkinProjection(fpoints2);
232  Blas::Dgemm('N', 'T', fnp0*fnp1, tnp2, fnp2, 1.0, from, fnp0*fnp1,
233  GP2->GetPtr().get(), tnp2, 0.0, wsp2.get(), fnp0*fnp1);
234 
235  GP1 = PointsManager()[tpoints1]->GetGalerkinProjection(fpoints1);
236  for(int i = 0; i < tnp2; i++)
237  {
238  Blas::Dgemm('N', 'T', fnp0, tnp1, fnp1, 1.0,
239  wsp2.get()+i*fnp0*fnp1,
240  fnp0, GP1->GetPtr().get(),tnp1, 0.0,
241  wsp1.get()+i*fnp0*tnp1,
242  fnp0);
243  }
244 
245  GP0 = PointsManager()[tpoints0]->GetGalerkinProjection(fpoints0);
246  Blas::Dgemm('N', 'N', tnp0, tnp1*tnp2, fnp0, 1.0,
247  GP0->GetPtr().get(), tnp0, wsp1.get(), fnp0, 0.0,
248  to, tnp0);
249  }
250 
251  } // end of namespace
252 } // end of namespace
253