Nektar++
NodalTetEvenlySpaced.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalTetEvenlySpaced.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: 3D Nodal Tetrahedron 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 tetrahedron
53  // edges and vertices are ordered as anticlockwise
54 
55  bool isVertex(int x, int y, int z, int npts){
56  return (x==0 && y==0 && z==0) || (x==(npts-1) && y==0 && z==0) || (x==0 && y==(npts-1) && z==0) || (x==0 && y==0 && z==(npts-1));
57  }
58 
59  bool isEdge_01(int x, int y, int z, int npts){ // edge 0
60  return y==0 && z==0;
61  }
62 
63  bool isEdge_12(int x, int y, int z, int npts){ // edge 1
64  return z==0 && x + y == npts -1;
65  }
66 
67  bool isEdge_20(int x, int y, int z, int npts){ // edge 2
68  return x==0 && z==0;
69  }
70 
71  bool isEdge_03(int x, int y, int z, int npts){ // edge 3
72  return x==0 && y==0;
73  }
74 
75  bool isEdge_13(int x, int y, int z, int npts){ // edge 4
76  return y==0 && x + z == npts -1;
77  }
78 
79  bool isEdge_23(int x, int y, int z, int npts){ // edge 5
80  return x==0 && y + z == npts -1;
81  }
82 
83  bool isEdge(int x, int y, int z, int npts){
84  return isEdge_01(x,y,z,npts)||isEdge_12(x,y,z,npts)||isEdge_20(x,y,z,npts)
85  ||isEdge_03(x,y,z,npts)||isEdge_13(x,y,z,npts)||isEdge_23(x,y,z,npts);
86  }
87 
88  bool isFace_012(int x, int y, int z, int npts){ // bottom face (face 0)
89  return z==0;
90  }
91 
92  bool isFace_013(int x, int y, int z, int npts){ // face 1
93  return y==0;
94  }
95 
96  bool isFace_123(int x, int y, int z, int npts){ // face 2
97  return x+y+z==npts-1;
98  }
99 
100  bool isFace_203(int x, int y, int z, int npts){ // face 3
101  return x==0;
102  }
103 
104  bool isFace(int x, int y, int z, int npts){
105  return isFace_012(x, y, z, npts) || isFace_013(x, y, z, npts)
106  || isFace_123(x, y, z, npts) || isFace_203(x, y, z, npts);
107  }
108  }
109 
110  // Calculate evenly spaced number of points
112  {
113  // Allocate the storage for points
115 
116  // Populate m_points
117  unsigned int npts = GetNumPoints();
118  NekDouble delta = 2.0/(npts - 1.0);
119  for(unsigned int z=0, index=0; z<npts; ++z){
120  for(int y=0; y<npts-z; ++y){
121  for(int x=0; x<npts-z-y; ++x, ++index){
122  NekDouble xi = -1.0 + x*delta;
123  NekDouble yi = -1.0 + y*delta;
124  NekDouble zi = -1.0 + z*delta;
125 
126  m_points[0][index] = xi;
127  m_points[1][index] = yi;
128  m_points[2][index] = zi;
129  }
130  }
131  }
132 
134  }
135 
137  {
138  unsigned int npts = GetNumPoints();
139  using std::vector;
140  vector<int> vertex;
141  vector<int> iEdge_01; // interior edge 0
142  vector<int> iEdge_12; // interior edge 1
143  vector<int> iEdge_20; // interior edge 2
144  vector<int> iEdge_03; // interior edge 3
145  vector<int> iEdge_13; // interior edge 4
146  vector<int> iEdge_23; // interior edge 5
147  vector<int> iFace_012; // interior face 0
148  vector<int> iFace_013; // interior face 1
149  vector<int> iFace_123; // interior face 2
150  vector<int> iFace_203; // interior face 3
151  vector<int> interiorVolumePoints; // interior volume points
152  vector<int> map;
153 
154  // Build the lattice tetrahedron left to right - bottom to top
155  for(int z=0, index=0; z<npts; ++z){
156  for(int y=0; y<npts-z; ++y){
157  for(int x=0; x<npts-z-y; ++x, ++index){
158 
159  if( isVertex(x,y,z,npts) ){ // vertex
160 
161  vertex.push_back(index);
162 
163  } else if( isEdge(x,y,z,npts) ){ // interior edge
164 
165  if( isEdge_01(x,y,z,npts) ){ // interior edge 0
166 
167  iEdge_01.push_back(index);
168 
169  } else if( isEdge_12(x,y,z,npts) ){ // interior edge 1
170 
171  iEdge_12.push_back(index);
172 
173  } else if( isEdge_20(x,y,z,npts) ){ // interior edge 2
174 
175  iEdge_20.insert(iEdge_20.begin(), index);
176 
177  } else if( isEdge_03(x,y,z,npts) ){ // interior edge 3
178 
179  iEdge_03.push_back(index);
180 
181  } else if( isEdge_13(x,y,z,npts) ){ // interior edge 4
182 
183  iEdge_13.push_back(index);
184 
185  } else if( isEdge_23(x,y,z,npts) ){ // interior edge 5
186 
187  iEdge_23.push_back(index);
188 
189  }
190 
191  } else if( isFace(x,y,z,npts) ) { // interior face
192 
193  if( isFace_012(x,y,z,npts) ){ // interior face 0
194 
195  iFace_012.push_back(index);
196 
197  } else if( isFace_013(x,y,z,npts) ){ // interior face 1
198 
199  iFace_013.push_back(index);
200 
201  } else if( isFace_123(x,y,z,npts) ){ // interior face 2
202 
203  iFace_123.push_back(index);
204 
205  } else if( isFace_203(x,y,z,npts) ){ // interior face 3
206 
207  iFace_203.push_back(index);
208 
209  }
210  } else { // interior volume points
211 
212  interiorVolumePoints.push_back(index);
213  }
214  }
215  }
216  }
217 
218  // Mapping the vertex, edges, faces, interior volume points using the permutation matrix,
219  // so the points are ordered anticlockwise.
220  for(unsigned int n=0; n<vertex.size(); ++n){
221 
222  map.push_back(vertex[n]);
223  }
224 
225  for(unsigned int n=0; n<iEdge_01.size(); ++n){
226 
227  map.push_back(iEdge_01[n]);
228  }
229 
230  for(unsigned int n=0; n<iEdge_12.size(); ++n){
231 
232  map.push_back(iEdge_12[n]);
233  }
234 
235  for(unsigned int n=0; n<iEdge_20.size(); ++n){
236 
237  map.push_back(iEdge_20[n]);
238  }
239 
240  for(unsigned int n=0; n<iEdge_03.size(); ++n){
241 
242  map.push_back(iEdge_03[n]);
243  }
244 
245  for(unsigned int n=0; n<iEdge_13.size(); ++n){
246 
247  map.push_back(iEdge_13[n]);
248  }
249 
250  for(unsigned int n=0; n<iEdge_23.size(); ++n){
251 
252  map.push_back(iEdge_23[n]);
253  }
254 
255  for(unsigned int n=0; n<iFace_012.size(); ++n){
256 
257  map.push_back(iFace_012[n]);
258  }
259 
260  for(unsigned int n=0; n<iFace_013.size(); ++n){
261 
262  map.push_back(iFace_013[n]);
263  }
264 
265  for(unsigned int n=0; n<iFace_123.size(); ++n){
266 
267  map.push_back(iFace_123[n]);
268  }
269 
270  for(unsigned int n=0; n<iFace_203.size(); ++n){
271 
272  map.push_back(iFace_203[n]);
273  }
274 
275  for(unsigned int n=0; n<interiorVolumePoints.size(); ++n){
276 
277  map.push_back(interiorVolumePoints[n]);
278  }
279 
280 
281  Array<OneD, NekDouble> points[3];
285  for(unsigned int index=0; index<map.size(); ++index){
286 
287  points[0][index] = m_points[0][index];
288  points[1][index] = m_points[1][index];
289  points[2][index] = m_points[2][index];
290 
291  }
292 
293  for(unsigned int index=0; index<map.size(); ++index){
294 
295  m_points[0][index] = points[0][map[index]];
296  m_points[1][index] = points[1][map[index]];
297  m_points[2][index] = points[2][map[index]];
298 
299  }
300 
301  }
302 
303 
304 
306  {
307  // Allocate storage for points
309 
310  typedef DataType T;
311 
312  // Solve the Vandermonde system of integrals for the weight vector
314 
315  m_weights = Array<OneD,T>( w.GetRows(), w.GetPtr() );
316 
317  }
318 
319 
320  // ////////////////////////////////////////
321  // CalculateInterpMatrix()
323  const Array<OneD, const NekDouble>& yia,
324  const Array<OneD, const NekDouble>& zia,
325  Array<OneD, NekDouble>& interp)
326  {
330  NekVector<NekDouble> xi( xia );
331  NekVector<NekDouble> yi( yia );
332  NekVector<NekDouble> zi( zia );
333  NekMatrix<NekDouble> interMat = GetTetInterpolationMatrix(x, y, z, xi, yi, zi);
334 
335  int rows = xi.GetRows(), cols = GetTotNumPoints();
336  for( int i = 0; i < rows; ++i ) {
337  for( int j = 0; j < cols; ++j ) {
338  interp[j + i*cols] = interMat(i,j);
339  }
340  }
341 
342  }
343 
344  // ////////////////////////////////////////
345  // CalculateDerivMatrix()
347  {
348 
349  // Allocate the derivative matrix.
351 
355  NekVector<NekDouble> xi = x;
356  NekVector<NekDouble> yi = y;
357  NekVector<NekDouble> zi = z;
358 
359  *m_derivmatrix[0] = *GetTetXDerivativeMatrix(x,y,z,xi,yi,zi);
360 
361  *m_derivmatrix[1] = *GetTetYDerivativeMatrix(x,y,z,xi,yi,zi);
362 
363  *m_derivmatrix[2] = *GetTetZDerivativeMatrix(x,y,z,xi,yi,zi);
364 
365  }
366 
367  boost::shared_ptr<PointsBaseType> NodalTetEvenlySpaced::Create(const PointsKey &key)
368  {
369  boost::shared_ptr<PointsBaseType> returnval(MemoryManager<NodalTetEvenlySpaced>::AllocateSharedPtr(key));
370 
371  returnval->Initialize();
372 
373  return returnval;
374  }
375 
376 
377 
378  } // end of namespace
379 } // end of namespace
380 
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, const Array< OneD, const NekDouble > &zi, Array< OneD, NekDouble > &interp)
NekMatrix< NekDouble > GetTetInterpolationMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
Definition: NodalUtil.cpp:594
NekVector< NekDouble > MakeTetWeights(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z)
Definition: NodalUtil.cpp:562
static boost::shared_ptr< PointsBaseType > Create(const PointsKey &key)
static std::string npts
Definition: InputFld.cpp:43
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
Points< NekDouble >::MatrixSharedPtrType GetTetZDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
Definition: NodalUtil.cpp:1195
Points< NekDouble >::MatrixSharedPtrType GetTetXDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
Definition: NodalUtil.cpp:806
unsigned int GetNumPoints() const
Definition: Points.h:246
unsigned int GetRows() const
Definition: NekVector.cpp:218
Array< OneD, DataType > m_points[3]
Definition: Points.h:349
MatrixSharedPtrType m_derivmatrix[3]
Definition: Points.h:351
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
Points< NekDouble >::MatrixSharedPtrType GetTetYDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
Definition: NodalUtil.cpp:1047