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