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