Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GaussPoints.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File GaussPoints.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: GaussPoints Definitions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/Foundations/Points.h>
37 #include <LibUtilities/Foundations/Foundations.hpp>
38 
39 #include <LibUtilities/BasicUtils/ErrorUtil.hpp>
40 #include <LibUtilities/Polylib/Polylib.h>
41 #include <LibUtilities/Foundations/GaussPoints.h>
42 #include <LibUtilities/Foundations/ManagerAccess.h>
43 #include <LibUtilities/Memory/NekMemoryManager.hpp>
44 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
45 
46 
47 namespace Nektar
48 {
49  namespace LibUtilities
50  {
51  void GaussPoints::CalculatePoints()
52  {
53  // Allocate the storage for points and weights
54  PointsBaseType::CalculatePoints();
55  PointsBaseType::CalculateWeights();
56 
57  int numpoints = m_pointsKey.GetNumPoints();
58 
59  switch(m_pointsKey.GetPointsType())
60  {
61  case eGaussGaussLegendre:
62  Polylib::zwgj(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
63  break;
64 
65  case eGaussRadauMLegendre:
66  Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
67  break;
68 
69  case eGaussRadauPLegendre:
70  Polylib::zwgrjp(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
71  break;
72 
73  case eGaussLobattoLegendre:
74  Polylib::zwglj(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
75  break;
76 
77  case eGaussGaussChebyshev:
78  Polylib::zwgj(m_points[0].data(),m_weights.data(),numpoints,-0.5,-0.5);
79  break;
80 
81  case eGaussRadauMChebyshev:
82  Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,-0.5,-0.5);
83  break;
84 
85  case eGaussRadauPChebyshev:
86  Polylib::zwgrjp(m_points[0].data(),m_weights.data(),numpoints,-0.5,-0.5);
87  break;
88 
89  case eGaussLobattoChebyshev:
90  Polylib::zwglj(m_points[0].data(),m_weights.data(),numpoints,-0.5,-0.5);
91  break;
92 
93  case eGaussRadauMAlpha0Beta1:
94  Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,0.0,1.0);
95  break;
96 
97  case eGaussRadauMAlpha0Beta2:
98  Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,0.0,2.0);
99  break;
100 
101  case eGaussRadauMAlpha1Beta0:
102  Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,1.0,0.0);
103  break;
104 
105  case eGaussRadauMAlpha2Beta0:
106  Polylib::zwgrjm(m_points[0].data(),m_weights.data(),numpoints,2.0,0.0);
107  break;
108 
109  case eGaussKronrodLegendre:
110  Polylib::zwgk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
111  break;
112 
113  case eGaussRadauKronrodMLegendre:
114  Polylib::zwrk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
115  break;
116 
117  case eGaussRadauKronrodMAlpha1Beta0:
118  Polylib::zwrk(m_points[0].data(),m_weights.data(),numpoints,1.0,0.0);
119  break;
120 
121  case eGaussLobattoKronrodLegendre:
122  Polylib::zwlk(m_points[0].data(),m_weights.data(),numpoints,0.0,0.0);
123  break;
124 
125  default:
126  ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
127  }
128  }
129 
130  void GaussPoints::CalculateWeights()
131  {
132  // For Gauss Quadrature, this is done as part of the points computation
133  }
134 
135  void GaussPoints::CalculateDerivMatrix()
136  {
137  // Allocate the derivative matrix
138  PointsBaseType::CalculateDerivMatrix();
139 
140  int numpoints = m_pointsKey.GetNumPoints();
141  int totpoints = m_pointsKey.GetTotNumPoints();
142  Array<OneD, NekDouble> dmtempSharedArray = Array<OneD, NekDouble>(totpoints*totpoints);
143  NekDouble *dmtemp = dmtempSharedArray.data();
144 
145  switch(m_pointsKey.GetPointsType())
146  {
147  case eGaussGaussLegendre:
148  Polylib::Dgj(dmtemp,m_points[0].data(),numpoints,0.0,0.0);
149  break;
150 
151  case eGaussRadauMLegendre:
152  Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,0.0,0.0);
153  break;
154 
155  case eGaussRadauPLegendre:
156  Polylib::Dgrjp(dmtemp,m_points[0].data(),numpoints,0.0,0.0);
157  break;
158 
159  case eGaussLobattoLegendre:
160  Polylib::Dglj(dmtemp,m_points[0].data(),numpoints,0.0,0.0);
161  break;
162 
163  case eGaussGaussChebyshev:
164  Polylib::Dgj(dmtemp,m_points[0].data(),numpoints,-0.5,-0.5);
165  break;
166 
167  case eGaussRadauMChebyshev:
168  Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,-0.5,-0.5);
169  break;
170 
171  case eGaussRadauPChebyshev:
172  Polylib::Dgrjp(dmtemp,m_points[0].data(),numpoints,-0.5,-0.5);
173  break;
174 
175  case eGaussLobattoChebyshev:
176  Polylib::Dglj(dmtemp,m_points[0].data(),numpoints,-0.5,-0.5);
177  break;
178 
179  case eGaussRadauMAlpha0Beta1:
180  Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,0.0,1.0);
181  break;
182 
183  case eGaussRadauMAlpha0Beta2:
184  Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,0.0,2.0);
185  break;
186 
187  case eGaussRadauMAlpha1Beta0:
188  Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,1.0,0.0);
189  break;
190 
191  case eGaussRadauMAlpha2Beta0:
192  Polylib::Dgrjm(dmtemp,m_points[0].data(),numpoints,2.0,0.0);
193  break;
194 
195  case eGaussKronrodLegendre:
196  case eGaussRadauKronrodMLegendre:
197  case eGaussRadauKronrodMAlpha1Beta0:
198  case eGaussLobattoKronrodLegendre:
199  {
200  for(unsigned int i=0;i<m_pointsKey.GetNumPoints();++i)
201  {
202  for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
203  {
204  (*m_derivmatrix[0])(i,j) = LagrangePolyDeriv(m_points[0][i],j,m_pointsKey.GetNumPoints(),m_points[0]);
205  }
206  }
207  return;
208  }
209  break;
210 
211  default:
212  ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
213  }
214 
215  std::copy(dmtemp,dmtemp+totpoints*totpoints,m_derivmatrix[0]->begin());
216  }
217 
218  void GaussPoints::CalculateInterpMatrix(unsigned int npts, const Array<OneD, const NekDouble>& xpoints, Array<OneD, NekDouble>& interp)
219  {
220  switch(m_pointsKey.GetPointsType())
221  {
222  case eGaussGaussLegendre:
223  Polylib::Imgj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
224  break;
225 
226  case eGaussRadauMLegendre:
227  Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
228  break;
229 
230  case eGaussRadauPLegendre:
231  Polylib::Imgrjp(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
232  break;
233 
234  case eGaussLobattoLegendre:
235  Polylib::Imglj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,0.0);
236  break;
237 
238  case eGaussGaussChebyshev:
239  Polylib::Imgj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
240  break;
241 
242  case eGaussRadauMChebyshev:
243  Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
244  break;
245 
246  case eGaussRadauPChebyshev:
247  Polylib::Imgrjp(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
248  break;
249 
250  case eGaussLobattoChebyshev:
251  Polylib::Imglj(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,-0.5,-0.5);
252  break;
253 
254  case eGaussRadauMAlpha0Beta1:
255  Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,1.0);
256  break;
257 
258  case eGaussRadauMAlpha0Beta2:
259  Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,0.0,2.0);
260  break;
261 
262  case eGaussRadauMAlpha1Beta0:
263  Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,1.0,0.0);
264  break;
265 
266  case eGaussRadauMAlpha2Beta0:
267  Polylib::Imgrjm(interp.data(),m_points[0].data(),xpoints.data(),GetNumPoints(),npts,2.0,0.0);
268  break;
269 
270  case eGaussKronrodLegendre:
271  case eGaussRadauKronrodMLegendre:
272  case eGaussRadauKronrodMAlpha1Beta0:
273  case eGaussLobattoKronrodLegendre:
274  {
275  for(unsigned int i=0;i<npts;++i)
276  {
277  for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
278  {
279  interp[i + j*npts] = LagrangePoly(xpoints[i],j,m_pointsKey.GetNumPoints(),m_points[0]);
280  }
281  }
282  }
283  break;
284 
285  default:
286  ASSERTL0(false, "Unknown Gauss quadrature point distribution requested");
287  }
288  }
289 
290  boost::shared_ptr<Points<NekDouble> > GaussPoints::Create(const PointsKey &pkey)
291  {
292  boost::shared_ptr< Points<NekDouble> > returnval(MemoryManager< GaussPoints >::AllocateSharedPtr(pkey));
293 
294  returnval->Initialize();
295 
296  return returnval;
297  }
298 
299  boost::shared_ptr< NekMatrix<NekDouble> > GaussPoints::CreateMatrix(const PointsKey &pkey)
300  {
301  int numpoints = pkey.GetNumPoints();
302  Array<OneD, const NekDouble> xpoints;
303 
304  PointsManager()[pkey]->GetPoints(xpoints);
305 
306  // Delegate to function below
307  return GetI(numpoints, xpoints);
308  }
309 
310 
311  const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetI(const PointsKey &pkey)
312  {
313  ASSERTL0(pkey.GetPointsDim()==1, "Gauss Points can only interp to other 1d point distributions");
314 
315  return m_InterpManager[pkey];
316  }
317 
318  const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetI(const Array<OneD, const NekDouble>& x)
319  {
320  int numpoints = 1;
321 
322  // Delegate to function below
323  return GetI(numpoints, x);
324  }
325 
326  const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetI(unsigned int numpoints, const Array<OneD, const NekDouble>& x)
327  {
328  Array<OneD, NekDouble> interp(GetNumPoints()*numpoints);
329 
330  CalculateInterpMatrix(numpoints, x, interp);
331 
332  NekDouble* t = interp.data();
333  unsigned int np = GetNumPoints();
334  boost::shared_ptr< NekMatrix<NekDouble> > returnval(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(numpoints,np,t));
335 
336  return returnval;
337  }
338 
339  NekDouble GaussPoints::LagrangeInterpolant(NekDouble x, int npts, const Array<OneD, const NekDouble>& xpts,
340  const Array<OneD, const NekDouble>& funcvals)
341  {
342  NekDouble sum = 0.0;
343 
344  for(int i=0;i<npts;++i)
345  {
346  sum += funcvals[i]*LagrangePoly(x,i,npts,xpts);
347  }
348  return sum;
349  }
350 
351 
352  NekDouble GaussPoints::LagrangePoly(NekDouble x, int pt, int npts, const Array<OneD, const NekDouble>& xpts)
353  {
354  NekDouble h=1.0;
355 
356  for(int i=0;i<pt; ++i)
357  {
358  h = h * (x - xpts[i])/(xpts[pt]-xpts[i]);
359  }
360 
361  for(int i=pt+1;i<npts;++i)
362  {
363  h = h * (x - xpts[i])/(xpts[pt]-xpts[i]);
364  }
365 
366  return h;
367  }
368 
369  NekDouble GaussPoints::LagrangePolyDeriv(NekDouble x, int pt, int npts, const Array<OneD, const NekDouble>& xpts)
370  {
371  NekDouble h;
372  NekDouble y=0.0;
373 
374  for(int j=0;j<npts;++j)
375  {
376  if(j!=pt)
377  {
378  h=1.0;
379  for(int i=0;i<npts;++i)
380  {
381  if(i!=pt)
382  {
383  if(i!=j)
384  {
385  h = h*(x-xpts[i]);
386  }
387  h = h/(xpts[pt]-xpts[i]);
388  }
389  }
390  y = y + h;
391  }
392  }
393  return y;
394  }
395 
396  const boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::GetGalerkinProjection(const PointsKey &pkey)
397  {
398  return m_GalerkinProjectionManager[pkey];
399  }
400 
401  boost::shared_ptr< NekMatrix<NekDouble> > GaussPoints::CreateGPMatrix(const PointsKey &pkey)
402  {
403  boost::shared_ptr< NekMatrix<NekDouble> > returnval = CalculateGalerkinProjectionMatrix(pkey);
404 
405  // Delegate to function below
406  return returnval;
407  }
408 
409  boost::shared_ptr<NekMatrix<NekDouble> > GaussPoints::CalculateGalerkinProjectionMatrix(const PointsKey &pkey)
410  {
411  int numpointsfrom = pkey.GetNumPoints();
412  int numpointsto = GetNumPoints();
413 
414  Array<OneD, const NekDouble> weightsfrom;
415 
416  weightsfrom = PointsManager()[pkey]->GetW();
417 
418  boost::shared_ptr< NekMatrix<NekDouble> > Interp = GetI(pkey);
419 
420  Array<OneD, NekDouble> GalProj(numpointsfrom*numpointsto);
421 
422  // set up inner product matrix and multiply by inverse of
423  // diagaonal mass matrix
424  for(int i = 0; i < numpointsto; ++i)
425  {
426  Vmath::Vmul(numpointsfrom,Interp->GetPtr().get() +i*numpointsfrom,1,
427  &weightsfrom[0],1,&GalProj[0] +i,numpointsto);
428  Vmath::Smul(numpointsfrom,1.0/m_weights[i],&GalProj[0]+i,numpointsto,
429  &GalProj[0]+i,numpointsto);
430  }
431 
432 
433  NekDouble* t = GalProj.data();
434  boost::shared_ptr< NekMatrix<NekDouble> > returnval(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(numpointsto,numpointsfrom,t));
435 
436  return returnval;
437  }
438 
439  } // end of namespace LibUtilities
440 } // end of namespace Nektar
441 
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2