Nektar++
NodalTriEvenlySpaced.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalTriEvenlySpaced.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: 2D Nodal Triangle Evenly Spaced Point Definitions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
42 #include <vector>
43 
44 
45 
46 namespace Nektar
47 {
48  namespace LibUtilities
49  {
50  namespace
51  {
52  // construct the geometory and set the coordinate of triangle
53  // edges and vertices are ordered as anticlockwise
54  bool isVertex(int i, int j, int npts){
55  return (i==0 && j==0) || (i==(npts-1) && j==0) || (i==0 && j==(npts-1));
56  }
57 
58  bool isEdge(int i, int j, int npts){
59  return i==0 || j==0 || i+j==npts-1; //i+j=tot num of steps
60  }
61 
62  bool isEdge_1(int i, int j, int npts){
63  return i==0;
64  }
65 
66  bool isEdge_2(int i, int j, int npts){
67  return i+j==npts-1;
68  }
69  }
70 
71 
73  {
74  // Allocate the storage for points
76 
77  // Populate m_points
78  unsigned int npts = GetNumPoints();
79  NekDouble delta = 2.0/(npts - 1.0);
80  for(int i=0, index=0; i<npts; ++i){ // y-direction
81  for(int j=0; j<npts-i; ++j,++index){ // x-direction
82  NekDouble x = -1.0 + j*delta;
83  NekDouble y = -1.0 + i*delta;
84  m_points[0][index] = x;
85  m_points[1][index] = y;
86  }
87  }
88 
90 
91 
92  }
93 
94 
96  {
97  // Allocate the storage for points
99 
100  typedef DataType T;
101 
102  // Solve the Vandermonde system of integrals for the weight vector
104 
105  m_weights = Array<OneD,T>( w.GetRows(), w.GetPtr() );
106 
107  }
108 
109 
110  // ////////////////////////////////////////
111  // CalculateInterpMatrix()
113  const Array<OneD, const NekDouble>& yia,
114  Array<OneD, NekDouble>& interp)
115  {
118  NekVector<NekDouble> xi( xia );
119  NekVector<NekDouble> yi( yia );
120  NekMatrix<NekDouble> interMat = GetInterpolationMatrix(x, y, xi, yi);
121 
122  int rows = xi.GetRows(), cols = GetTotNumPoints();
123  for( int i = 0; i < rows; ++i ) {
124  for( int j = 0; j < cols; ++j ) {
125  interp[j + i*cols] = interMat(i,j);
126  }
127  }
128  }
129 
130  // ////////////////////////////////////////
131  // CalculateDerivMatrix()
133  {
134 
135  // Allocate the derivative matrix.
137 
140  NekVector<NekDouble> xi = x;
141  NekVector<NekDouble> yi = y;
142 
143  m_derivmatrix[0] = GetXDerivativeMatrix(x,y,xi,yi);
144  m_derivmatrix[1] = GetYDerivativeMatrix(x,y,xi,yi);
145 
146  }
147 
148  boost::shared_ptr<PointsBaseType> NodalTriEvenlySpaced::Create(const PointsKey &key)
149  {
150  boost::shared_ptr<PointsBaseType> returnval(MemoryManager<NodalTriEvenlySpaced>::AllocateSharedPtr(key));
151 
152  returnval->Initialize();
153 
154  return returnval;
155  }
156 
158  {
159  unsigned int npts = GetNumPoints();
160  using std::vector;
161  vector<int> vertex;
162  vector<int> iEdge_1; // interior edge points on the bottom triangle edge
163  vector<int> iEdge_2; // interior edge points on the right triangle edge
164  vector<int> iEdge_3; // interior edge points on the left triangle edge
165  vector<int> interiorPoints;
166  vector<int> map;
167 
168  // Build the lattice triangle left to right - bottom to top
169  for(int i=0, index=0; i<npts; ++i){ // y-direction
170  for(int j=0; j<npts-i; ++j,++index){ // x-direction
171 
172  if( isVertex(i,j,npts) ) {
173 
174  vertex.push_back(index);
175 
176  } else if( isEdge(i,j,npts) ) { // interior edge
177 
178  if(isEdge_1(i,j,npts)){ // bottom edge
179 
180  iEdge_1.push_back(index);
181 
182  }else if(isEdge_2(i,j,npts)){ // right edge
183 
184  iEdge_2.push_back(index);
185 
186  }else // left edge
187  {
188  // Add backwards. This is because of counter clockwise.
189  iEdge_3.insert(iEdge_3.begin(), index);
190  }
191 
192  } else { // Interior points
193 
194  interiorPoints.push_back(index);
195  }
196  }
197  }
198 
199  // Mapping the vertex, edges, and interior points using the permutation matrix,
200  // so the points are ordered anticlockwise.
201  for(unsigned int k=0; k<vertex.size(); ++k){
202 
203  map.push_back(vertex[k]);
204  }
205 
206  for(unsigned int k=0; k<iEdge_1.size(); ++k){
207 
208  map.push_back(iEdge_1[k]);
209  }
210 
211  for(unsigned int k=0; k<iEdge_2.size(); ++k){
212 
213  map.push_back(iEdge_2[k]);
214  }
215 
216  for(unsigned int k=0; k<iEdge_3.size(); ++k){
217 
218  map.push_back(iEdge_3[k]);
219  }
220 
221  for(unsigned int k=0; k<interiorPoints.size(); ++k){
222 
223  map.push_back(interiorPoints[k]);
224  }
225 
226 
227  Array<OneD,NekDouble> points[2];
228  points[0] = Array<OneD,NekDouble>(GetTotNumPoints());
229  points[1] = Array<OneD,NekDouble>(GetTotNumPoints());
230  for(unsigned int index=0; index<map.size(); ++index){
231  points[0][index] = m_points[0][index];
232  points[1][index] = m_points[1][index];
233  }
234 
235  for(unsigned int index=0; index<map.size(); ++index){
236  m_points[0][index] = points[0][map[index]];
237  m_points[1][index] = points[1][map[index]];
238  }
239 
240  }
241 
242  } // end of namespace
243 } // end of namespace
244 
unsigned int GetTotNumPoints() const
Definition: Points.h:251
Array< OneD, DataType > m_weights
Definition: Points.h:350
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void CalculateInterpMatrix(const Array< OneD, const NekDouble > &xi, const Array< OneD, const NekDouble > &yi, Array< OneD, NekDouble > &interp)
NekVector< NekDouble > MakeQuadratureWeights(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y)
Definition: NodalUtil.cpp:579
static std::string npts
Definition: InputFld.cpp:43
Points< NekDouble >::MatrixSharedPtrType GetXDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi)
Definition: NodalUtil.cpp:852
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
unsigned int GetNumPoints() const
Definition: Points.h:246
unsigned int GetRows() const
Definition: NekVector.cpp:218
NekMatrix< NekDouble > GetInterpolationMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi)
Definition: NodalUtil.cpp:609
Array< OneD, DataType > m_points[3]
Definition: Points.h:349
MatrixSharedPtrType m_derivmatrix[3]
Definition: Points.h:351
Points< NekDouble >::MatrixSharedPtrType GetYDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi)
Definition: NodalUtil.cpp:1303
static boost::shared_ptr< PointsBaseType > Create(const PointsKey &key)
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230