Intrepid
Intrepid_HGRAD_PYR_C1_FEMDef.hpp
Go to the documentation of this file.
1 #ifndef INTREPID_HGRAD_PYR_C1_FEMDEF_HPP
2 #define INTREPID_HGRAD_PYR_C1_FEMDEF_HPP
3 
4 #include <limits>
5 
6 // @HEADER
7 // ************************************************************************
8 //
9 // Intrepid Package
10 // Copyright (2007) Sandia Corporation
11 //
12 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
13 // license for use of this work by or on behalf of the U.S. Government.
14 //
15 // Redistribution and use in source and binary forms, with or without
16 // modification, are permitted provided that the following conditions are
17 // met:
18 //
19 // 1. Redistributions of source code must retain the above copyright
20 // notice, this list of conditions and the following disclaimer.
21 //
22 // 2. Redistributions in binary form must reproduce the above copyright
23 // notice, this list of conditions and the following disclaimer in the
24 // documentation and/or other materials provided with the distribution.
25 //
26 // 3. Neither the name of the Corporation nor the names of the
27 // contributors may be used to endorse or promote products derived from
28 // this software without specific prior written permission.
29 //
30 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
31 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
32 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
33 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
34 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
35 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
36 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
37 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
38 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
39 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
40 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
41 //
42 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
43 // Denis Ridzal (dridzal@sandia.gov), or
44 // Kara Peterson (kjpeter@sandia.gov)
45 //
46 // ************************************************************************
47 // @HEADER
48 
54 namespace Intrepid {
55 
56  template<class Scalar, class ArrayScalar>
58  {
59  this -> basisCardinality_ = 5;
60  this -> basisDegree_ = 1;
61  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() );
62  this -> basisType_ = BASIS_FEM_DEFAULT;
63  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
64  this -> basisTagsAreSet_ = false;
65  }
66 
67 
68 template<class Scalar, class ArrayScalar>
70 
71  // Basis-dependent intializations
72  int tagSize = 4; // size of DoF tag
73  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
74  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
75  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
76 
77  // An array with local DoF tags assigned to basis functions, in the order of their local enumeration
78  int tags[] = { 0, 0, 0, 1,
79  0, 1, 0, 1,
80  0, 2, 0, 1,
81  0, 3, 0, 1,
82  0, 4, 0, 1 };
83 
84  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
85  Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
86  this -> ordinalToTag_,
87  tags,
88  this -> basisCardinality_,
89  tagSize,
90  posScDim,
91  posScOrd,
92  posDfOrd);
93 }
94 
95 
96 
97 template<class Scalar, class ArrayScalar>
99  const ArrayScalar & inputPoints,
100  const EOperator operatorType) const {
101 
102  // Verify arguments
103 #ifdef HAVE_INTREPID_DEBUG
104  Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
105  inputPoints,
106  operatorType,
107  this -> getBaseCellTopology(),
108  this -> getCardinality() );
109 #endif
110 
111  // Number of evaluation points = dim 0 of inputPoints
112  int dim0 = inputPoints.dimension(0);
113 
114  // Temporaries: (x,y,z) coordinates of the evaluation point
115  Scalar x = 0.0;
116  Scalar y = 0.0;
117  Scalar z = 0.0;
118  const Scalar eps = std::numeric_limits<Scalar>::epsilon( );
119 
120  switch (operatorType) {
121 
122  case OPERATOR_VALUE:
123  for (int i0 = 0; i0 < dim0; i0++) {
124  x = inputPoints(i0, 0);
125  y = inputPoints(i0, 1);
126  z = inputPoints(i0, 2);
127 
128  //be sure that the basis functions are defined when z is very close to 1.
129  if(fabs(z-1.0) < eps) {
130  if(z <= 1.0) z = 1.0-eps;
131  else z = 1.0+eps;
132  }
133 
134 
135  Scalar zTerm = 0.25/(1.0 - z);
136 
137  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
138  outputValues(0, i0) = (1.0 - x - z) * (1.0 - y - z) * zTerm;
139  outputValues(1, i0) = (1.0 + x - z) * (1.0 - y - z) * zTerm;
140  outputValues(2, i0) = (1.0 + x - z) * (1.0 + y - z) * zTerm;
141  outputValues(3, i0) = (1.0 - x - z) * (1.0 + y - z) * zTerm;
142  outputValues(4, i0) = z;
143  }
144  break;
145 
146  case OPERATOR_GRAD:
147  case OPERATOR_D1:
148  for (int i0 = 0; i0 < dim0; i0++) {
149 
150  x = inputPoints(i0, 0);
151  y = inputPoints(i0, 1);
152  z = inputPoints(i0, 2);
153 
154 
155  //be sure that the basis functions are defined when z is very close to 1.
156  //warning, the derivatives are discontinuous in (0, 0, 1)
157  if(fabs(z-1.0) < eps) {
158  if(z <= 1.0) z = 1.0-eps;
159  else z = 1.0+eps;
160  }
161 
162 
163  Scalar zTerm = 0.25/(1.0 - z);
164  Scalar zTerm2 = 4.0 * zTerm * zTerm;
165 
166  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
167  outputValues(0, i0, 0) = (y + z - 1.0) * zTerm;
168  outputValues(0, i0, 1) = (x + z - 1.0) * zTerm;
169  outputValues(0, i0, 2) = x * y * zTerm2 - 0.25;
170 
171  outputValues(1, i0, 0) = (1.0 - y - z) * zTerm;
172  outputValues(1, i0, 1) = (z - x - 1.0) * zTerm;
173  outputValues(1, i0, 2) = - x*y * zTerm2 - 0.25;
174 
175  outputValues(2, i0, 0) = (1.0 + y - z) * zTerm;
176  outputValues(2, i0, 1) = (1.0 + x - z) * zTerm;
177  outputValues(2, i0, 2) = x * y * zTerm2 - 0.25;
178 
179  outputValues(3, i0, 0) = (z - y - 1.0) * zTerm;
180  outputValues(3, i0, 1) = (1.0 - x - z) * zTerm;
181  outputValues(3, i0, 2) = - x*y * zTerm2 - 0.25;
182 
183  outputValues(4, i0, 0) = 0.0;
184  outputValues(4, i0, 1) = 0.0;
185  outputValues(4, i0, 2) = 1;
186  }
187  break;
188 
189  case OPERATOR_CURL:
190  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
191  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D");
192  break;
193 
194  case OPERATOR_DIV:
195  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
196  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D");
197  break;
198 
199  case OPERATOR_D2:
200  for (int i0 = 0; i0 < dim0; i0++) {
201  x = inputPoints(i0,0);
202  y = inputPoints(i0,1);
203  z = inputPoints(i0,2);
204 
205  //be sure that the basis functions are defined when z is very close to 1.
206  //warning, the derivatives are discontinuous in (0, 0, 1)
207  if(fabs(z-1.0) < eps) {
208  if(z <= 1.0) z = 1.0-eps;
209  else z = 1.0+eps;
210  }
211 
212 
213  Scalar zTerm = 0.25/(1.0 - z);
214  Scalar zTerm2 = 4.0 * zTerm * zTerm;
215  Scalar zTerm3 = 8.0 * zTerm * zTerm2;
216 
217  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6)
218  outputValues(0, i0, 0) = 0.0; // {2, 0, 0}
219  outputValues(0, i0, 1) = zTerm; // {1, 1, 0}
220  outputValues(0, i0, 2) = y*zTerm2; // {1, 0, 1}
221  outputValues(0, i0, 3) = 0.0; // {0, 2, 0}
222  outputValues(0, i0, 4) = x*zTerm2; // {0, 1, 1}
223  outputValues(0, i0, 5) = x*y*zTerm3; // {0, 0, 2}
224 
225  outputValues(1, i0, 0) = 0.0; // {2, 0, 0}
226  outputValues(1, i0, 1) = -zTerm; // {1, 1, 0}
227  outputValues(1, i0, 2) = -y*zTerm2; // {1, 0, 1}
228  outputValues(1, i0, 3) = 0.0; // {0, 2, 0}
229  outputValues(1, i0, 4) = -x*zTerm2; // {0, 1, 1}
230  outputValues(1, i0, 5) = -x*y*zTerm3; // {0, 0, 2}
231 
232  outputValues(2, i0, 0) = 0.0; // {2, 0, 0}
233  outputValues(2, i0, 1) = zTerm; // {1, 1, 0}
234  outputValues(2, i0, 2) = y*zTerm2; // {1, 0, 1}
235  outputValues(2, i0, 3) = 0.0; // {0, 2, 0}
236  outputValues(2, i0, 4) = x*zTerm2; // {0, 1, 1}
237  outputValues(2, i0, 5) = x*y*zTerm3; // {0, 0, 2}
238 
239  outputValues(3, i0, 0) = 0.0; // {2, 0, 0}
240  outputValues(3, i0, 1) = -zTerm; // {1, 1, 0}
241  outputValues(3, i0, 2) = -y*zTerm2; // {1, 0, 1}
242  outputValues(3, i0, 3) = 0.0; // {0, 2, 0}
243  outputValues(3, i0, 4) = -x*zTerm2; // {0, 1, 1}
244  outputValues(3, i0, 5) = -x*y*zTerm3; // {0, 0, 2}
245 
246  outputValues(4, i0, 0) = 0.0; // {2, 0, 0}
247  outputValues(4, i0, 1) = 0.0; // {1, 1, 0}
248  outputValues(4, i0, 2) = 0.0; // {1, 0, 1}
249  outputValues(4, i0, 3) = 0.0; // {0, 2, 0}
250  outputValues(4, i0, 4) = 0.0; // {0, 1, 1}
251  outputValues(4, i0, 5) = 0.0; // {0, 0, 2}
252  }
253  break;
254 
255  case OPERATOR_D3:
256  case OPERATOR_D4:
257  case OPERATOR_D5:
258  case OPERATOR_D6:
259  case OPERATOR_D7:
260  case OPERATOR_D8:
261  case OPERATOR_D9:
262  case OPERATOR_D10:
263  {
264  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
265  int DkCardinality = Intrepid::getDkCardinality(operatorType,
266  this -> basisCellTopology_.getDimension() );
267  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
268  for (int i0 = 0; i0 < dim0; i0++) {
269  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
270  outputValues(dofOrd, i0, dkOrd) = 0.0;
271  }
272  }
273  }
274  }
275  break;
276  default:
277  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
278  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): Invalid operator type");
279  }
280 }
281 
282 
283 
284 template<class Scalar, class ArrayScalar>
286  const ArrayScalar & inputPoints,
287  const ArrayScalar & cellVertices,
288  const EOperator operatorType) const {
289  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
290  ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): FEM Basis calling an FVD member function");
291 }
292 }// namespace Intrepid
293 #endif
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
int isValidOperator(const EOperator operatorType)
Verifies validity of an operator enum.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Pyramid cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.