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