Intrepid
Intrepid_CellTools_Kokkos.hpp
1 #ifndef INTREPID_CELTOOLS_KOKKOS_HPP
2 #define INTREPID_CELTOOLS_KOKKOS_HPP
3 
4 #ifdef INTREPID_OLD_KOKKOS_CODE
5 
6 #include "Intrepid_CellTools.hpp"
7 /*
8 namespace Intrepid{
9 
10  template<class Scalar>
11  template<class Scalar1,class Scalar2,class Scalar3,class Layout,class MemorySpace>
12  void CellTools<Scalar>::setJacobianTemp(Kokkos::View<Scalar1,Layout,MemorySpace> & jacobian,
13  const Kokkos::View<Scalar2,Layout,MemorySpace> & points,
14  const Kokkos::View<Scalar3,Layout,MemorySpace> & cellWorkset,
15  const shards::CellTopology & cellTopo,
16  const int & whichCell){
17 
18  CellTools<Scalar>::setJacobianTempSpecKokkos<Kokkos::View<Scalar1,Layout,MemorySpace>,Kokkos::View<Scalar2,Layout,MemorySpace>, Kokkos::View<Scalar3,Layout,MemorySpace>, Rank<Kokkos::View<Scalar2,Layout,MemorySpace> >::value ,Rank<Kokkos::View<Scalar1,Layout,MemorySpace> >::value>(jacobian, points, cellWorkset, cellTopo, whichCell);
19  }
20 
21 template<class ViewType1, class ViewType2>
22 struct points_temppointsCopy2 {
23  ViewType1 tempPoints;
24  typename ViewType2::const_type points;
25 
26  points_temppointsCopy2(ViewType1 tempPoints_, ViewType2 points_):tempPoints(tempPoints_),points(points_) {}
27 
28  KOKKOS_INLINE_FUNCTION
29  void operator() (int pt) const {
30  for(int dm = 0; dm < points.dimension(1) ; dm++){
31  tempPoints(pt, dm) = points(pt, dm);
32  }
33  }
34 };
35 template<class ViewType1, class ViewType2,class ViewType3>
36 struct calculateJacobian2_4 {
37  ViewType1 jacobian;
38  typename ViewType2::const_type cellWorkset;
39  typename ViewType3::const_type basisGrads;
40  calculateJacobian2_4(ViewType1 jacobian_, ViewType2 cellWorkset_, ViewType3 basisGrads_):jacobian(jacobian_),cellWorkset(cellWorkset_),basisGrads(basisGrads_) {}
41 
42  KOKKOS_INLINE_FUNCTION
43  void operator() (int cellOrd, int numPoints, int spaceDim, int basisCardinality) const {
44  for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
45  for(int row = 0; row < spaceDim; row++){
46  for(int col = 0; col < spaceDim; col++){
47  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
48  jacobian(cellOrd, pointOrd, row, col) += cellWorkset(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
49  } // bfOrd
50  } // col
51  } // row
52  } // pointOrd
53  }
54 };
55  template<class Scalar>
56  template< class ArrayJac, class ArrayPoint, class ArrayCell>
57  struct CellTools<Scalar>::setJacobianTempSpecKokkos<ArrayJac,ArrayPoint, ArrayCell, 2 ,4>{
58  setJacobianTempSpecKokkos(ArrayJac & jacobian,
59  const ArrayPoint & points,
60  const ArrayCell & cellWorkset,
61  const shards::CellTopology & cellTopo,
62  const int & whichCell){
63  INTREPID_VALIDATE( validateArguments_setJacobian(jacobian, points, cellWorkset, whichCell, cellTopo) );
64 
65  int spaceDim = (int)cellTopo.getDimension();
66  int numCells = cellWorkset.dimension(0);
67  //points can be rank-2 (P,D), or rank-3 (C,P,D)
68  int numPoints = points.dimension(0);
69 
70  // Jacobian is computed using gradients of an appropriate H(grad) basis function: define RCP to the base class
71  Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
72 
73  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
74  switch( cellTopo.getKey() ){
75 
76  // Standard Base topologies (number of cellWorkset = number of vertices)
77  case shards::Line<2>::key:
78  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
79  break;
80 
81  case shards::Triangle<3>::key:
82  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
83  break;
84 
85  case shards::Quadrilateral<4>::key:
86  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
87  break;
88 
89  case shards::Tetrahedron<4>::key:
90  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
91  break;
92 
93  case shards::Hexahedron<8>::key:
94  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
95  break;
96 
97  case shards::Wedge<6>::key:
98  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
99  break;
100 
101  case shards::Pyramid<5>::key:
102  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
103  break;
104 
105  // Standard Extended topologies
106  case shards::Triangle<6>::key:
107  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
108  break;
109  case shards::Quadrilateral<9>::key:
110  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
111  break;
112 
113  case shards::Tetrahedron<10>::key:
114  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
115  break;
116 
117  case shards::Tetrahedron<11>::key:
118  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
119  break;
120 
121  case shards::Hexahedron<20>::key:
122  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
123  break;
124 
125  case shards::Hexahedron<27>::key:
126  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
127  break;
128 
129  case shards::Wedge<15>::key:
130  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
131  break;
132 
133  case shards::Wedge<18>::key:
134  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
135  break;
136 
137  case shards::Pyramid<13>::key:
138  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
139  break;
140 
141  // These extended topologies are not used for mapping purposes
142  case shards::Quadrilateral<8>::key:
143  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
144  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
145  break;
146 
147  // Base and Extended Line, Beam and Shell topologies
148  case shards::Line<3>::key:
149  case shards::Beam<2>::key:
150  case shards::Beam<3>::key:
151  case shards::ShellLine<2>::key:
152  case shards::ShellLine<3>::key:
153  case shards::ShellTriangle<3>::key:
154  case shards::ShellTriangle<6>::key:
155  case shards::ShellQuadrilateral<4>::key:
156  case shards::ShellQuadrilateral<8>::key:
157  case shards::ShellQuadrilateral<9>::key:
158  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
159  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
160  break;
161  default:
162  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
163  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported.");
164  }// switch
165 
166  // Temp (F,P,D) array for the values of basis functions gradients at the reference points
167  int basisCardinality = HGRAD_Basis -> getCardinality();
168  // FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
169  Kokkos::View<Scalar***> basisGrads(basisCardinality, numPoints, spaceDim);
170  // Initialize jacobian
171 
172  Kokkos::deep_copy( jacobian, 0.0);
173 
174  // getValues requires rank-2 (P,D) input array, but points cannot be passed directly as argument because they are a user type
175  Kokkos::View<Scalar**>tempPoints(points.dimension(0), points.dimension(1));
176  // FieldContainer<Scalar> tempPoints( points.dimension(0), points.dimension(1) );
177  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
178 
179  parallel_for(points.dimension(0),points_temppointsCopy2<Kokkos::View<Scalar**>,ArrayPoint>(tempPoints,points));
180  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
181 
182  // The outer loops select the multi-index of the Jacobian entry: cell, point, row, col
183  // If whichCell = -1, all jacobians are computed, otherwise a single cell jacobian is computed
184  int cellLoop = numCells ;
185 
186  parallel_for(cellLoop,calculateJacobian2_4<ArrayJac,ArrayCell, Kokkos::View<Scalar***> >(jacobian, cellWorkset, basisGrads),numPoints,spaceDim,basisCardinality);
187  }
188  };
189 
190 template<class ViewType1, class ViewType2,class ViewType3>
191 struct calculateJacobian2_3 {
192  ViewType1 jacobian;
193  typename ViewType2::const_type cellWorkset;
194  typename ViewType3::const_type basisGrads;
195  calculateJacobian2_3(ViewType1 jacobian_, ViewType2 cellWorkset_, ViewType3 basisGrads_):jacobian(jacobian_),cellWorkset(cellWorkset_),basisGrads(basisGrads_) {}
196 
197  KOKKOS_INLINE_FUNCTION
198  void operator() (int pointOrd, int numPoints, int spaceDim, int basisCardinality,int whichCell) const {
199 
200  for(int row = 0; row < spaceDim; row++){
201  for(int col = 0; col < spaceDim; col++){
202  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
203  jacobian(pointOrd, row, col) += cellWorkset(whichCell, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
204  } // bfOrd
205  } // col
206  } // row
207 
208  }
209 };
210  template<class Scalar>
211  template< class ArrayJac, class ArrayPoint, class ArrayCell>
212  struct CellTools<Scalar>::setJacobianTempSpecKokkos<ArrayJac,ArrayPoint, ArrayCell, 2 ,3>{
213  setJacobianTempSpecKokkos(ArrayJac & jacobian,
214  const ArrayPoint & points,
215  const ArrayCell & cellWorkset,
216  const shards::CellTopology & cellTopo,
217  const int & whichCell){
218  INTREPID_VALIDATE( validateArguments_setJacobian(jacobian, points, cellWorkset, whichCell, cellTopo) );
219 
220  int spaceDim = (int)cellTopo.getDimension();
221  int numCells = cellWorkset.dimension(0);
222  //points can be rank-2 (P,D), or rank-3 (C,P,D)
223  int numPoints = points.dimension(0);
224 
225  // Jacobian is computed using gradients of an appropriate H(grad) basis function: define RCP to the base class
226  Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
227 
228  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
229  switch( cellTopo.getKey() ){
230 
231  // Standard Base topologies (number of cellWorkset = number of vertices)
232  case shards::Line<2>::key:
233  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
234  break;
235 
236  case shards::Triangle<3>::key:
237  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
238  break;
239 
240  case shards::Quadrilateral<4>::key:
241  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
242  break;
243 
244  case shards::Tetrahedron<4>::key:
245  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
246  break;
247 
248  case shards::Hexahedron<8>::key:
249  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
250  break;
251 
252  case shards::Wedge<6>::key:
253  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
254  break;
255 
256  case shards::Pyramid<5>::key:
257  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
258  break;
259 
260  // Standard Extended topologies
261  case shards::Triangle<6>::key:
262  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
263  break;
264  case shards::Quadrilateral<9>::key:
265  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
266  break;
267 
268  case shards::Tetrahedron<10>::key:
269  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
270  break;
271 
272  case shards::Tetrahedron<11>::key:
273  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
274  break;
275 
276  case shards::Hexahedron<20>::key:
277  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
278  break;
279 
280  case shards::Hexahedron<27>::key:
281  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
282  break;
283 
284  case shards::Wedge<15>::key:
285  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
286  break;
287 
288  case shards::Wedge<18>::key:
289  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
290  break;
291 
292  case shards::Pyramid<13>::key:
293  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
294  break;
295 
296  // These extended topologies are not used for mapping purposes
297  case shards::Quadrilateral<8>::key:
298  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
299  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
300  break;
301 
302  // Base and Extended Line, Beam and Shell topologies
303  case shards::Line<3>::key:
304  case shards::Beam<2>::key:
305  case shards::Beam<3>::key:
306  case shards::ShellLine<2>::key:
307  case shards::ShellLine<3>::key:
308  case shards::ShellTriangle<3>::key:
309  case shards::ShellTriangle<6>::key:
310  case shards::ShellQuadrilateral<4>::key:
311  case shards::ShellQuadrilateral<8>::key:
312  case shards::ShellQuadrilateral<9>::key:
313  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
314  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
315  break;
316  default:
317  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
318  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported.");
319  }// switch
320 
321  // Temp (F,P,D) array for the values of basis functions gradients at the reference points
322  int basisCardinality = HGRAD_Basis -> getCardinality();
323  // FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
324  Kokkos::View<Scalar***> basisGrads(basisCardinality, numPoints, spaceDim);
325  // Initialize jacobian
326 
327  Kokkos::deep_copy( jacobian, 0.0);
328 
329  // getValues requires rank-2 (P,D) input array, but points cannot be passed directly as argument because they are a user type
330  Kokkos::View<Scalar**>tempPoints(points.dimension(0), points.dimension(1));
331  // FieldContainer<Scalar> tempPoints( points.dimension(0), points.dimension(1) );
332  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
333 
334  parallel_for(points.dimension(0),points_temppointsCopy2<Kokkos::View<Scalar**>,ArrayPoint>(tempPoints,points));
335  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
336 
337  // The outer loops select the multi-index of the Jacobian entry: point, row, col
338  parallel_for(numPoints,calculateJacobian2_3<ArrayJac,ArrayCell, Kokkos::View<Scalar***> >(jacobian, cellWorkset, basisGrads),numPoints,spaceDim,basisCardinality,whichCell);
339  }
340  };
341 
342 template<class ViewType1, class ViewType2,class ViewType3>
343 struct calculateJacobian3_4 {
344  ViewType1 jacobian;
345  typename ViewType2::const_type cellWorkset;
346  typename ViewType3::const_type basisGrads;
347  calculateJacobian3_4(ViewType1 jacobian_, ViewType2 cellWorkset_, ViewType3 basisGrads_):jacobian(jacobian_),cellWorkset(cellWorkset_),basisGrads(basisGrads_) {}
348 
349  KOKKOS_INLINE_FUNCTION
350  void operator() (int pointOrd, int numPoints, int spaceDim, int basisCardinality,int whichCell) const {
351 
352  for(int row = 0; row < spaceDim; row++){
353  for(int col = 0; col < spaceDim; col++){
354  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
355  jacobian(pointOrd, row, col) += cellWorkset(whichCell, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
356  } // bfOrd
357  } // col
358  } // row
359 
360  }
361 };
362  template<class Scalar>
363  template<class ArrayJac, class ArrayPoint, class ArrayCell>
364  struct CellTools<Scalar>::setJacobianTempSpecKokkos<ArrayJac,ArrayPoint, ArrayCell, 3 , 4>{
365  setJacobianTempSpecKokkos(ArrayJac & jacobian,
366  const ArrayPoint & points,
367  const ArrayCell & cellWorkset,
368  const shards::CellTopology & cellTopo,
369  const int & whichCell){
370  INTREPID_VALIDATE( validateArguments_setJacobian(jacobian, points, cellWorkset, whichCell, cellTopo) );
371 
372  int spaceDim = (int)cellTopo.getDimension();
373  int numCells = cellWorkset.dimension(0);
374 
375  int numPoints = points.dimension(1);
376 
377  // Jacobian is computed using gradients of an appropriate H(grad) basis function: define RCP to the base class
378  Teuchos::RCP< Basis< Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
379 
380  // Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
381  switch( cellTopo.getKey() ){
382 
383  // Standard Base topologies (number of cellWorkset = number of vertices)
384  case shards::Line<2>::key:
385  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
386  break;
387 
388  case shards::Triangle<3>::key:
389  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
390  break;
391 
392  case shards::Quadrilateral<4>::key:
393  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
394  break;
395 
396  case shards::Tetrahedron<4>::key:
397  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
398  break;
399 
400  case shards::Hexahedron<8>::key:
401  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
402  break;
403 
404  case shards::Wedge<6>::key:
405  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
406  break;
407 
408  case shards::Pyramid<5>::key:
409  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
410  break;
411 
412  // Standard Extended topologies
413  case shards::Triangle<6>::key:
414  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
415  break;
416  case shards::Quadrilateral<9>::key:
417  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
418  break;
419 
420  case shards::Tetrahedron<10>::key:
421  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
422  break;
423 
424  case shards::Tetrahedron<11>::key:
425  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
426  break;
427 
428  case shards::Hexahedron<20>::key:
429  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_I2_FEM<Scalar, FieldContainer<Scalar> >() );
430  break;
431 
432  case shards::Hexahedron<27>::key:
433  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
434  break;
435 
436  case shards::Wedge<15>::key:
437  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_I2_FEM<Scalar, FieldContainer<Scalar> >() );
438  break;
439 
440  case shards::Wedge<18>::key:
441  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
442  break;
443 
444  case shards::Pyramid<13>::key:
445  HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_I2_FEM<Scalar, FieldContainer<Scalar> >() );
446  break;
447 
448  // These extended topologies are not used for mapping purposes
449  case shards::Quadrilateral<8>::key:
450  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
451  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
452  break;
453 
454  // Base and Extended Line, Beam and Shell topologies
455  case shards::Line<3>::key:
456  case shards::Beam<2>::key:
457  case shards::Beam<3>::key:
458  case shards::ShellLine<2>::key:
459  case shards::ShellLine<3>::key:
460  case shards::ShellTriangle<3>::key:
461  case shards::ShellTriangle<6>::key:
462  case shards::ShellQuadrilateral<4>::key:
463  case shards::ShellQuadrilateral<8>::key:
464  case shards::ShellQuadrilateral<9>::key:
465  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
466  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported. ");
467  break;
468  default:
469  TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
470  ">>> ERROR (Intrepid::CellTools::setJacobian): Cell topology not supported.");
471  }// switch
472 
473  // Temp (F,P,D) array for the values of basis functions gradients at the reference points
474  int basisCardinality = HGRAD_Basis -> getCardinality();
475  FieldContainer<Scalar> basisGrads(basisCardinality, numPoints, spaceDim);
476 
477  // Initialize jacobian
478  for(int i = 0; i < jacobian.size(); i++){
479  jacobian[i] = 0.0;
480  }
481 
482  // getValues requires rank-2 (P,D) input array, refPoints cannot be used as argument: need temp (P,D) array
483  FieldContainer<Scalar> tempPoints( points.dimension(1), points.dimension(2) );
484 
485  for(int cellOrd = 0; cellOrd < numCells; cellOrd++) {
486 
487  // Copy point set corresponding to this cell oridinal to the temp (P,D) array
488  for(int pt = 0; pt < points.dimension(1); pt++){
489  for(int dm = 0; dm < points.dimension(2) ; dm++){
490  tempPoints(pt, dm) = points(cellOrd, pt, dm);
491  }//dm
492  }//pt
493 
494  // Compute gradients of basis functions at this set of ref. points
495  HGRAD_Basis -> getValues(basisGrads, tempPoints, OPERATOR_GRAD);
496 
497  // Compute jacobians for the point set corresponding to the current cellordinal
498  for(int pointOrd = 0; pointOrd < numPoints; pointOrd++) {
499  for(int row = 0; row < spaceDim; row++){
500  for(int col = 0; col < spaceDim; col++){
501 
502  // The entry is computed by contracting the basis index. Number of basis functions and vertices must be the same
503  for(int bfOrd = 0; bfOrd < basisCardinality; bfOrd++){
504  jacobian(cellOrd, pointOrd, row, col) += cellWorkset(cellOrd, bfOrd, row)*basisGrads(bfOrd, pointOrd, col);
505  } // bfOrd
506  } // col
507  } // row
508  } // pointOrd
509  }//cellOrd
510  }//end function
511  };
512 }//end namespace Intrepid_Kokkos
513 */
514 #endif
515 
516 #endif
Intrepid_CellTools.hpp
Header file for the Intrepid::CellTools class.
Generated by   doxygen 1.8.5