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 
52 namespace Intrepid {
53 
54 template<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 
67 template<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 
97 
98 template<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 
144 template<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 
155 template<typename Scalar, typename ArrayScalar>
156 void 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
167  int idx_cur = TabulatorTri<Scalar,ArrayScalar,0>::idx(0,0);
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++) {
194  idx_curp1 = TabulatorTri<Scalar,ArrayScalar,0>::idx(p+1,0);
195  idx_curm1 = TabulatorTri<Scalar,ArrayScalar,0>::idx(p-1,0);
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++) {
219  int idxpqp1=TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q+1);
221  int idxpqm1=TabulatorTri<Scalar,ArrayScalar,0>::idx(p,q-1);
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 
247 template<typename Scalar, typename ArrayScalar>
248 void 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  }
261  FieldContainer<Sacado::Fad::DFad<Scalar> > dResult(card,np);
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 
281 template<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
318 
319 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
320 #ifdef __GNUC__
321 #warning "The Intrepid package is deprecated"
322 #endif
323 #endif
324 
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.
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...
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...