Intrepid2
Intrepid2_HGRAD_TRI_Cn_FEM_ORTHDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
50 #ifndef __INTREPID2_HGRAD_TRI_CN_FEM_ORTH_DEF_HPP__
51 #define __INTREPID2_HGRAD_TRI_CN_FEM_ORTH_DEF_HPP__
52 
53 namespace Intrepid2 {
54 // -------------------------------------------------------------------------------------
55 
56 namespace Impl {
57 
58 template<typename OutputViewType,
59 typename inputViewType,
60 typename workViewType,
61 bool hasDeriv>
62 KOKKOS_INLINE_FUNCTION
63 void OrthPolynomialTri<OutputViewType,inputViewType,workViewType,hasDeriv,0>::generate(
64  OutputViewType output,
65  const inputViewType input,
66  workViewType /*work*/,
67  const ordinal_type order ) {
68 
69  constexpr ordinal_type spaceDim = 2;
70  constexpr ordinal_type maxNumPts = Parameters::MaxNumPtsPerBasisEval;
71 
72  typedef typename OutputViewType::value_type value_type;
73 
74  auto output0 = (hasDeriv) ? Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL(),0) : Kokkos::subview(output, Kokkos::ALL(), Kokkos::ALL());
75 
76  const ordinal_type
77  npts = input.extent(0);
78 
79  const auto z = input;
80 
81  // each point needs to be transformed from Pavel's element
82  // z(i,0) --> (2.0 * z(i,0) - 1.0)
83  // z(i,1) --> (2.0 * z(i,1) - 1.0)
84 
85 
86  // set D^{0,0} = 1.0
87  {
88  const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(0,0);
89  for (ordinal_type i=0;i<npts;++i) {
90  output0(loc, i) = 1.0;
91  if(hasDeriv) {
92  output.access(loc,i,1) = 0;
93  output.access(loc,i,2) = 0;
94  }
95  }
96  }
97 
98  if (order > 0) {
99  value_type f1[maxNumPts]={},f2[maxNumPts]={}, df2_1[maxNumPts]={};
100  value_type df1_0, df1_1;
101 
102  for (ordinal_type i=0;i<npts;++i) {
103  f1[i] = 0.5 * (1.0+2.0*(2.0*z(i,0)-1.0)+(2.0*z(i,1)-1.0)); // \eta_1 * (1 - \eta_2)/2
104  f2[i] = std::pow(z(i,1)-1,2); //( (1 - \eta_2)/2 )^2
105  if(hasDeriv) {
106  df1_0 = 2.0;
107  df1_1 = 1.0;
108  df2_1[i] = 2.0*(z(i,1)-1);
109  }
110  }
111 
112  // set D^{1,0} = f1
113  {
114  const ordinal_type loc = Intrepid2::getPnEnumeration<spaceDim>(1,0);
115  for (ordinal_type i=0;i<npts;++i) {
116  output0(loc, i) = f1[i];
117  if(hasDeriv) {
118  output.access(loc,i,1) = df1_0;
119  output.access(loc,i,2) = df1_1;
120  }
121  }
122  }
123 
124  // recurrence in p
125  for (ordinal_type p=1;p<order;p++) {
126  const ordinal_type
127  loc = Intrepid2::getPnEnumeration<spaceDim>(p,0),
128  loc_p1 = Intrepid2::getPnEnumeration<spaceDim>(p+1,0),
129  loc_m1 = Intrepid2::getPnEnumeration<spaceDim>(p-1,0);
130 
131  const value_type
132  a = (2.0*p+1.0)/(1.0+p),
133  b = p / (p+1.0);
134 
135  for (ordinal_type i=0;i<npts;++i) {
136  output0(loc_p1,i) = ( a * f1[i] * output0(loc,i) -
137  b * f2[i] * output0(loc_m1,i) );
138  if(hasDeriv) {
139  output.access(loc_p1,i,1) = a * (f1[i] * output.access(loc,i,1) + df1_0 * output0(loc,i)) -
140  b * f2[i] * output.access(loc_m1,i,1) ;
141  output.access(loc_p1,i,2) = a * (f1[i] * output.access(loc,i,2) + df1_1 * output0(loc,i)) -
142  b * (df2_1[i] * output0(loc_m1,i) + f2[i] * output.access(loc_m1,i,2)) ;
143  }
144  }
145  }
146 
147  // D^{p,1}
148  for (ordinal_type p=0;p<order;++p) {
149  const ordinal_type
150  loc_p_0 = Intrepid2::getPnEnumeration<spaceDim>(p,0),
151  loc_p_1 = Intrepid2::getPnEnumeration<spaceDim>(p,1);
152 
153  for (ordinal_type i=0;i<npts;++i) {
154  output0(loc_p_1,i) = output0(loc_p_0,i)*0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
155  if(hasDeriv) {
156  output.access(loc_p_1,i,1) = output.access(loc_p_0,i,1)*0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0));
157  output.access(loc_p_1,i,2) = output.access(loc_p_0,i,2)*0.5*(1.0+2.0*p+(3.0+2.0*p)*(2.0*z(i,1)-1.0)) + output0(loc_p_0,i)*(3.0+2.0*p);
158  }
159  }
160  }
161 
162 
163  // recurrence in q
164  for (ordinal_type p=0;p<order-1;++p)
165  for (ordinal_type q=1;q<order-p;++q) {
166  const ordinal_type
167  loc_p_qp1 = Intrepid2::getPnEnumeration<spaceDim>(p,q+1),
168  loc_p_q = Intrepid2::getPnEnumeration<spaceDim>(p,q),
169  loc_p_qm1 = Intrepid2::getPnEnumeration<spaceDim>(p,q-1);
170 
171  value_type a,b,c;
172  Intrepid2::getJacobyRecurrenceCoeffs(a,b,c, 2*p+1,0,q);
173  for (ordinal_type i=0;i<npts;++i) {
174  output0(loc_p_qp1,i) = (a*(2.0*z(i,1)-1.0)+b)*output0(loc_p_q,i)
175  - c*output0(loc_p_qm1,i) ;
176  if(hasDeriv) {
177  output.access(loc_p_qp1,i,1) = (a*(2.0*z(i,1)-1.0)+b)*output.access(loc_p_q,i,1)
178  - c*output.access(loc_p_qm1,i,1) ;
179  output.access(loc_p_qp1,i,2) = (a*(2.0*z(i,1)-1.0)+b)*output.access(loc_p_q,i,2) +2*a*output0(loc_p_q,i)
180  - c*output.access(loc_p_qm1,i,2) ;
181  }
182  }
183  }
184  }
185 
186  // orthogonalize
187  for (ordinal_type p=0;p<=order;++p)
188  for (ordinal_type q=0;q<=order-p;++q)
189  for (ordinal_type i=0;i<npts;++i) {
190  output0(Intrepid2::getPnEnumeration<spaceDim>(p,q),i) *= std::sqrt( (p+0.5)*(p+q+1.0));
191  if(hasDeriv) {
192  output.access(Intrepid2::getPnEnumeration<spaceDim>(p,q),i,1) *= std::sqrt( (p+0.5)*(p+q+1.0));
193  output.access(Intrepid2::getPnEnumeration<spaceDim>(p,q),i,2) *= std::sqrt( (p+0.5)*(p+q+1.0));
194  }
195  }
196 }
197 
198 template<typename OutputViewType,
199 typename inputViewType,
200 typename workViewType,
201 bool hasDeriv>
202 KOKKOS_INLINE_FUNCTION
203 void OrthPolynomialTri<OutputViewType,inputViewType,workViewType,hasDeriv,1>::generate(
204  OutputViewType output,
205  const inputViewType input,
206  workViewType work,
207  const ordinal_type order ) {
208  constexpr ordinal_type spaceDim = 2;
209  const ordinal_type
210  npts = input.extent(0),
211  card = output.extent(0);
212 
213  workViewType dummyView;
214  OrthPolynomialTri<workViewType,inputViewType,workViewType,hasDeriv,0>::generate(work, input, dummyView, order);
215  for (ordinal_type i=0;i<card;++i)
216  for (ordinal_type j=0;j<npts;++j)
217  for (ordinal_type k=0;k<spaceDim;++k)
218  output.access(i,j,k) = work(i,j,k+1);
219 }
220 
221 
222 // when n >= 2, use recursion
223 template<typename OutputViewType,
224 typename inputViewType,
225 typename workViewType,
226 bool hasDeriv,
227 ordinal_type n>
228 KOKKOS_INLINE_FUNCTION
229 void OrthPolynomialTri<OutputViewType,inputViewType,workViewType,hasDeriv,n>::generate(
230  OutputViewType /* output */,
231  const inputViewType /* input */,
232  workViewType /* work */,
233  const ordinal_type /* order */ ) {
234 #if 0 //#ifdef HAVE_INTREPID2_SACADO
235 
236 constexpr ordinal_type spaceDim = 2;
237 constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
238 
239 typedef typename OutputViewType::value_type value_type;
240 typedef Sacado::Fad::SFad<value_type,spaceDim> fad_type;
241 
242 const ordinal_type
243 npts = input.extent(0),
244 card = output.extent(0);
245 
246 // use stack buffer
247 fad_type inBuf[Parameters::MaxNumPtsPerBasisEval][spaceDim],
248 outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][n];
249 
250 typedef typename inputViewType::memory_space memory_space;
251 typedef typename Kokkos::View<fad_type***, memory_space> outViewType;
252 typedef typename Kokkos::View<fad_type**, memory_space> inViewType;
253 auto vcprop = Kokkos::common_view_alloc_prop(input);
254 
255 inViewType in(Kokkos::view_wrap((value_type*)&inBuf[0][0], vcprop), npts, spaceDim);
256 outViewType out(Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, n);
257 
258 for (ordinal_type i=0;i<npts;++i)
259  for (ordinal_type j=0;j<spaceDim;++j) {
260  in(i,j) = input(i,j);
261  in(i,j).diff(j,spaceDim);
262  }
263 
264 typedef typename Kokkos::DynRankView<fad_type, memory_space> outViewType_;
265 outViewType_ workView;
266 if (n==2) {
267  //char outBuf[bufSize*sizeof(typename inViewType::value_type)];
268  fad_type outBuf[maxCard][Parameters::MaxNumPtsPerBasisEval][spaceDim+1];
269  auto vcprop = Kokkos::common_view_alloc_prop(in);
270  workView = outViewType_( Kokkos::view_wrap((value_type*)&outBuf[0][0][0], vcprop), card, npts, spaceDim+1);
271 }
272 OrthPolynomialTri<outViewType,inViewType,outViewType_,hasDeriv,n-1>::generate(out, in, workView, order);
273 
274 for (ordinal_type i=0;i<card;++i)
275  for (ordinal_type j=0;j<npts;++j) {
276  for (ordinal_type i_dx = 0; i_dx <= n; ++i_dx) {
277  ordinal_type i_dy = n-i_dx;
278  ordinal_type i_Dn = i_dy;
279  if(i_dx > 0) {
280  //n=2: (f_x)_x, (f_y)_x
281  //n=3: (f_xx)_x, (f_xy)_x, (f_yy)_x
282  ordinal_type i_Dnm1 = i_dy;
283  output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(0);
284  }
285  else {
286  //n=2: (f_y)_y, (f_z)_y
287  //n=3: (f_yy)_y
288  ordinal_type i_Dnm1 = i_dy-1;
289  output.access(i,j,i_Dn) = out(i,j,i_Dnm1).dx(1);
290  }
291  }
292  }
293 #else
294 INTREPID2_TEST_FOR_ABORT( true,
295  ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::OrthPolynomialTri) Computing of second and higher-order derivatives is not currently supported");
296 #endif
297 }
298 
299 
300 
301 template<EOperator opType>
302 template<typename OutputViewType,
303 typename inputViewType,
304 typename workViewType>
305 KOKKOS_INLINE_FUNCTION
306 void
307 Basis_HGRAD_TRI_Cn_FEM_ORTH::Serial<opType>::
308 getValues( OutputViewType output,
309  const inputViewType input,
310  workViewType work,
311  const ordinal_type order) {
312  switch (opType) {
313  case OPERATOR_VALUE: {
314  OrthPolynomialTri<OutputViewType,inputViewType,workViewType,false,0>::generate( output, input, work, order );
315  break;
316  }
317  case OPERATOR_GRAD:
318  case OPERATOR_D1: {
319  OrthPolynomialTri<OutputViewType,inputViewType,workViewType,true,1>::generate( output, input, work, order );
320  break;
321  }
322  case OPERATOR_D2: {
323  OrthPolynomialTri<OutputViewType,inputViewType,workViewType,true,2>::generate( output, input, work, order );
324  break;
325  }
326  default: {
327  INTREPID2_TEST_FOR_ABORT( true,
328  ">>> ERROR: (Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::Serial::getValues) operator is not supported");
329  }
330  }
331 }
332 
333 template<typename SpT, ordinal_type numPtsPerEval,
334 typename outputValueValueType, class ...outputValueProperties,
335 typename inputPointValueType, class ...inputPointProperties>
336 void
337 Basis_HGRAD_TRI_Cn_FEM_ORTH::
338 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
339  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
340  const ordinal_type order,
341  const EOperator operatorType ) {
342  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
343  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
344  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
345 
346  // loopSize corresponds to the # of points
347  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
348  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
349  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
350  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
351 
352  typedef typename inputPointViewType::value_type inputPointType;
353  const ordinal_type cardinality = outputValues.extent(0);
354  const ordinal_type spaceDim = 2;
355 
356  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
357  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
358 
359  switch (operatorType) {
360  case OPERATOR_VALUE: {
361  workViewType dummyWorkView;
362  typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_VALUE,numPtsPerEval> FunctorType;
363  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, dummyWorkView, order) );
364  break;
365  }
366  case OPERATOR_GRAD:
367  case OPERATOR_D1: {
368  workViewType work(Kokkos::view_alloc("Basis_HGRAD_TRI_In_FEM_ORTH::getValues::work", vcprop), cardinality, inputPoints.extent(0), spaceDim+1);
369  typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D1,numPtsPerEval> FunctorType;
370  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, work, order) );
371  break;
372  }
373  case OPERATOR_D2:{
374  workViewType dummyWorkView;
375  typedef Functor<outputValueViewType,inputPointViewType,workViewType,OPERATOR_D2,numPtsPerEval> FunctorType;
376  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints ,dummyWorkView, order) );
377  break;
378  }
379  default: {
380  INTREPID2_TEST_FOR_EXCEPTION( !Intrepid2::isValidOperator(operatorType), std::invalid_argument,
381  ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM_ORTH): invalid operator type");
382  }
383  }
384 }
385 }
386 
387 
388 // -------------------------------------------------------------------------------------
389 
390 template<typename SpT, typename OT, typename PT>
391 Basis_HGRAD_TRI_Cn_FEM_ORTH<SpT,OT,PT>::
392 Basis_HGRAD_TRI_Cn_FEM_ORTH( const ordinal_type order ) {
393 
394  constexpr ordinal_type spaceDim = 2;
395  this->basisCardinality_ = Intrepid2::getPnCardinality<spaceDim>(order);
396  this->basisDegree_ = order;
397  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
398  this->basisType_ = BASIS_FEM_HIERARCHICAL;
399  this->basisCoordinates_ = COORDINATES_CARTESIAN;
400  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
401 
402  // initialize tags
403  {
404  // Basis-dependent initializations
405  constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
406  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
407  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
408  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
409 
410  constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
411  ordinal_type tags[maxCard][tagSize];
412  const ordinal_type card = this->basisCardinality_;
413  for (ordinal_type i=0;i<card;++i) {
414  tags[i][0] = 2; // these are all "internal" i.e. "volume" DoFs
415  tags[i][1] = 0; // there is only one line
416  tags[i][2] = i; // local DoF id
417  tags[i][3] = card; // total number of DoFs
418  }
419 
420  OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
421 
422  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
423  // tags are constructed on host
424  this->setOrdinalTagData(this->tagToOrdinal_,
425  this->ordinalToTag_,
426  tagView,
427  this->basisCardinality_,
428  tagSize,
429  posScDim,
430  posScOrd,
431  posDfOrd);
432  }
433 
434  // dof coords is not applicable to hierarchical functions
435 }
436 
437 }
438 #endif
Intrepid2::Basis::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecSpaceType::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:135
Intrepid2::Parameters::MaxNumPtsPerBasisEval
static constexpr ordinal_type MaxNumPtsPerBasisEval
The maximum number of points to eval in serial mode.
Definition: Intrepid2_Types.hpp:121
Intrepid2::Basis_HGRAD_TRI_Cn_FEM_ORTH::Basis_HGRAD_TRI_Cn_FEM_ORTH
Basis_HGRAD_TRI_Cn_FEM_ORTH(const ordinal_type order)
Constructor.
Definition: Intrepid2_HGRAD_TRI_Cn_FEM_ORTHDef.hpp:392
Generated by   doxygen 1.8.5