Intrepid2
Intrepid2_CellToolsDefJacobian.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10
16#ifndef __INTREPID2_CELLTOOLS_DEF_HPP__
17#define __INTREPID2_CELLTOOLS_DEF_HPP__
18
19// disable clang warnings
20#if defined (__clang__) && !defined (__INTEL_COMPILER)
21#pragma clang system_header
22#endif
23
24#include "Intrepid2_Kernels.hpp"
25#include "Intrepid2_DataTools.hpp"
26
27namespace Intrepid2 {
28
29
30 //============================================================================================//
31 // //
32 // Jacobian, inverse Jacobian and Jacobian determinant //
33 // //
34 //============================================================================================//
35
36 namespace FunctorCellTools {
40 template<typename jacobianViewType,
41 typename worksetCellType,
42 typename basisGradType,
43 bool squareJacobian = false>
44 struct F_setJacobian {
45 jacobianViewType _jacobian;
46 const worksetCellType _worksetCells;
47 const basisGradType _basisGrads;
48 const int _startCell;
49
50 KOKKOS_INLINE_FUNCTION
51 F_setJacobian( jacobianViewType jacobian_,
52 worksetCellType worksetCells_,
53 basisGradType basisGrads_,
54 const int startCell_)
55 : _jacobian(jacobian_), _worksetCells(worksetCells_), _basisGrads(basisGrads_),
56 _startCell(startCell_) {}
57
58 KOKKOS_INLINE_FUNCTION
59 void operator()(const ordinal_type cell,
60 const ordinal_type point) const {
61
62 const ordinal_type phys_dim = _jacobian.extent(2); // dim2 and dim3 should match
63
64 const ordinal_type gradRank = rank(_basisGrads);
65
66 if constexpr (squareJacobian)
67 {
68 const ordinal_type dim = phys_dim;
69
70 if ( gradRank == 3)
71 {
72 const ordinal_type cardinality = _basisGrads.extent(0);
73 for (ordinal_type i=0;i<dim;++i)
74 for (ordinal_type j=0;j<dim;++j) {
75 _jacobian(cell, point, i, j) = 0;
76 for (ordinal_type bf=0;bf<cardinality;++bf)
77 _jacobian(cell, point, i, j) += _worksetCells(cell+_startCell, bf, i) * _basisGrads(bf, point, j);
78 }
79 }
80 else
81 {
82 const ordinal_type cardinality = _basisGrads.extent(1);
83 for (ordinal_type i=0;i<dim;++i)
84 for (ordinal_type j=0;j<dim;++j) {
85 _jacobian(cell, point, i, j) = 0;
86 for (ordinal_type bf=0;bf<cardinality;++bf)
87 _jacobian(cell, point, i, j) += _worksetCells(cell+_startCell, bf, i) * _basisGrads(cell, bf, point, j);
88 }
89 }
90 }
91 else
92 {
93 const ordinal_type ref_dim = _basisGrads.extent_int(gradRank-1);
94
95 if ( gradRank == 3)
96 {
97 const ordinal_type cardinality = _basisGrads.extent(0);
98 for (ordinal_type i=0;i<phys_dim;++i)
99 for (ordinal_type j=0;j<ref_dim;++j) {
100 _jacobian(cell, point, i, j) = 0;
101 for (ordinal_type bf=0;bf<cardinality;++bf)
102 _jacobian(cell, point, i, j) += _worksetCells(cell+_startCell, bf, i) * _basisGrads(bf, point, j);
103 }
104 }
105 else
106 {
107 const ordinal_type cardinality = _basisGrads.extent(1);
108 for (ordinal_type i=0;i<phys_dim;++i)
109 for (ordinal_type j=0;j<ref_dim;++j) {
110 _jacobian(cell, point, i, j) = 0;
111 for (ordinal_type bf=0;bf<cardinality;++bf)
112 _jacobian(cell, point, i, j) += _worksetCells(cell+_startCell, bf, i) * _basisGrads(cell, bf, point, j);
113 }
114 }
115 }
116 }
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 = Impl::createMatchingDynRankView(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 = Impl::createMatchingDynRankView(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 = Impl::createMatchingDynRankView(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 = Impl::createMatchingDynRankView(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 = Impl::createMatchingDynRankView(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 = Impl::createMatchingDynRankView(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 = Impl::createMatchingDynRankView(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::SpaceAccessibility<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 SquareFunctorType = FunctorCellTools::F_setJacobian<JacobianViewType,WorksetType,BasisGradientsType,true> ;
814 using RectangularFunctorType = FunctorCellTools::F_setJacobian<JacobianViewType,WorksetType,BasisGradientsType,false> ;
815
816 // resolve the -1 default argument for endCell into the true end cell index
817 // In some cases, endCell may be smaller than worksetCell.extent(0).
818 // This can occur when cells are split across worksets and the final
819 // workset is not completely filled with valid cells.
820 ordinal_type endCellResolved = (endCell == -1) ? worksetCell.extent_int(0) : endCell;
821
822 using range_policy_type = Kokkos::MDRangePolicy
823 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
824 range_policy_type policy( { 0, 0 },
825 { std::min<ordinal_type>(static_cast<ordinal_type>(jacobian.extent_int(0)), endCellResolved-startCell), jacobian.extent_int(1) } );
826
827 const ordinal_type gradRank = rank(gradients);
828 const ordinal_type physDim = jacobian.extent_int(2);
829 const ordinal_type refDim = gradients.extent_int(gradRank-1);
830
831 if (physDim == refDim)
832 {
833 Kokkos::parallel_for( policy, SquareFunctorType(jacobian, worksetCell, gradients, startCell) );
834 }
835 else
836 {
837 Kokkos::parallel_for( policy, RectangularFunctorType(jacobian, worksetCell, gradients, startCell) );
838 }
839 }
840
841
842 template<typename DeviceType>
843 template<typename JacobianViewType,
844 typename PointViewType,
845 typename WorksetType,
846 typename HGradBasisType>
847 void
848 CellTools<DeviceType>::
849 setJacobian( JacobianViewType jacobian,
850 const PointViewType points,
851 const WorksetType worksetCell,
852 const Teuchos::RCP<HGradBasisType> basis,
853 const int startCell, const int endCell) {
854 constexpr bool are_accessible =
855 Kokkos::SpaceAccessibility<MemSpaceType,
856 typename decltype(jacobian)::memory_space>::accessible &&
857 Kokkos::SpaceAccessibility<MemSpaceType,
858 typename decltype(points)::memory_space>::accessible;
859 static_assert(are_accessible, "CellTools<DeviceType>::setJacobian(..): input/output views' memory spaces are not compatible with DeviceType");
860
861#ifdef HAVE_INTREPID2_DEBUG
862 CellTools_setJacobianArgs(jacobian, points, worksetCell, basis->getBaseCellTopology(), startCell, endCell);
863 //static_assert(std::is_same( pointValueType, decltype(basis->getDummyOutputValue()) ));
864#endif
865 const ordinal_type numCells = jacobian.extent(0);
866
867 //points can be rank-2 (P,D), or rank-3 (C,P,D)
868 const ordinal_type pointRank = points.rank();
869 const ordinal_type numPoints = (pointRank == 2 ? points.extent(0) : points.extent(1));
870 const ordinal_type basisCardinality = basis->getCardinality();
871
872 using GradViewType = Kokkos::DynRankView<typename decltype(points)::value_type,DeviceType>;
873
874 GradViewType grads;
875
876 switch (pointRank) {
877 case 2: {
878 // For most FEMs
879 grads = Impl::createMatchingView<GradViewType>(points, "CellTools::setJacobian::grads", basisCardinality, numPoints, basis->getDomainDimension());
880 basis->getValues(grads,
881 points,
882 OPERATOR_GRAD);
883 break;
884 }
885 case 3: {
886 // For CVFEM or map to reference frames for point inclusion
887 grads = Impl::createMatchingView<GradViewType>(points, "CellTools::setJacobian::grads", numCells, basisCardinality, numPoints, basis->getDomainDimension());
888 getHGradValues<OPERATOR_GRAD>(grads,points,basis.get());
889 break;
890 }
891 }
892
893 setJacobian(jacobian, worksetCell, grads, startCell, endCell);
894 }
895
896 template<typename DeviceType>
897 template<typename JacobianInvViewType,
898 typename JacobianViewType>
899 void
900 CellTools<DeviceType>::
901 setJacobianInv( JacobianInvViewType jacobianInv,
902 const JacobianViewType jacobian ) {
903#ifdef HAVE_INTREPID2_DEBUG
904 CellTools_setJacobianInvArgs(jacobianInv, jacobian);
905#endif
906 RealSpaceTools<DeviceType>::inverse(jacobianInv, jacobian);
907 }
908
909 template<typename DeviceType>
910 template<typename JacobianDetViewType,
911 typename JacobianViewType>
912 void
913 CellTools<DeviceType>::
914 setJacobianDet( JacobianDetViewType jacobianDet,
915 const JacobianViewType jacobian ) {
916#ifdef HAVE_INTREPID2_DEBUG
917 CellTools_setJacobianDetArgs(jacobianDet, jacobian);
918#endif
919 RealSpaceTools<DeviceType>::det(jacobianDet, jacobian);
920 }
921
922 template<typename DeviceType>
923 template<typename Scalar>
924 void
925 CellTools<DeviceType>::
926 setJacobianDividedByDet( Data<Scalar,DeviceType> & jacobianDividedByDet,
927 const Data<Scalar,DeviceType> & jacobian,
928 const Data<Scalar,DeviceType> & jacobianDetInv)
929 {
930 DataTools::multiplyByCPWeights(jacobianDividedByDet,jacobian,jacobianDetInv);
931 }
932}
933
934#endif
Intrepid2::CellTools_setJacobianInvArgs
void CellTools_setJacobianInvArgs(const jacobianInvViewType jacobianInv, const jacobianViewType jacobian)
Validates arguments to Intrepid2::CellTools::setJacobianInv.
Definition Intrepid2_CellToolsDefValidateArguments.hpp:116
Intrepid2::CellTools_setJacobianDetArgs
void CellTools_setJacobianDetArgs(const jacobianDetViewType jacobianDet, const jacobianViewType jacobian)
Validates arguments to Intrepid2::CellTools::setJacobianDet.
Definition Intrepid2_CellToolsDefValidateArguments.hpp:148
Intrepid2_DataTools.hpp
Utility methods for manipulating Intrepid2::Data objects.
Intrepid2::CONSTANT
@ CONSTANT
does not vary
Definition Intrepid2_DataVariationType.hpp:40
Intrepid2::BLOCK_PLUS_DIAGONAL
@ BLOCK_PLUS_DIAGONAL
one of two dimensions in a matrix; bottom-right part of matrix is diagonal
Definition Intrepid2_DataVariationType.hpp:42
Intrepid2_Kernels.hpp
Header file for small functions used in Intrepid2.
INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE
#define INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(test, x, msg)
Definition Intrepid2_Utils.hpp:70
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,...
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::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:914
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:901
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::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:849
Intrepid2::CellTools::setJacobianDetInv
static void setJacobianDetInv(Data< PointScalar, DeviceType > &jacobianDetInv, const Data< PointScalar, DeviceType > &jacobian)
Computes reciprocals of determinants corresponding to the Jacobians in the Data container provided.
Definition Intrepid2_CellToolsDefJacobian.hpp:465
Intrepid2::Data
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
Definition Intrepid2_Data.hpp:79
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:1216
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:991
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:1005
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:998
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:858
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:1763
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:1168
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:984
Intrepid2::Data::allocateConstantData
Data< DataScalar, DeviceType > allocateConstantData(const DataScalar &value)
Definition Intrepid2_Data.hpp:1248
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:865
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:1049
Intrepid2::DataTools::multiplyByCPWeights
static void multiplyByCPWeights(Data< Scalar, DeviceType > &resultMatrixData, const Data< Scalar, DeviceType > &matrixDataIn, const Data< Scalar, DeviceType > &scalarDataIn)
Definition Intrepid2_DataTools.hpp:31
Intrepid2::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition Intrepid2_RealSpaceTools.hpp:46
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::FunctorCellTools::F_setJacobian
Functor for calculation of Jacobian on cell workset see Intrepid2::CellTools for more.
Definition Intrepid2_CellToolsDefJacobian.hpp:44
Generated by doxygen 1.9.8