Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PolyEPoints.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PolyEPoints.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: 1D Evenly-Spaced Points
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
41 
42 namespace Nektar
43 {
44  namespace LibUtilities
45  {
47  {
48  // Allocate the storage for points
50 
51  unsigned int npts = m_pointsKey.GetNumPoints();
52  if(npts==1)
53  {
54  m_points[0][0] = 0.0;
55  }
56  else
57  {
58  NekDouble dx = 2.0/(NekDouble)(npts-1);
59  for(unsigned int i=0;i<npts;++i)
60  {
61  m_points[0][i] = -1.0 + i*dx;
62  }
63  }
64  }
65 
67  {
68  // Allocate the storage for points
70 
71  unsigned int npts = m_pointsKey.GetNumPoints();
72  if(npts==1)
73  {
74  m_weights[0] = 2.0; //midpoint rule
75  }
76  else
77  {
78  PointsKey gaussKey(npts, eGaussLobattoLegendre);
79  boost::shared_ptr<PointsBaseType> ptr = PointsManager()[gaussKey];
80  Array<OneD, const NekDouble> z;
81  Array<OneD, const NekDouble> w;
82 
83  ptr->GetZW(z,w);
84  for(unsigned int i=0; i<npts;++i)
85  {
86  m_weights[i] = 0.0;
87  for(unsigned j=0;j<npts;++j)
88  {
89  m_weights[i] += w[j]*LagrangePoly(z[j],i,npts,m_points[0]);
90  }
91  }
92  }
93  }
94 
96  {
97  // Allocate the derivative matrix.
99 
100  for(unsigned int i=0;i<m_pointsKey.GetNumPoints();++i)
101  {
102  for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
103  {
105  }
106  }
107  }
108 
109  void PolyEPoints::CalculateInterpMatrix(unsigned int npts, const Array<OneD, const NekDouble>& xpoints, Array<OneD, NekDouble>& interp)
110  {
111  for(unsigned int i=0;i<npts;++i)
112  {
113  for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
114  {
115  interp[i + j*npts] = LagrangePoly(xpoints[i],j,m_pointsKey.GetNumPoints(),m_points[0]);
116  }
117  }
118  }
119 
120  boost::shared_ptr<PolyEPoints::PointsBaseType> PolyEPoints::Create(const PointsKey &key)
121  {
122  boost::shared_ptr<PointsBaseType> returnval(MemoryManager<PolyEPoints>::AllocateSharedPtr(key));
123 
124  returnval->Initialize();
125 
126  return returnval;
127  }
128 
129  const boost::shared_ptr<NekMatrix<NekDouble> > PolyEPoints::GetI(const PointsKey &pkey)
130  {
131  ASSERTL0(pkey.GetPointsDim()==1, "Gauss Points can only interp to other 1d point distributions");
132 
133  int numpoints = pkey.GetNumPoints();
134  Array<OneD, const NekDouble> xpoints;
135 
136  PointsManager()[pkey]->GetPoints(xpoints);
137 
138  return GetI(numpoints, xpoints);
139  }
140 
141  const boost::shared_ptr<NekMatrix<NekDouble> > PolyEPoints::GetI(const Array<OneD, const NekDouble>& x)
142  {
143  int numpoints = 1;
144 
145  /// Delegate to function below.
146  return GetI(numpoints, x);
147  }
148 
149  const boost::shared_ptr<NekMatrix<NekDouble> > PolyEPoints::GetI(unsigned int numpoints, const Array<OneD, const NekDouble>& x)
150  {
151  Array<OneD, NekDouble> interp(GetNumPoints()*numpoints);
152 
153  CalculateInterpMatrix(numpoints, x, interp);
154 
155  unsigned int np = GetTotNumPoints();
156  NekDouble* d = interp.data();
157  boost::shared_ptr< NekMatrix<NekDouble> > returnval(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(numpoints,np, d));
158 
159  return returnval;
160  }
161 
162 
163  NekDouble PolyEPoints::LagrangeInterpolant(NekDouble x, int npts, const Array<OneD, const NekDouble>& xpts,
164  const Array<OneD, const NekDouble>& funcvals)
165  {
166  NekDouble sum = 0.0;
167 
168  for(int i=0;i<npts;++i)
169  {
170  sum += funcvals[i]*LagrangePoly(x,i,npts,xpts);
171  }
172  return sum;
173  }
174 
175 
176  NekDouble PolyEPoints::LagrangePoly(NekDouble x, int pt, int npts, const Array<OneD, const NekDouble>& xpts)
177  {
178  NekDouble h=1.0;
179 
180  for(int i=0;i<pt; ++i)
181  {
182  h = h * (x - xpts[i])/(xpts[pt]-xpts[i]);
183  }
184 
185  for(int i=pt+1;i<npts;++i)
186  {
187  h = h * (x - xpts[i])/(xpts[pt]-xpts[i]);
188  }
189 
190  return h;
191  }
192 
193  NekDouble PolyEPoints::LagrangePolyDeriv(NekDouble x, int pt, int npts, const Array<OneD, const NekDouble>& xpts)
194  {
195  NekDouble h;
196  NekDouble y=0.0;
197 
198  for(int j=0;j<npts;++j)
199  {
200  if(j!=pt)
201  {
202  h=1.0;
203  for(int i=0;i<npts;++i)
204  {
205  if(i!=pt)
206  {
207  if(i!=j)
208  {
209  h = h*(x-xpts[i]);
210  }
211  h = h/(xpts[pt]-xpts[i]);
212  }
213  }
214  y = y + h;
215  }
216  }
217  return y;
218  }
219  } // end of namespace LibUtilities
220 } // end of namespace Nektar
221 
222 
223