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