Intrepid2
Intrepid2_CellToolsDefJacobian.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 
43 
49 #ifndef __INTREPID2_CELLTOOLS_DEF_HPP__
50 #define __INTREPID2_CELLTOOLS_DEF_HPP__
51 
52 // disable clang warnings
53 #if defined (__clang__) && !defined (__INTEL_COMPILER)
54 #pragma clang system_header
55 #endif
56 
57 #include "Intrepid2_Kernels.hpp"
58 
59 namespace Intrepid2 {
60 
61 
62  //============================================================================================//
63  // //
64  // Jacobian, inverse Jacobian and Jacobian determinant //
65  // //
66  //============================================================================================//
67 
68  namespace FunctorCellTools {
72  template<typename jacobianViewType,
73  typename worksetCellType,
74  typename basisGradType>
75  struct F_setJacobian {
76  jacobianViewType _jacobian;
77  const worksetCellType _worksetCells;
78  const basisGradType _basisGrads;
79  const int _startCell;
80  const int _endCell;
81 
82  KOKKOS_INLINE_FUNCTION
83  F_setJacobian( jacobianViewType jacobian_,
84  worksetCellType worksetCells_,
85  basisGradType basisGrads_,
86  const int startCell_,
87  const int endCell_)
88  : _jacobian(jacobian_), _worksetCells(worksetCells_), _basisGrads(basisGrads_),
89  _startCell(startCell_), _endCell(endCell_) {}
90 
91  KOKKOS_INLINE_FUNCTION
92  void operator()(const ordinal_type cell,
93  const ordinal_type point) const {
94 
95  const ordinal_type dim = _jacobian.extent(2); // dim2 and dim3 should match
96 
97  const ordinal_type gradRank = _basisGrads.rank();
98  if ( gradRank == 3)
99  {
100  const ordinal_type cardinality = _basisGrads.extent(0);
101  for (ordinal_type i=0;i<dim;++i)
102  for (ordinal_type j=0;j<dim;++j) {
103  _jacobian(cell, point, i, j) = 0;
104  for (ordinal_type bf=0;bf<cardinality;++bf)
105  _jacobian(cell, point, i, j) += _worksetCells(cell+_startCell, bf, i) * _basisGrads(bf, point, j);
106  }
107  }
108  else
109  {
110  const ordinal_type cardinality = _basisGrads.extent(1);
111  for (ordinal_type i=0;i<dim;++i)
112  for (ordinal_type j=0;j<dim;++j) {
113  _jacobian(cell, point, i, j) = 0;
114  for (ordinal_type bf=0;bf<cardinality;++bf)
115  _jacobian(cell, point, i, j) += _worksetCells(cell+_startCell, bf, i) * _basisGrads(cell, bf, point, j);
116  }
117  }
118  }
119  };
120  }
121 
122  template<typename DeviceType>
123  template<class PointScalar>
125  {
126  auto extents = jacobian.getExtents(); // C,P,D,D, which we reduce to C,P
127  auto variationTypes = jacobian.getVariationTypes();
128  const bool cellVaries = (variationTypes[0] != CONSTANT);
129  const bool pointVaries = (variationTypes[1] != CONSTANT);
130 
131  extents[2] = 1;
132  extents[3] = 1;
133  variationTypes[2] = CONSTANT;
134  variationTypes[3] = CONSTANT;
135 
136  if ( cellVaries && pointVaries )
137  {
138  auto data = jacobian.getUnderlyingView4();
139  auto detData = getMatchingViewWithLabel(data, "Jacobian det data", data.extent_int(0), data.extent_int(1));
140  return Data<PointScalar,DeviceType>(detData,2,extents,variationTypes);
141  }
142  else if (cellVaries || pointVaries)
143  {
144  auto data = jacobian.getUnderlyingView3();
145  auto detData = getMatchingViewWithLabel(data, "Jacobian det data", data.extent_int(0));
146  return Data<PointScalar,DeviceType>(detData,2,extents,variationTypes);
147  }
148  else
149  {
150  auto data = jacobian.getUnderlyingView1();
151  auto detData = getMatchingViewWithLabel(data, "Jacobian det data", 1);
152  return Data<PointScalar,DeviceType>(detData,2,extents,variationTypes);
153  }
154  }
155 
156  template<typename DeviceType>
157  template<class PointScalar>
159  {
160  auto extents = jacobian.getExtents(); // C,P,D,D
161  auto variationTypes = jacobian.getVariationTypes();
162  int jacDataRank = jacobian.getUnderlyingViewRank();
163 
164  if ( jacDataRank == 4 )
165  {
166  auto jacData = jacobian.getUnderlyingView4();
167  auto invData = getMatchingViewWithLabel(jacData, "Jacobian inv data",jacData.extent(0),jacData.extent(1),jacData.extent(2),jacData.extent(3));
168  return Data<PointScalar,DeviceType>(invData,4,extents,variationTypes);
169  }
170  else if (jacDataRank == 3)
171  {
172  auto jacData = jacobian.getUnderlyingView3();
173  auto invData = getMatchingViewWithLabel(jacData, "Jacobian inv data",jacData.extent(0),jacData.extent(1),jacData.extent(2));
174  return Data<PointScalar,DeviceType>(invData,4,extents,variationTypes);
175  }
176  else if (jacDataRank == 2)
177  {
178  auto jacData = jacobian.getUnderlyingView2();
179  auto invData = getMatchingViewWithLabel(jacData, "Jacobian inv data",jacData.extent(0),jacData.extent(1));
180  return Data<PointScalar,DeviceType>(invData,4,extents,variationTypes);
181  }
182  else if (jacDataRank == 1)
183  {
184  auto jacData = jacobian.getUnderlyingView1();
185  auto invData = getMatchingViewWithLabel(jacData, "Jacobian inv data",jacData.extent(0));
186  return Data<PointScalar,DeviceType>(invData,4,extents,variationTypes);
187  }
188  else
189  {
190  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "allocateJacobianInv requires jacobian to vary in *some* dimension…");
191  return Data<PointScalar,DeviceType>(); // unreachable statement; this line added to avoid compiler warning on CUDA
192  }
193  }
194 
195  template<typename DeviceType>
196  template<class PointScalar>
198  {
199  auto variationTypes = jacobian.getVariationTypes();
200 
201  const int CELL_DIM = 0;
202  const int POINT_DIM = 1;
203  const int D1_DIM = 2;
204  const bool cellVaries = (variationTypes[CELL_DIM] != CONSTANT);
205  const bool pointVaries = (variationTypes[POINT_DIM] != CONSTANT);
206 
207  auto det2 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c, const PointScalar &d) -> PointScalar
208  {
209  return a * d - b * c;
210  };
211 
212  auto det3 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c,
213  const PointScalar &d, const PointScalar &e, const PointScalar &f,
214  const PointScalar &g, const PointScalar &h, const PointScalar &i) -> PointScalar
215  {
216  return a * det2(e,f,h,i) - b * det2(d,f,g,i) + c * det2(d,e,g,h);
217  };
218 
219  auto det4 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c, const PointScalar &d,
220  const PointScalar &e, const PointScalar &f, const PointScalar &g, const PointScalar &h,
221  const PointScalar &i, const PointScalar &j, const PointScalar &k, const PointScalar &l,
222  const PointScalar &m, const PointScalar &n, const PointScalar &o, const PointScalar &p) -> PointScalar
223  {
224  return a * det3(f,g,h,j,k,l,n,o,p) - b * det3(e,g,h,i,k,l,m,o,p) + c * det3(e,f,h,i,j,l,m,n,p) - d * det3(e,f,g,i,j,k,m,n,o);
225  };
226 
227  if (variationTypes[D1_DIM] == BLOCK_PLUS_DIAGONAL)
228  {
229  if (cellVaries && pointVaries)
230  {
231  auto data = jacobian.getUnderlyingView3();
232  auto detData = jacobianDet.getUnderlyingView2();
233 
234  Kokkos::parallel_for(
235  Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>>({0,0},{data.extent_int(0),data.extent_int(1)}),
236  KOKKOS_LAMBDA (int cellOrdinal, int pointOrdinal) {
237  const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
238  const int numDiagonals = data.extent_int(2) - blockWidth * blockWidth;
239  const int spaceDim = blockWidth + numDiagonals;
240 
241  PointScalar det = 1.0;
242  switch (blockWidth)
243  {
244  case 0:
245  det = 1.0;
246  break;
247  case 1:
248  {
249  det = data(cellOrdinal,pointOrdinal,0);
250  break;
251  }
252  case 2:
253  {
254  const auto & a = data(cellOrdinal,pointOrdinal,0);
255  const auto & b = data(cellOrdinal,pointOrdinal,1);
256  const auto & c = data(cellOrdinal,pointOrdinal,2);
257  const auto & d = data(cellOrdinal,pointOrdinal,3);
258  det = det2(a,b,c,d);
259 
260  break;
261  }
262  case 3:
263  {
264  const auto & a = data(cellOrdinal,pointOrdinal,0);
265  const auto & b = data(cellOrdinal,pointOrdinal,1);
266  const auto & c = data(cellOrdinal,pointOrdinal,2);
267  const auto & d = data(cellOrdinal,pointOrdinal,3);
268  const auto & e = data(cellOrdinal,pointOrdinal,4);
269  const auto & f = data(cellOrdinal,pointOrdinal,5);
270  const auto & g = data(cellOrdinal,pointOrdinal,6);
271  const auto & h = data(cellOrdinal,pointOrdinal,7);
272  const auto & i = data(cellOrdinal,pointOrdinal,8);
273  det = det3(a,b,c,d,e,f,g,h,i);
274 
275  break;
276  }
277  case 4:
278  {
279  const auto & a = data(cellOrdinal,pointOrdinal,0);
280  const auto & b = data(cellOrdinal,pointOrdinal,1);
281  const auto & c = data(cellOrdinal,pointOrdinal,2);
282  const auto & d = data(cellOrdinal,pointOrdinal,3);
283  const auto & e = data(cellOrdinal,pointOrdinal,4);
284  const auto & f = data(cellOrdinal,pointOrdinal,5);
285  const auto & g = data(cellOrdinal,pointOrdinal,6);
286  const auto & h = data(cellOrdinal,pointOrdinal,7);
287  const auto & i = data(cellOrdinal,pointOrdinal,8);
288  const auto & j = data(cellOrdinal,pointOrdinal,9);
289  const auto & k = data(cellOrdinal,pointOrdinal,10);
290  const auto & l = data(cellOrdinal,pointOrdinal,11);
291  const auto & m = data(cellOrdinal,pointOrdinal,12);
292  const auto & n = data(cellOrdinal,pointOrdinal,13);
293  const auto & o = data(cellOrdinal,pointOrdinal,14);
294  const auto & p = data(cellOrdinal,pointOrdinal,15);
295  det = det4(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p);
296 
297  break;
298  }
299  default:
300  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 4 not supported for BLOCK_PLUS_DIAGONAL");
301  }
302  // incorporate the remaining, diagonal entries
303  const int offset = blockWidth * blockWidth;
304 
305  for (int d=blockWidth; d<spaceDim; d++)
306  {
307  const int index = d-blockWidth+offset;
308  det *= data(cellOrdinal,pointOrdinal,index);
309  }
310  detData(cellOrdinal,pointOrdinal) = det;
311  });
312  }
313  else if (cellVaries || pointVaries) // exactly one of cell,point vary -- whichever it is will be in rank 0 of data, invData
314  {
315  auto data = jacobian.getUnderlyingView2();
316  auto detData = jacobianDet.getUnderlyingView1();
317 
318  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,data.extent_int(0)),
319  KOKKOS_LAMBDA (const int &cellPointOrdinal) {
320  const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
321  const int numDiagonals = data.extent_int(2) - blockWidth * blockWidth;
322  const int spaceDim = blockWidth + numDiagonals;
323 
324  PointScalar det = 1.0;
325  switch (blockWidth)
326  {
327  case 0:
328  det = 1.0;
329  break;
330  case 1:
331  {
332  det = data(cellPointOrdinal,0);
333  break;
334  }
335  case 2:
336  {
337  const auto & a = data(cellPointOrdinal,0);
338  const auto & b = data(cellPointOrdinal,1);
339  const auto & c = data(cellPointOrdinal,2);
340  const auto & d = data(cellPointOrdinal,3);
341  det = det2(a,b,c,d);
342 
343  break;
344  }
345  case 3:
346  {
347  const auto & a = data(cellPointOrdinal,0);
348  const auto & b = data(cellPointOrdinal,1);
349  const auto & c = data(cellPointOrdinal,2);
350  const auto & d = data(cellPointOrdinal,3);
351  const auto & e = data(cellPointOrdinal,4);
352  const auto & f = data(cellPointOrdinal,5);
353  const auto & g = data(cellPointOrdinal,6);
354  const auto & h = data(cellPointOrdinal,7);
355  const auto & i = data(cellPointOrdinal,8);
356  det = det3(a,b,c,d,e,f,g,h,i);
357 
358  break;
359  }
360  default:
361  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL");
362  }
363  // incorporate the remaining, diagonal entries
364  const int offset = blockWidth * blockWidth;
365 
366  for (int d=blockWidth; d<spaceDim; d++)
367  {
368  const int index = d-blockWidth+offset;
369  det *= data(cellPointOrdinal,index);
370  }
371  detData(cellPointOrdinal) = det;
372  });
373  }
374  else // neither cell nor point varies
375  {
376  auto data = jacobian.getUnderlyingView1();
377  auto detData = jacobianDet.getUnderlyingView1();
378  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,1),
379  KOKKOS_LAMBDA (const int &dummyArg) {
380  const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
381  const int numDiagonals = jacobian.extent_int(2) - blockWidth * blockWidth;
382  const int spaceDim = blockWidth + numDiagonals;
383 
384  PointScalar det = 1.0;
385  switch (blockWidth)
386  {
387  case 0:
388  det = 1.0;
389  break;
390  case 1:
391  {
392  det = data(0);
393  break;
394  }
395  case 2:
396  {
397  const auto & a = data(0);
398  const auto & b = data(1);
399  const auto & c = data(2);
400  const auto & d = data(3);
401  det = det2(a,b,c,d);
402 
403  break;
404  }
405  case 3:
406  {
407  const auto & a = data(0);
408  const auto & b = data(1);
409  const auto & c = data(2);
410  const auto & d = data(3);
411  const auto & e = data(4);
412  const auto & f = data(5);
413  const auto & g = data(6);
414  const auto & h = data(7);
415  const auto & i = data(8);
416  det = det3(a,b,c,d,e,f,g,h,i);
417 
418  break;
419  }
420  default:
421  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL");
422  }
423  // incorporate the remaining, diagonal entries
424  const int offset = blockWidth * blockWidth;
425 
426  for (int d=blockWidth; d<spaceDim; d++)
427  {
428  const int index = d-blockWidth+offset;
429  det *= data(index);
430  }
431  detData(0) = det;
432  });
433  }
434  }
435  else if ( jacobian.getUnderlyingViewRank() == 4 )
436  {
437  auto data = jacobian.getUnderlyingView4();
438  auto detData = jacobianDet.getUnderlyingView2();
440  }
441  else if ( jacobian.getUnderlyingViewRank() == 3 )
442  {
443  auto data = jacobian.getUnderlyingView3();
444  auto detData = jacobianDet.getUnderlyingView1();
446  }
447  else if ( jacobian.getUnderlyingViewRank() == 2 )
448  {
449  auto data = jacobian.getUnderlyingView2();
450  auto detData = jacobianDet.getUnderlyingView1();
451  Kokkos::parallel_for("fill jacobian det", Kokkos::RangePolicy<ExecSpaceType>(0,1), KOKKOS_LAMBDA(const int &i)
452  {
453  detData(0) = Intrepid2::Kernels::Serial::determinant(data);
454  });
455  }
456  else
457  {
458  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "jacobian's underlying view must have rank 2,3, or 4");
459  }
460  }
461 
462  template<typename DeviceType>
463  template<class PointScalar>
465  {
466  auto variationTypes = jacobian.getVariationTypes();
467 
468  const int CELL_DIM = 0;
469  const int POINT_DIM = 1;
470  const int D1_DIM = 2;
471 
472  auto det2 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c, const PointScalar &d) -> PointScalar
473  {
474  return a * d - b * c;
475  };
476 
477  auto det3 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c,
478  const PointScalar &d, const PointScalar &e, const PointScalar &f,
479  const PointScalar &g, const PointScalar &h, const PointScalar &i) -> PointScalar
480  {
481  return a * det2(e,f,h,i) - b * det2(d,f,g,i) + c * det2(d,e,g,h);
482  };
483 
484  if (variationTypes[D1_DIM] == BLOCK_PLUS_DIAGONAL)
485  {
486  const bool cellVaries = variationTypes[CELL_DIM] != CONSTANT;
487  const bool pointVaries = variationTypes[POINT_DIM] != CONSTANT;
488  if (cellVaries && pointVaries)
489  {
490  auto data = jacobian.getUnderlyingView3();
491  auto invData = jacobianInv.getUnderlyingView3();
492 
493  Kokkos::parallel_for(
494  Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>>({0,0},{data.extent_int(0),data.extent_int(1)}),
495  KOKKOS_LAMBDA (int cellOrdinal, int pointOrdinal) {
496  const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
497  const int numDiagonals = data.extent_int(2) - blockWidth * blockWidth;
498  const int spaceDim = blockWidth + numDiagonals;
499 
500  switch (blockWidth)
501  {
502  case 0:
503  break;
504  case 1:
505  {
506  invData(cellOrdinal,pointOrdinal,0) = 1.0 / data(cellOrdinal,pointOrdinal,0);
507  break;
508  }
509  case 2:
510  {
511  const auto & a = data(cellOrdinal,pointOrdinal,0);
512  const auto & b = data(cellOrdinal,pointOrdinal,1);
513  const auto & c = data(cellOrdinal,pointOrdinal,2);
514  const auto & d = data(cellOrdinal,pointOrdinal,3);
515  const PointScalar det = det2(a,b,c,d);
516 
517  invData(cellOrdinal,pointOrdinal,0) = d/det;
518  invData(cellOrdinal,pointOrdinal,1) = - b/det;
519  invData(cellOrdinal,pointOrdinal,2) = - c/det;
520  invData(cellOrdinal,pointOrdinal,3) = a/det;
521  break;
522  }
523  case 3:
524  {
525  const auto & a = data(cellOrdinal,pointOrdinal,0);
526  const auto & b = data(cellOrdinal,pointOrdinal,1);
527  const auto & c = data(cellOrdinal,pointOrdinal,2);
528  const auto & d = data(cellOrdinal,pointOrdinal,3);
529  const auto & e = data(cellOrdinal,pointOrdinal,4);
530  const auto & f = data(cellOrdinal,pointOrdinal,5);
531  const auto & g = data(cellOrdinal,pointOrdinal,6);
532  const auto & h = data(cellOrdinal,pointOrdinal,7);
533  const auto & i = data(cellOrdinal,pointOrdinal,8);
534  const PointScalar det = det3(a,b,c,d,e,f,g,h,i);
535 
536  {
537  const auto val0 = e*i - h*f;
538  const auto val1 = - d*i + g*f;
539  const auto val2 = d*h - g*e;
540 
541  invData(cellOrdinal,pointOrdinal,0) = val0/det;
542  invData(cellOrdinal,pointOrdinal,1) = val1/det;
543  invData(cellOrdinal,pointOrdinal,2) = val2/det;
544  }
545  {
546  const auto val0 = h*c - b*i;
547  const auto val1 = a*i - g*c;
548  const auto val2 = - a*h + g*b;
549 
550  invData(cellOrdinal,pointOrdinal,3) = val0/det;
551  invData(cellOrdinal,pointOrdinal,4) = val1/det;
552  invData(cellOrdinal,pointOrdinal,5) = val2/det;
553  }
554  {
555  const auto val0 = b*f - e*c;
556  const auto val1 = - a*f + d*c;
557  const auto val2 = a*e - d*b;
558 
559  invData(cellOrdinal,pointOrdinal,6) = val0/det;
560  invData(cellOrdinal,pointOrdinal,7) = val1/det;
561  invData(cellOrdinal,pointOrdinal,8) = val2/det;
562  }
563  break;
564  }
565  default:
566  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL");
567  }
568  // handle the remaining, diagonal entries
569  const int offset = blockWidth * blockWidth;
570 
571  for (int d=blockWidth; d<spaceDim; d++)
572  {
573  const int index = d-blockWidth+offset;
574  invData(cellOrdinal,pointOrdinal,index) = 1.0 / data(cellOrdinal,pointOrdinal,index);
575  }
576  });
577  }
578  else if (cellVaries || pointVaries) // exactly one of cell,point vary -- whichever it is will be in rank 0 of data, invData
579  {
580  auto data = jacobian.getUnderlyingView2();
581  auto invData = jacobianInv.getUnderlyingView2();
582 
583  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,data.extent_int(0)),
584  KOKKOS_LAMBDA (const int &cellPointOrdinal) {
585  const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
586  const int numDiagonals = data.extent_int(1) - blockWidth * blockWidth;
587  const int spaceDim = blockWidth + numDiagonals;
588 
589  switch (blockWidth)
590  {
591  case 0:
592  break;
593  case 1:
594  {
595  invData(cellPointOrdinal,0) = 1.0 / data(cellPointOrdinal,0);
596  break;
597  }
598  case 2:
599  {
600  const auto & a = data(cellPointOrdinal,0);
601  const auto & b = data(cellPointOrdinal,1);
602  const auto & c = data(cellPointOrdinal,2);
603  const auto & d = data(cellPointOrdinal,3);
604  const PointScalar det = det2(a,b,c,d);
605 
606  invData(cellPointOrdinal,0) = d/det;
607  invData(cellPointOrdinal,1) = - b/det;
608  invData(cellPointOrdinal,2) = - c/det;
609  invData(cellPointOrdinal,3) = a/det;
610  break;
611  }
612  case 3:
613  {
614  const auto & a = data(cellPointOrdinal,0);
615  const auto & b = data(cellPointOrdinal,1);
616  const auto & c = data(cellPointOrdinal,2);
617  const auto & d = data(cellPointOrdinal,3);
618  const auto & e = data(cellPointOrdinal,4);
619  const auto & f = data(cellPointOrdinal,5);
620  const auto & g = data(cellPointOrdinal,6);
621  const auto & h = data(cellPointOrdinal,7);
622  const auto & i = data(cellPointOrdinal,8);
623  const PointScalar det = det3(a,b,c,d,e,f,g,h,i);
624 
625  {
626  const auto val0 = e*i - h*f;
627  const auto val1 = - d*i + g*f;
628  const auto val2 = d*h - g*e;
629 
630  invData(cellPointOrdinal,0) = val0/det;
631  invData(cellPointOrdinal,1) = val1/det;
632  invData(cellPointOrdinal,2) = val2/det;
633  }
634  {
635  const auto val0 = h*c - b*i;
636  const auto val1 = a*i - g*c;
637  const auto val2 = - a*h + g*b;
638 
639  invData(cellPointOrdinal,3) = val0/det;
640  invData(cellPointOrdinal,4) = val1/det;
641  invData(cellPointOrdinal,5) = val2/det;
642  }
643  {
644  const auto val0 = b*f - e*c;
645  const auto val1 = - a*f + d*c;
646  const auto val2 = a*e - d*b;
647 
648  invData(cellPointOrdinal,6) = val0/det;
649  invData(cellPointOrdinal,7) = val1/det;
650  invData(cellPointOrdinal,8) = val2/det;
651  }
652  break;
653  }
654  default:
655  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL in setJacobianInv()");
656  }
657  // handle the remaining, diagonal entries
658  const int offset = blockWidth * blockWidth;
659 
660  for (int d=blockWidth; d<spaceDim; d++)
661  {
662  const int index = d-blockWidth+offset;
663  invData(cellPointOrdinal,index) = 1.0 / data(cellPointOrdinal,index);
664  }
665  });
666  }
667  else // neither cell nor point varies
668  {
669  auto data = jacobian.getUnderlyingView1();
670  auto invData = jacobianInv.getUnderlyingView1();
671 
672  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,1),
673  KOKKOS_LAMBDA (const int &dummyArg) {
674  const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
675  const int numDiagonals = data.extent_int(0) - blockWidth * blockWidth;
676  const int spaceDim = blockWidth + numDiagonals;
677 
678  switch (blockWidth)
679  {
680  case 0:
681  break;
682  case 1:
683  {
684  invData(0) = 1.0 / data(0);
685  break;
686  }
687  case 2:
688  {
689  const auto & a = data(0);
690  const auto & b = data(1);
691  const auto & c = data(2);
692  const auto & d = data(3);
693  const PointScalar det = det2(a,b,c,d);
694 
695  invData(0) = d/det;
696  invData(1) = - b/det;
697  invData(2) = - c/det;
698  invData(3) = a/det;
699  break;
700  }
701  case 3:
702  {
703  const auto & a = data(0);
704  const auto & b = data(1);
705  const auto & c = data(2);
706  const auto & d = data(3);
707  const auto & e = data(4);
708  const auto & f = data(5);
709  const auto & g = data(6);
710  const auto & h = data(7);
711  const auto & i = data(8);
712  const PointScalar det = det3(a,b,c,d,e,f,g,h,i);
713 
714  {
715  const auto val0 = e*i - h*f;
716  const auto val1 = - d*i + g*f;
717  const auto val2 = d*h - g*e;
718 
719  invData(0) = val0/det;
720  invData(1) = val1/det;
721  invData(2) = val2/det;
722  }
723  {
724  const auto val0 = h*c - b*i;
725  const auto val1 = a*i - g*c;
726  const auto val2 = - a*h + g*b;
727 
728  invData(3) = val0/det;
729  invData(4) = val1/det;
730  invData(5) = val2/det;
731  }
732  {
733  const auto val0 = b*f - e*c;
734  const auto val1 = - a*f + d*c;
735  const auto val2 = a*e - d*b;
736 
737  invData(6) = val0/det;
738  invData(7) = val1/det;
739  invData(8) = val2/det;
740  }
741  break;
742  }
743  default:
744  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL in setJacobianInv()");
745  }
746  // handle the remaining, diagonal entries
747  const int offset = blockWidth * blockWidth;
748 
749  for (int d=blockWidth; d<spaceDim; d++)
750  {
751  const int index = d-blockWidth+offset;
752  invData(index) = 1.0 / data(index);
753  }
754  });
755  }
756  }
757  else if ( jacobian.getUnderlyingViewRank() == 4 )
758  {
759  auto data = jacobian.getUnderlyingView4();
760  auto invData = jacobianInv.getUnderlyingView4();
761 
763  }
764  else if ( jacobian.getUnderlyingViewRank() == 3 )
765  {
766  auto data = jacobian.getUnderlyingView3();
767  auto invData = jacobianInv.getUnderlyingView3();
768 
770  }
771  else if ( jacobian.getUnderlyingViewRank() == 2 ) // Cell, point constant; D1, D2 vary normally
772  {
773  auto data = jacobian.getUnderlyingView2();
774  auto invData = jacobianInv.getUnderlyingView2();
775 
776  Kokkos::parallel_for("fill jacobian inverse", Kokkos::RangePolicy<ExecSpaceType>(0,1), KOKKOS_LAMBDA(const int &i)
777  {
778  Intrepid2::Kernels::Serial::inverse(invData,data);
779  });
780  }
781  else
782  {
783  INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "jacobian's underlying view must have rank 2,3, or 4");
784  }
785  }
786 
787  template<typename DeviceType>
788  template<typename jacobianValueType, class ...jacobianProperties,
789  typename BasisGradientsType,
790  typename WorksetType>
791  void
793  setJacobian( Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian,
794  const WorksetType worksetCell,
795  const BasisGradientsType gradients, const int startCell, const int endCell)
796  {
797  constexpr bool is_accessible =
798  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
799  typename decltype(jacobian)::memory_space>::accessible;
800  static_assert(is_accessible, "CellTools<DeviceType>::setJacobian(..): output view's memory space is not compatible with DeviceType");
801 
802  using JacobianViewType = Kokkos::DynRankView<jacobianValueType,jacobianProperties...>;
804 
805  // resolve the -1 default argument for endCell into the true end cell index
806  int endCellResolved = (endCell == -1) ? worksetCell.extent_int(0) : endCell;
807 
808  using range_policy_type = Kokkos::Experimental::MDRangePolicy
809  < ExecSpaceType, Kokkos::Experimental::Rank<2>, Kokkos::IndexType<ordinal_type> >;
810  range_policy_type policy( { 0, 0 },
811  { jacobian.extent(0), jacobian.extent(1) } );
812  Kokkos::parallel_for( policy, FunctorType(jacobian, worksetCell, gradients, startCell, endCellResolved) );
813  }
814 
815  template<typename DeviceType>
816  template<typename jacobianValueType, class ...jacobianProperties,
817  typename pointValueType, class ...pointProperties,
818  typename WorksetType,
819  typename HGradBasisType>
820  void
822  setJacobian( Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian,
823  const Kokkos::DynRankView<pointValueType,pointProperties...> points,
824  const WorksetType worksetCell,
825  const Teuchos::RCP<HGradBasisType> basis,
826  const int startCell, const int endCell) {
827  constexpr bool are_accessible =
828  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
829  typename decltype(jacobian)::memory_space>::accessible &&
830  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
831  typename decltype(points)::memory_space>::accessible;
832  static_assert(are_accessible, "CellTools<DeviceType>::setJacobian(..): input/output views' memory spaces are not compatible with DeviceType");
833 
834 #ifdef HAVE_INTREPID2_DEBUG
835  CellTools_setJacobianArgs(jacobian, points, worksetCell, basis->getBaseCellTopology(), startCell, endCell);
836  //static_assert(std::is_same( pointValueType, decltype(basis->getDummyOutputValue()) ));
837 #endif
838  const auto cellTopo = basis->getBaseCellTopology();
839  const ordinal_type spaceDim = cellTopo.getDimension();
840  const ordinal_type numCells = jacobian.extent(0);
841 
842  //points can be rank-2 (P,D), or rank-3 (C,P,D)
843  const ordinal_type pointRank = points.rank();
844  const ordinal_type numPoints = (pointRank == 2 ? points.extent(0) : points.extent(1));
845  const ordinal_type basisCardinality = basis->getCardinality();
846 
847  // the following does not work for RCP; its * operator returns reference not the object
848  //typedef typename decltype(*basis)::output_value_type gradValueType;
849  //typedef Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),DeviceType> gradViewType;
850 
851  auto vcprop = Kokkos::common_view_alloc_prop(points);
852  using GradViewType = Kokkos::DynRankView<typename decltype(vcprop)::value_type,DeviceType>;
853 
854  GradViewType grads;
855 
856  switch (pointRank) {
857  case 2: {
858  // For most FEMs
859  grads = GradViewType(Kokkos::view_alloc("CellTools::setJacobian::grads", vcprop),basisCardinality, numPoints, spaceDim);
860  basis->getValues(grads,
861  points,
862  OPERATOR_GRAD);
863  break;
864  }
865  case 3: {
866  // For CVFEM
867  grads = GradViewType(Kokkos::view_alloc("CellTools::setJacobian::grads", vcprop), numCells, basisCardinality, numPoints, spaceDim);
868  for (ordinal_type cell=0;cell<numCells;++cell)
869  basis->getValues(Kokkos::subview( grads, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() ),
870  Kokkos::subview( points, cell, Kokkos::ALL(), Kokkos::ALL() ),
871  OPERATOR_GRAD);
872  break;
873  }
874  }
875 
876  setJacobian(jacobian, worksetCell, grads, startCell, endCell);
877  }
878 
879  template<typename DeviceType>
880  template<typename jacobianInvValueType, class ...jacobianInvProperties,
881  typename jacobianValueType, class ...jacobianProperties>
882  void
884  setJacobianInv( Kokkos::DynRankView<jacobianInvValueType,jacobianInvProperties...> jacobianInv,
885  const Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian ) {
886 #ifdef HAVE_INTREPID2_DEBUG
887  CellTools_setJacobianInvArgs(jacobianInv, jacobian);
888 #endif
889  RealSpaceTools<DeviceType>::inverse(jacobianInv, jacobian);
890  }
891 
892  template<typename DeviceType>
893  template<typename jacobianDetValueType, class ...jacobianDetProperties,
894  typename jacobianValueType, class ...jacobianProperties>
895  void
897  setJacobianDet( Kokkos::DynRankView<jacobianDetValueType,jacobianDetProperties...> jacobianDet,
898  const Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian ) {
899 #ifdef HAVE_INTREPID2_DEBUG
900  CellTools_setJacobianDetArgs(jacobianDet, jacobian);
901 #endif
902  RealSpaceTools<DeviceType>::det(jacobianDet, jacobian);
903  }
904 
905 }
906 
907 #endif
static void det(DeterminantArrayViewType detArray, const MatrixViewType inMats)
Computes determinants of matrices stored in an array of total rank 3 (array of matrices), indexed by (i0, D, D), or 4 (array of arrays of matrices), indexed by (i0, i1, D, D).
KOKKOS_INLINE_FUNCTION Kokkos::Array< int, 7 > getExtents() const
Returns an array containing the logical extents in each dimension.
KOKKOS_INLINE_FUNCTION ordinal_type getUnderlyingViewRank() const
returns the rank of the View that stores the unique data
KOKKOS_INLINE_FUNCTION const int & blockPlusDiagonalLastNonDiagonal() const
For a Data object containing data with variation type BLOCK_PLUS_DIAGONAL, returns the row and column...
one of two dimensions in a matrix; bottom-right part of matrix is diagonal
Functor for calculation of Jacobian on cell workset see Intrepid2::CellTools for more.
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar *, DeviceType > & getUnderlyingView1() const
returns the View that stores the unique data. For rank-1 underlying containers.
static void inverse(InverseMatrixViewType inverseMats, MatrixViewType inMats)
Computes inverses of nonsingular matrices stored in an array of total rank 2 (single matrix)...
void CellTools_setJacobianDetArgs(const jacobianDetViewType jacobianDet, const jacobianViewType jacobian)
Validates arguments to Intrepid2::CellTools::setJacobianDet.
KOKKOS_INLINE_FUNCTION int extent_int(const int &r) const
Returns the logical extent in the specified dimension.
static Data< PointScalar, DeviceType > allocateJacobianInv(const Data< PointScalar, DeviceType > &jacobian)
Allocates and returns a Data container suitable for storing inverses corresponding to the Jacobians i...
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar ****, DeviceType > & getUnderlyingView4() const
returns the View that stores the unique data. For rank-4 underlying containers.
Kokkos::DynRankView< typename ViewType::value_type, typename DeduceLayout< ViewType >::result_layout, typename ViewType::device_type > getMatchingViewWithLabel(const ViewType &view, const std::string &label, DimArgs... dims)
Creates and returns a view that matches the provided view in Kokkos Layout.
void CellTools_setJacobianInvArgs(const jacobianInvViewType jacobianInv, const jacobianViewType jacobian)
Validates arguments to Intrepid2::CellTools::setJacobianInv.
static void setJacobianInv(Kokkos::DynRankView< jacobianInvValueType, jacobianInvProperties... > jacobianInv, const Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F...
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar **, DeviceType > & getUnderlyingView2() const
returns the View that stores the unique data. For rank-2 underlying containers.
Header file for small functions used in Intrepid2.
static Data< PointScalar, DeviceType > allocateJacobianDet(const Data< PointScalar, DeviceType > &jacobian)
Allocates and returns a Data container suitable for storing determinants corresponding to the Jacobia...
static void setJacobianDet(Kokkos::DynRankView< jacobianDetValueType, jacobianDetProperties... > jacobianDet, const Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F...
KOKKOS_INLINE_FUNCTION const Kokkos::View< DataScalar ***, DeviceType > & getUnderlyingView3() const
returns the View that stores the unique data. For rank-3 underlying containers.
static void setJacobian(Kokkos::DynRankView< jacobianValueType, jacobianProperties... > jacobian, const Kokkos::DynRankView< pointValueType, pointProperties... > points, const WorksetType worksetCell, const Teuchos::RCP< HGradBasisType > basis, const int startCell=0, const int endCell=-1)
Computes the Jacobian matrix DF of the reference-to-physical frame map F.
KOKKOS_INLINE_FUNCTION const Kokkos::Array< DataVariationType, 7 > & getVariationTypes() const
Returns an array with the variation types in each logical dimension.