Intrepid
Intrepid_HGRAD_TRI_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
1#ifndef INTREPID_HGRAD_TRI_CN_FEM_ORTHDEF_HPP
2#define INTREPID_HGRAD_TRI_CN_FEM_ORTHDEF_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
52namespace Intrepid {
53
54template<class Scalar, class ArrayScalar>
56{
57 this -> basisCardinality_ = (degree+1)*(degree+2)/2;
58 this -> basisDegree_ = degree;
59 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
60 this -> basisType_ = BASIS_FEM_HIERARCHICAL;
61 this -> basisCoordinates_ = COORDINATES_CARTESIAN;
62 this -> basisTagsAreSet_ = false;
63}
64
65
66
67template<class Scalar, class ArrayScalar>
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 = new int[tagSize * this->getCardinality()];
78 for (int i=0;i<this->getCardinality();i++) {
79 tags[4*i] = 2;
80 tags[4*i+1] = 0;
81 tags[4*i+2] = i;
82 tags[4*i+3] = this->getCardinality();
83 }
84
85 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
86 Intrepid::setOrdinalTagData(this -> tagToOrdinal_,
87 this -> ordinalToTag_,
88 tags,
89 this -> basisCardinality_,
90 tagSize,
91 posScDim,
92 posScOrd,
93 posDfOrd);
94}
95
96
98template<class Scalar, class ArrayScalar>
100 const ArrayScalar & inputPoints,
101 const EOperator operatorType) const {
102
103 // Verify arguments
104#ifdef HAVE_INTREPID_DEBUG
105 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
106 inputPoints,
107 operatorType,
108 this -> getBaseCellTopology(),
109 this -> getCardinality() );
110#endif
111 const int deg = this->getDegree();
112
113 // add more here and put in appropriate extra case statements below to enable higher derivatives.
114 void (*tabulators[])(ArrayScalar &, const int, const ArrayScalar &)
118
119
120 switch (operatorType) {
121 case OPERATOR_VALUE:
122 tabulators[0]( outputValues , deg , inputPoints );
123 break;
124 case OPERATOR_GRAD:
125 tabulators[1]( outputValues , deg , inputPoints );
126 break;
127 case OPERATOR_D1:
128 case OPERATOR_D2:
129 // add more case OPEATOR_Dn statements if you've added more items to the
130 // array above.
131 tabulators[operatorType-OPERATOR_D1+1]( outputValues , deg , inputPoints );
132 break;
133 default:
134 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument,
135 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid or unsupported operator" );
136
137 }
138
139 return;
140}
141
142
143
144template<class Scalar, class ArrayScalar>
146 const ArrayScalar & inputPoints,
147 const ArrayScalar & cellVertices,
148 const EOperator operatorType) const {
149 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
150 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function");
151}
152
153
154
155template<typename Scalar, typename ArrayScalar>
156void TabulatorTri<Scalar,ArrayScalar,0>::tabulate(ArrayScalar &outputValues ,
157 const int deg ,
158 const ArrayScalar &z )
159{
160 const int np = z.dimension( 0 );
161
162 // each point needs to be transformed from Pavel's element
163 // z(i,0) --> (2.0 * z(i,0) - 1.0)
164 // z(i,1) --> (2.0 * z(i,1) - 1.0)
165
166 // set up constant term
168 int idx_curp1,idx_curm1;
169
170 // set D^{0,0} = 1.0
171 for (int i=0;i<np;i++) {
172 outputValues(idx_cur,i) = Scalar( 1.0 ) + z(i,0) - z(i,0) + z(i,1) - z(i,1);
173 }
174
175
176 if (deg > 0) {
177 Teuchos::Array<Scalar> f1(np),f2(np),f3(np);
178
179 for (int i=0;i<np;i++) {
180 f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0));
181 f2[i] = 0.5 * (1.0-(2.0*z(i,1)-1.0));
182 f3[i] = f2[i] * f2[i];
183 }
184
185 // set D^{1,0} = f1
187 for (int i=0;i<np;i++) {
188 outputValues(idx_cur,i) = f1[i];
189 }
190
191 // recurrence in p
192 for (int p=1;p<deg;p++) {
196 Scalar a = (2.0*p+1.0)/(1.0+p);
197 Scalar b = p / (p+1.0);
198
199 for (int i=0;i<np;i++) {
200 outputValues(idx_curp1,i) = a * f1[i] * outputValues(idx_cur,i)
201 - b * f3[i] * outputValues(idx_curm1,i);
202 }
203 }
204
205 // D^{p,1}
206 for (int p=0;p<deg;p++) {
209 for (int i=0;i<np;i++) {
210 outputValues(idxp1,i) = outputValues(idxp0,i)
211 *0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
212 }
213 }
214
215
216 // recurrence in q
217 for (int p=0;p<deg-1;p++) {
218 for (int q=1;q<deg-p;q++) {
222 Scalar a,b,c;
223 TabulatorTri<Scalar,ArrayScalar,0>::jrc((Scalar)(2*p+1),(Scalar)0,q,a,b,c);
224 for (int i=0;i<np;i++) {
225 outputValues(idxpqp1,i)
226 = (a*(2.0*z(i,1)-1.0)+b)*outputValues(idxpq,i)
227 - c*outputValues(idxpqm1,i);
228 }
229 }
230 }
231 }
232
233 // orthogonalize
234 for (int p=0;p<=deg;p++) {
235 for (int q=0;q<=deg-p;q++) {
236 for (int i=0;i<np;i++) {
237 outputValues(TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q),i) *= sqrt( (p+0.5)*(p+q+1.0));
238 }
239 }
240 }
241
242 return;
243}
244
245
246
247template<typename Scalar, typename ArrayScalar>
248void TabulatorTri<Scalar,ArrayScalar,1>::tabulate(ArrayScalar &outputValues ,
249 const int deg ,
250 const ArrayScalar &z )
251{
252 const int np = z.dimension(0);
253 const int card = outputValues.dimension(0);
254 FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) );
255 for (int i=0;i<np;i++) {
256 for (int j=0;j<2;j++) {
257 dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) );
258 dZ(i,j).diff(j,2);
259 }
260 }
262
264 deg ,
265 dZ );
266
267 for (int i=0;i<card;i++) {
268 for (int j=0;j<np;j++) {
269 for (int k=0;k<2;k++) {
270 outputValues(i,j,k) = dResult(i,j).dx(k);
271 }
272 }
273 }
274
275 return;
276
277}
278
279
280
281template<typename Scalar, typename ArrayScalar, unsigned derivOrder>
283 const int deg ,
284 const ArrayScalar &z )
285{
286 const int np = z.dimension(0);
287 const int card = outputValues.dimension(0);
288 FieldContainer<Sacado::Fad::DFad<Scalar> > dZ( z.dimension(0) , z.dimension(1) );
289 for (int i=0;i<np;i++) {
290 for (int j=0;j<2;j++) {
291 dZ(i,j) = Sacado::Fad::DFad<Scalar>( z(i,j) );
292 dZ(i,j).diff(j,2);
293 }
294 }
295 FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np,derivOrder+1);
296
298 deg ,
299 dZ );
300
301 for (int i=0;i<card;i++) {
302 for (int j=0;j<np;j++) {
303 outputValues(i,j,0) = dResult(i,j,0).dx(0);
304 for (unsigned k=0;k<derivOrder;k++) {
305 outputValues(i,j,k+1) = dResult(i,j,k).dx(1);
306 }
307 }
308 }
309
310 return;
311
312
313}
314
315
316}// namespace Intrepid
317#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.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Triangle cell.
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
This is an internal class with a static member function for tabulating derivatives of orthogonal expa...
static void tabulate(ArrayScalar &outputValues, const int deg, const ArrayScalar &inputPoints)
basic tabulate mathod evaluates the derivOrder^th derivatives of the basis functions at inputPoints i...