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 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...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Wedge cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
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.