Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FourierPoints.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FourierPoints.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 
39 
40 
42 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
43 
44 namespace Nektar
45 {
46  namespace LibUtilities
47  {
48 
50  {
51  // Allocate the storage for points
53 
54  unsigned int npts = m_pointsKey.GetNumPoints();
55  ASSERTL0(!(npts%2), "Fourier points need to be of even order");
56 
57  if(npts==1)
58  {
59  m_points[0][0] = 0.0;
60  }
61  else
62  {
63  NekDouble dx = 2.0/(NekDouble)(npts);
64  for(unsigned int i=0;i<npts;++i)
65  {
66  m_points[0][i] = -1.0 + i*dx;
67  }
68  }
69  }
70 
72  {
73  // Allocate the storage for points
75 
76  unsigned int npts = m_pointsKey.GetNumPoints();
77  if(npts==1)
78  {
79  m_weights[0] = 1.0; //midpoint rule
80  }
81  else
82  {
83  for(unsigned int i=0; i<npts; ++i)
84  {
85  m_weights[i] = 1.0/(NekDouble)npts;
86  }
87  }
88  }
89 
91  {
92  //// Allocate the derivative matrix.
94 
95  unsigned int npts = m_pointsKey.GetNumPoints();
96 
97  for(unsigned int i=1;i<npts;++i)
98  {
99  m_derivmatrix[0]->SetValue(i,i, 0.0);
100  }
101 
102  for(unsigned int i=1;i<npts;++i)
103  {
104  m_derivmatrix[0]->SetValue(0,i, -0.5*M_PI*pow(-1.0,NekDouble(i))*cos(M_PI*i/npts)/sin(M_PI*i/npts));
105  }
106 
107  for(unsigned int i=1;i<npts;++i)
108  {
109  for(unsigned int j=0;j<npts;++j)
110  {
111  m_derivmatrix[0]->SetValue(i,j, (*m_derivmatrix[0])(0,(j-i+npts)%npts));
112  }
113  }
114  }
115 
116  boost::shared_ptr<Points<NekDouble> > FourierPoints::Create(const PointsKey &key)
117  {
118  boost::shared_ptr<Points<NekDouble> > returnval(MemoryManager<FourierPoints>::AllocateSharedPtr(key));
119 
120  returnval->Initialize();
121 
122  return returnval;
123  }
124 
125 
126  boost::shared_ptr< NekMatrix<NekDouble> > FourierPoints::CreateMatrix(const PointsKey &pkey)
127  {
128  int numpoints = pkey.GetNumPoints();
130 
131  PointsManager()[pkey]->GetPoints(xpoints);
132 
133  /// Delegate to function below.
134  return GetI(numpoints, xpoints);
135  }
136 
137  const boost::shared_ptr<NekMatrix<NekDouble> > FourierPoints::GetI(const PointsKey& pkey)
138  {
139  ASSERTL0(pkey.GetPointsDim()==1, "Fourier Points can only interp to other 1d point distributions");
140 
141  return m_InterpManager[pkey];
142  }
143 
144  const boost::shared_ptr<NekMatrix<NekDouble> > FourierPoints::GetI(const Array<OneD, const NekDouble>& x)
145  {
146  int numpoints = 1;
147 
148  /// Delegate to function below.
149  return GetI(numpoints, x);
150  }
151 
152  const boost::shared_ptr<NekMatrix<NekDouble> > FourierPoints::GetI(unsigned int numpoints, const Array<OneD, const NekDouble>& x)
153  {
154  Array<OneD, NekDouble> interp(GetNumPoints()*numpoints);
155 
156  CalculateInterpMatrix(numpoints, x, interp);
157 
158  NekDouble* t = interp.data();
159  unsigned int np = GetNumPoints();
160  boost::shared_ptr< NekMatrix<NekDouble> > returnval(MemoryManager<NekMatrix<NekDouble> >::AllocateSharedPtr(numpoints,np,t));
161 
162  return returnval;
163  }
164 
166  {
167  const NekDouble h = 2.0/m_pointsKey.GetNumPoints();
168  for(unsigned int i=0;i<npts;++i)
169  {
170  for(unsigned int j=0;j<m_pointsKey.GetNumPoints();++j)
171  {
172  interp[i*m_pointsKey.GetNumPoints()+j] = PeriodicSincFunction(M_PI*(xpoints[i]-m_points[0][j]),h);
173  }
174  }
175  }
176 
178  {
179  // This formula is based upon a mapped version of
180  // the periodic sinc presented in Trefethen's "Spectral Methods
181  // in Matlab"
182 
183  NekDouble y = 1.0;
184 
185  if(fabs(x)>1.0e-12)
186  {
187  y = sin(M_PI*x/(M_PI*h))/((2.0/h)*tan(0.5*x));
188  }
189 
190  return y;
191  }
192 
193  } // end of namespace LibUtilities
194 } // end of namespace Nektar
195 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
MatrixSharedPtrType m_derivmatrix[3]
Definition: Points.h:351
Array< OneD, DataType > m_points[3]
Definition: Points.h:349
boost::shared_ptr< NekMatrix< NekDouble > > CreateMatrix(const PointsKey &pkey)
void CalculateInterpMatrix(unsigned int npts, const Array< OneD, const NekDouble > &xpoints, Array< OneD, NekDouble > &interp)
unsigned int GetNumPoints() const
Definition: Points.h:246
static std::string npts
Definition: InputFld.cpp:43
unsigned int GetPointsDim() const
Definition: Points.h:149
PointsManagerT & PointsManager(void)
Array< OneD, DataType > m_weights
Definition: Points.h:350
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
static boost::shared_ptr< PointsBaseType > Create(const PointsKey &key)
virtual void CalculateDerivMatrix()
Definition: Points.h:371
const MatrixSharedPtrType GetI(const PointsKey &pkey)
unsigned int GetNumPoints() const
Definition: Points.h:106
NekDouble PeriodicSincFunction(const NekDouble x, const NekDouble h)
NekManager< PointsKey, NekMatrix< DataType >, PointsKey::opLess > m_InterpManager
Definition: Points.h:352