Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ShapeType.hpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ShapeType.h
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 ////////////////////////////////////////////////////////////////////////////////
33 
34 
35 #ifndef NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHAPE_TYPE_H
36 #define NEKTAR_LIB_UTILITIES_BASIC_UTILS_SHAPE_TYPE_H
37 
38 #include <algorithm>
39 
41 
42 #ifdef min
43 #undef min
44 #endif
45 
46 namespace Nektar
47 {
48  namespace LibUtilities
49  {
50 
51  // Types of geometry types.
52  enum ShapeType
53  {
64  };
65 
66  const char* const ShapeTypeMap[] =
67  {
68  "NoGeomShapeType",
69  "Point",
70  "Segment",
71  "Triangle",
72  "Quadrilateral",
73  "Tetrahedron",
74  "Pyramid",
75  "Prism",
76  "Hexahedron"
77  };
78 
79 
80  // Hold the dimension of each of the types of shapes.
81  const unsigned int ShapeTypeDimMap[SIZE_ShapeType] =
82  {
83  0, // Unknown
84  0, // ePoint
85  1, // eSegment
86  2, // eTriangle
87  2, // eQuadrilateral
88  3, // eTetrahedron
89  3, // ePyramid
90  3, // ePrism
91  3, // eHexahedron
92  };
93 
94 
95  namespace StdSegData
96  {
97  inline int getNumberOfCoefficients(int Na)
98  {
99  return Na;
100  }
101 
102  inline int getNumberOfBndCoefficients(int Na)
103  {
104  return 2;
105  }
106  }
107 
108  // Dimensions of coefficients for each space
109  namespace StdTriData
110  {
111  inline int getNumberOfCoefficients(int Na, int Nb)
112  {
113  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
114  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
115  ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
116  "than order in 'b' direction");
117  return Na*(Na+1)/2 + Na*(Nb-Na);
118  }
119 
120  inline int getNumberOfBndCoefficients(int Na, int Nb)
121  {
122  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
123  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
124  ASSERTL1(Na <= Nb, "order in 'a' direction is higher "
125  "than order in 'b' direction");
126  return (Na-1) + 2*(Nb-1);
127  }
128  }
129 
130  namespace StdQuadData
131  {
132  inline int getNumberOfCoefficients(int Na, int Nb)
133  {
134  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
135  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
136  return Na*Nb;
137  }
138 
139  inline int getNumberOfBndCoefficients(int Na, int Nb)
140  {
141  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
142  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
143  return 2*(Na-1) + 2*(Nb-1);
144  }
145  }
146 
147 
148 
149  namespace StdHexData
150  {
151  inline int getNumberOfCoefficients( int Na, int Nb, int Nc )
152  {
153  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
154  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
155  ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
156  return Na*Nb*Nc;
157  }
158 
159  inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
160  {
161  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
162  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
163  ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
164  return 2*Na*Nb + 2*Na*Nc + 2*Nb*Nc
165  - 4*(Na + Nb + Nc) + 8;
166  }
167  }
168 
169  namespace StdTetData
170  {
171  /**
172  * Adds up the number of cells in a truncated Nc by Nc by Nc
173  * pyramid, where the longest Na rows and longest Nb columns are
174  * kept. Example: (Na, Nb, Nc) = (3, 4, 5); The number of
175  * coefficients is the sum of the elements of the following
176  * matrix:
177  *
178  * |5 4 3 2 0|
179  * |4 3 2 0 |
180  * |3 2 0 |
181  * |0 0 |
182  * |0 |
183  *
184  * Sum = 28 = number of tet coefficients.
185  */
186  inline int getNumberOfCoefficients(int Na, int Nb, int Nc)
187  {
188  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
189  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
190  ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
191  ASSERTL1(Na <= Nc, "order in 'a' direction is higher "
192  "than order in 'c' direction");
193  ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
194  "than order in 'c' direction");
195  int nCoef = 0;
196  for (int a = 0; a < Na; ++a)
197  {
198  for (int b = 0; b < Nb - a; ++b)
199  {
200  for (int c = 0; c < Nc - a - b; ++c)
201  {
202  ++nCoef;
203  }
204  }
205  }
206  return nCoef;
207  }
208 
209  inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
210  {
211  ASSERTL2(Na > 1, "Order in 'a' direction must be > 1.");
212  ASSERTL2(Nb > 1, "Order in 'b' direction must be > 1.");
213  ASSERTL2(Nc > 1, "Order in 'c' direction must be > 1.");
214  ASSERTL1(Na <= Nc, "order in 'a' direction is higher "
215  "than order in 'c' direction");
216  ASSERTL1(Nb <= Nc, "order in 'b' direction is higher "
217  "than order in 'c' direction");
218 
219  int nCoef = Na*(Na+1)/2 + (Nb-Na)*Na // base
220  + Na*(Na+1)/2 + (Nc-Na)*Na // front
221  + 2*(Nb*(Nb+1)/2 + (Nc-Nb)*Nb)// 2 other sides
222  - Na - 2*Nb - 3*Nc // less edges
223  + 4; // plus vertices
224 
225  return nCoef;
226  }
227  }
228 
229 
230  namespace StdPyrData
231  {
232  inline int getNumberOfCoefficients(int Na, int Nb, int Nc)
233  {
234  ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
235  ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
236  ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
237  ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
238  "than order in 'c' direction.");
239  ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
240  "than order in 'c' direction.");
241 
242  // Count number of coefficients explicitly.
243  const int Pi = Na - 2, Qi = Nb - 2, Ri = Nc - 2;
244  int nCoeff =
245  5 + // vertices
246  Pi * 2 + Qi * 2 + Ri * 4 + // base edges
247  Pi * Qi + // base quad
248  Pi * (2*Ri - Pi - 1) + // p-r triangles;
249  Qi * (2*Ri - Qi - 1); // q-r triangles;
250 
251  // Count number of interior tet modes
252  for (int a = 0; a < Pi - 1; ++a)
253  {
254  for (int b = 0; b < Qi - a - 1; ++b)
255  {
256  for (int c = 0; c < Ri - a - b -1; ++c)
257  {
258  ++nCoeff;
259  }
260  }
261  }
262 
263  return nCoeff;
264  }
265 
266  inline int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
267  {
268  ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
269  ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
270  ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
271  ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
272  "than order in 'c' direction.");
273  ASSERTL1(Nb <= Nc, "Order in 'b' direction is higher "
274  "than order in 'c' direction.");
275 
276  return Na*Nb // base
277  + 2*(Na*(Na+1)/2 + (Nc-Na)*Na) // front and back
278  + 2*(Nb*(Nb+1)/2 + (Nc-Nb)*Nb) // sides
279  - 2*Na - 2*Nb - 4*Nc // less edges
280  + 5; // plus vertices
281  }
282  }
283 
284  namespace StdPrismData
285  {
286  inline int getNumberOfCoefficients( int Na, int Nb, int Nc )
287  {
288  ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
289  ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
290  ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
291  ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
292  "than order in 'c' direction.");
293 
294  return Nb*StdTriData::getNumberOfCoefficients(Na,Nc);
295  }
296 
297  inline int getNumberOfBndCoefficients( int Na, int Nb, int Nc)
298  {
299  ASSERTL1(Na > 1, "Order in 'a' direction must be > 1.");
300  ASSERTL1(Nb > 1, "Order in 'b' direction must be > 1.");
301  ASSERTL1(Nc > 1, "Order in 'c' direction must be > 1.");
302  ASSERTL1(Na <= Nc, "Order in 'a' direction is higher "
303  "than order in 'c' direction.");
304 
305  return Na*Nb + 2*Nb*Nc // rect faces
306  + 2*( Na*(Na+1)/2 + (Nc-Na)*Na ) // tri faces
307  - 2*Na - 3*Nb - 4*Nc // less edges
308  + 6; // plus vertices
309  }
310  }
311 
312  inline int GetNumberOfCoefficients(ShapeType shape, std::vector<unsigned int> &modes, int offset)
313  {
314  int returnval = 0;
315  switch(shape)
316  {
317  case eSegment:
318  returnval = modes[offset];
319  break;
320  case eTriangle:
321  returnval = StdTriData::getNumberOfCoefficients(modes[offset],modes[offset+1]);
322  break;
323  case eQuadrilateral:
324  returnval = modes[offset]*modes[offset+1];
325  break;
326  case eTetrahedron:
327  returnval = StdTetData::getNumberOfCoefficients(modes[offset],modes[offset+1],modes[offset+2]);
328  break;
329  case ePyramid:
330  returnval = StdPyrData::getNumberOfCoefficients(modes[offset],modes[offset+1],modes[offset+2]);
331  break;
332  case ePrism:
333  returnval = StdPrismData::getNumberOfCoefficients(modes[offset],modes[offset+1],modes[offset+2]);
334  break;
335  case eHexahedron:
336  returnval = modes[offset]*modes[offset+1]*modes[offset+2];
337  break;
338  default:
339  ASSERTL0(false,"Unknown Shape Type");
340  break;
341  }
342 
343  return returnval;
344  }
345  }
346 }
347 
348 #endif
349