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