Intrepid
Intrepid_HGRAD_TRI_C2_FEMDef.hpp
1 #ifndef INTREPID_HGRAD_TRI_C2_FEMDEF_HPP
2 #define INTREPID_HGRAD_TRI_C2_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 
54 template<class Scalar, class ArrayScalar>
55 Basis_HGRAD_TRI_C2_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TRI_C2_FEM()
56  {
57  this -> basisCardinality_ = 6;
58  this -> basisDegree_ = 2;
59  this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
60  this -> basisType_ = BASIS_FEM_DEFAULT;
61  this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62  this -> basisTagsAreSet_ = false;
63  }
64 
65 
66 
67 template<class Scalar, class ArrayScalar>
68 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::initializeTags() {
69 
70  // Basis-dependent initializations
71  int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
72  int posScDim = 0; // position in the tag, counting from 0, of the subcell dim
73  int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
74  int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
75 
76  // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
77  int tags[] = { 0, 0, 0, 1,
78  0, 1, 0, 1,
79  0, 2, 0, 1,
80  1, 0, 0, 1,
81  1, 1, 0, 1,
82  1, 2, 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>
98 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues,
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) coordinates of the evaluation point
115  Scalar x = 0.0;
116  Scalar y = 0.0;
117 
118  switch (operatorType) {
119 
120  case OPERATOR_VALUE:
121  for (int i0 = 0; i0 < dim0; i0++) {
122  x = inputPoints(i0, 0);
123  y = inputPoints(i0, 1);
124 
125  // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0)
126  outputValues(0, i0) = (x + y - 1.0)*(2.0*x + 2.0*y - 1.0);
127  outputValues(1, i0) = x*(2.0*x - 1.0);
128  outputValues(2, i0) = y*(2.0*y - 1.0);
129  outputValues(3, i0) = -4.0*x*(x + y - 1.0);
130  outputValues(4, i0) = 4.0*x*y;
131  outputValues(5, i0) = -4.0*y*(x + y - 1.0);
132 
133  }
134  break;
135 
136  case OPERATOR_GRAD:
137  case OPERATOR_D1:
138  for (int i0 = 0; i0 < dim0; i0++) {
139  x = inputPoints(i0, 0);
140  y = inputPoints(i0, 1);
141 
142  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim)
143  outputValues(0, i0, 0) = 4.0*x + 4.0*y - 3.0;
144  outputValues(0, i0, 1) = 4.0*x + 4.0*y - 3.0;
145 
146  outputValues(1, i0, 0) = 4.0*x - 1.0;
147  outputValues(1, i0, 1) = 0.0;
148 
149  outputValues(2, i0, 0) = 0.0;
150  outputValues(2, i0, 1) = 4.0*y - 1.0;
151 
152  outputValues(3, i0, 0) = -4.0*(2.0*x + y - 1.0);
153  outputValues(3, i0, 1) = -4.0*x;
154 
155  outputValues(4, i0, 0) = 4.0*y;
156  outputValues(4, i0, 1) = 4.0*x;
157 
158  outputValues(5, i0, 0) = -4.0*y;
159  outputValues(5, i0, 1) = -4.0*(x + 2.0*y - 1.0);
160  }
161  break;
162 
163  case OPERATOR_CURL:
164  for (int i0 = 0; i0 < dim0; i0++) {
165  x = inputPoints(i0, 0);
166  y = inputPoints(i0, 1);
167 
168  // CURL(u) = (u_y, -u_x), is rotated GRAD
169  outputValues(0, i0, 1) =-(4.0*x + 4.0*y - 3.0);
170  outputValues(0, i0, 0) = 4.0*x + 4.0*y - 3.0;
171 
172  outputValues(1, i0, 1) =-(4.0*x - 1.0);
173  outputValues(1, i0, 0) = 0.0;
174 
175  outputValues(2, i0, 1) = 0.0;
176  outputValues(2, i0, 0) = 4.0*y - 1.0;
177 
178  outputValues(3, i0, 1) = 4.0*(2.0*x + y - 1.0);
179  outputValues(3, i0, 0) = -4.0*x;
180 
181  outputValues(4, i0, 1) = -4.0*y;
182  outputValues(4, i0, 0) = 4.0*x;
183 
184  outputValues(5, i0, 1) = 4.0*y;
185  outputValues(5, i0, 0) = -4.0*(x + 2.0*y - 1.0);
186  }
187  break;
188 
189  case OPERATOR_DIV:
190  TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
191  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): DIV is invalid operator for rank-0 (scalar) fields in 2D.");
192  break;
193 
194  case OPERATOR_D2:
195  for (int i0 = 0; i0 < dim0; i0++) {
196  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
197  // D2 -> (2,0) -> dx^2.
198  outputValues(0, i0, 0) = 4.0;
199  outputValues(1, i0, 0) = 4.0;
200  outputValues(2, i0, 0) = 0.0;
201  outputValues(3, i0, 0) =-8.0;
202  outputValues(4, i0, 0) = 0.0;
203  outputValues(5, i0, 0) = 0.0;
204 
205  // D2 -> (1,1) -> dx dy
206  outputValues(0, i0, 1) = 4.0;
207  outputValues(1, i0, 1) = 0.0;
208  outputValues(2, i0, 1) = 0.0;
209  outputValues(3, i0, 1) =-4.0;
210  outputValues(4, i0, 1) = 4.0;
211  outputValues(5, i0, 1) =-4.0;
212 
213  // D2 -> (0,2) -> dy^2
214  outputValues(0, i0, 2) = 4.0;
215  outputValues(1, i0, 2) = 0.0;
216  outputValues(2, i0, 2) = 4.0;
217  outputValues(3, i0, 2) = 0.0;
218  outputValues(4, i0, 2) = 0.0;
219  outputValues(5, i0, 2) =-8.0;
220  }// for i0
221  break;
222 
223  case OPERATOR_D3:
224  case OPERATOR_D4:
225  case OPERATOR_D5:
226  case OPERATOR_D6:
227  case OPERATOR_D7:
228  case OPERATOR_D8:
229  case OPERATOR_D9:
230  case OPERATOR_D10:
231  {
232  // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality)
233  int DkCardinality = Intrepid::getDkCardinality(operatorType,
234  this -> basisCellTopology_.getDimension() );
235  for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) {
236  for (int i0 = 0; i0 < dim0; i0++) {
237  for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){
238  outputValues(dofOrd, i0, dkOrd) = 0.0;
239  }
240  }
241  }
242  }
243  break;
244 
245  default:
246  TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument,
247  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): Invalid operator type");
248  }
249 }
250 
251 
252 
253 template<class Scalar, class ArrayScalar>
254 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues,
255  const ArrayScalar & inputPoints,
256  const ArrayScalar & cellVertices,
257  const EOperator operatorType) const {
258  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
259  ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): FEM Basis calling an FVD member function");
260 }
261 
262 template<class Scalar, class ArrayScalar>
263 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getDofCoords(ArrayScalar & DofCoords) const {
264 #ifdef HAVE_INTREPID_DEBUG
265  // Verify rank of output array.
266  TEUCHOS_TEST_FOR_EXCEPTION( !(DofCoords.rank() == 2), std::invalid_argument,
267  ">>> ERROR: (Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords) rank = 2 required for DofCoords array");
268  // Verify 0th dimension of output array.
269  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(0) == this -> basisCardinality_ ), std::invalid_argument,
270  ">>> ERROR: (Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords) mismatch in number of DoF and 0th dimension of DofCoords array");
271  // Verify 1st dimension of output array.
272  TEUCHOS_TEST_FOR_EXCEPTION( !( DofCoords.dimension(1) == (int)(this -> basisCellTopology_.getDimension()) ), std::invalid_argument,
273  ">>> ERROR: (Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in DofCoords array");
274 #endif
275 
276  DofCoords(0,0) = 0.0; DofCoords(0,1) = 0.0;
277  DofCoords(1,0) = 1.0; DofCoords(1,1) = 0.0;
278  DofCoords(2,0) = 0.0; DofCoords(2,1) = 1.0;
279  DofCoords(3,0) = 0.5; DofCoords(3,1) = 0.0;
280  DofCoords(4,0) = 0.5; DofCoords(4,1) = 0.5;
281  DofCoords(5,0) = 0.0; DofCoords(5,1) = 0.5;
282 }
283 
284 }// namespace Intrepid
285 #endif
286 
287 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
288 #ifdef __GNUC__
289 #warning "The Intrepid package is deprecated"
290 #endif
291 #endif
292 
Intrepid::Basis_HGRAD_TRI_C2_FEM::initializeTags
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Definition: Intrepid_HGRAD_TRI_C2_FEMDef.hpp:68
Intrepid::Basis_HGRAD_TRI_C2_FEM::getValues
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
Definition: Intrepid_HGRAD_TRI_C2_FEMDef.hpp:98
Intrepid::Basis_HGRAD_TRI_C2_FEM::Basis_HGRAD_TRI_C2_FEM
Basis_HGRAD_TRI_C2_FEM()
Constructor.
Definition: Intrepid_HGRAD_TRI_C2_FEMDef.hpp:55
Intrepid::Basis_HGRAD_TRI_C2_FEM::getDofCoords
void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference Quadrilateral.
Definition: Intrepid_HGRAD_TRI_C2_FEMDef.hpp:263
Generated by   doxygen 1.8.5