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>
75 struct F_setJacobian {
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>
124 Data<PointScalar,DeviceType> CellTools<DeviceType>::allocateJacobianDet( const Data<PointScalar,DeviceType> & jacobian )
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>
158 Data<PointScalar,DeviceType> CellTools<DeviceType>::allocateJacobianInv( const Data<PointScalar,DeviceType> & jacobian )
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>
197 void CellTools<DeviceType>::setJacobianDet( Data<PointScalar,DeviceType> &jacobianDet, const Data<PointScalar,DeviceType> & jacobian )
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();
439 Intrepid2::RealSpaceTools<ExecSpaceType>::det(detData, data);
440 }
441 else if ( jacobian.getUnderlyingViewRank() == 3 )
442 {
443 auto data = jacobian.getUnderlyingView3();
444 auto detData = jacobianDet.getUnderlyingView1();
445 Intrepid2::RealSpaceTools<ExecSpaceType>::det(detData, data);
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>
464 void CellTools<DeviceType>::setJacobianInv( Data<PointScalar,DeviceType> &jacobianInv, const Data<PointScalar,DeviceType> & jacobian )
465 {
466 auto variationTypes = jacobian.getVariationTypes();
467
468 const int CELL_DIM = 0;
469 const int POINT_DIM = 1;
470 const int D1_DIM = 2;
471
472 auto det2 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c, const PointScalar &d) -> PointScalar
473 {
474 return a * d - b * c;
475 };
476
477 auto det3 = KOKKOS_LAMBDA (const PointScalar &a, const PointScalar &b, const PointScalar &c,
478 const PointScalar &d, const PointScalar &e, const PointScalar &f,
479 const PointScalar &g, const PointScalar &h, const PointScalar &i) -> PointScalar
480 {
481 return a * det2(e,f,h,i) - b * det2(d,f,g,i) + c * det2(d,e,g,h);
482 };
483
484 if (variationTypes[D1_DIM] == BLOCK_PLUS_DIAGONAL)
485 {
486 const bool cellVaries = variationTypes[CELL_DIM] != CONSTANT;
487 const bool pointVaries = variationTypes[POINT_DIM] != CONSTANT;
488 if (cellVaries && pointVaries)
489 {
490 auto data = jacobian.getUnderlyingView3();
491 auto invData = jacobianInv.getUnderlyingView3();
492
493 Kokkos::parallel_for(
494 Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>>({0,0},{data.extent_int(0),data.extent_int(1)}),
495 KOKKOS_LAMBDA (int cellOrdinal, int pointOrdinal) {
496 const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
497 const int numDiagonals = data.extent_int(2) - blockWidth * blockWidth;
498 const int spaceDim = blockWidth + numDiagonals;
499
500 switch (blockWidth)
501 {
502 case 0:
503 break;
504 case 1:
505 {
506 invData(cellOrdinal,pointOrdinal,0) = 1.0 / data(cellOrdinal,pointOrdinal,0);
507 break;
508 }
509 case 2:
510 {
511 const auto & a = data(cellOrdinal,pointOrdinal,0);
512 const auto & b = data(cellOrdinal,pointOrdinal,1);
513 const auto & c = data(cellOrdinal,pointOrdinal,2);
514 const auto & d = data(cellOrdinal,pointOrdinal,3);
515 const PointScalar det = det2(a,b,c,d);
516
517 invData(cellOrdinal,pointOrdinal,0) = d/det;
518 invData(cellOrdinal,pointOrdinal,1) = - b/det;
519 invData(cellOrdinal,pointOrdinal,2) = - c/det;
520 invData(cellOrdinal,pointOrdinal,3) = a/det;
521 break;
522 }
523 case 3:
524 {
525 const auto & a = data(cellOrdinal,pointOrdinal,0);
526 const auto & b = data(cellOrdinal,pointOrdinal,1);
527 const auto & c = data(cellOrdinal,pointOrdinal,2);
528 const auto & d = data(cellOrdinal,pointOrdinal,3);
529 const auto & e = data(cellOrdinal,pointOrdinal,4);
530 const auto & f = data(cellOrdinal,pointOrdinal,5);
531 const auto & g = data(cellOrdinal,pointOrdinal,6);
532 const auto & h = data(cellOrdinal,pointOrdinal,7);
533 const auto & i = data(cellOrdinal,pointOrdinal,8);
534 const PointScalar det = det3(a,b,c,d,e,f,g,h,i);
535
536 {
537 const auto val0 = e*i - h*f;
538 const auto val1 = - d*i + g*f;
539 const auto val2 = d*h - g*e;
540
541 invData(cellOrdinal,pointOrdinal,0) = val0/det;
542 invData(cellOrdinal,pointOrdinal,1) = val1/det;
543 invData(cellOrdinal,pointOrdinal,2) = val2/det;
544 }
545 {
546 const auto val0 = h*c - b*i;
547 const auto val1 = a*i - g*c;
548 const auto val2 = - a*h + g*b;
549
550 invData(cellOrdinal,pointOrdinal,3) = val0/det;
551 invData(cellOrdinal,pointOrdinal,4) = val1/det;
552 invData(cellOrdinal,pointOrdinal,5) = val2/det;
553 }
554 {
555 const auto val0 = b*f - e*c;
556 const auto val1 = - a*f + d*c;
557 const auto val2 = a*e - d*b;
558
559 invData(cellOrdinal,pointOrdinal,6) = val0/det;
560 invData(cellOrdinal,pointOrdinal,7) = val1/det;
561 invData(cellOrdinal,pointOrdinal,8) = val2/det;
562 }
563 break;
564 }
565 default:
566 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL");
567 }
568 // handle the remaining, diagonal entries
569 const int offset = blockWidth * blockWidth;
570
571 for (int d=blockWidth; d<spaceDim; d++)
572 {
573 const int index = d-blockWidth+offset;
574 invData(cellOrdinal,pointOrdinal,index) = 1.0 / data(cellOrdinal,pointOrdinal,index);
575 }
576 });
577 }
578 else if (cellVaries || pointVaries) // exactly one of cell,point vary -- whichever it is will be in rank 0 of data, invData
579 {
580 auto data = jacobian.getUnderlyingView2();
581 auto invData = jacobianInv.getUnderlyingView2();
582
583 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,data.extent_int(0)),
584 KOKKOS_LAMBDA (const int &cellPointOrdinal) {
585 const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
586 const int numDiagonals = data.extent_int(1) - blockWidth * blockWidth;
587 const int spaceDim = blockWidth + numDiagonals;
588
589 switch (blockWidth)
590 {
591 case 0:
592 break;
593 case 1:
594 {
595 invData(cellPointOrdinal,0) = 1.0 / data(cellPointOrdinal,0);
596 break;
597 }
598 case 2:
599 {
600 const auto & a = data(cellPointOrdinal,0);
601 const auto & b = data(cellPointOrdinal,1);
602 const auto & c = data(cellPointOrdinal,2);
603 const auto & d = data(cellPointOrdinal,3);
604 const PointScalar det = det2(a,b,c,d);
605
606 invData(cellPointOrdinal,0) = d/det;
607 invData(cellPointOrdinal,1) = - b/det;
608 invData(cellPointOrdinal,2) = - c/det;
609 invData(cellPointOrdinal,3) = a/det;
610 break;
611 }
612 case 3:
613 {
614 const auto & a = data(cellPointOrdinal,0);
615 const auto & b = data(cellPointOrdinal,1);
616 const auto & c = data(cellPointOrdinal,2);
617 const auto & d = data(cellPointOrdinal,3);
618 const auto & e = data(cellPointOrdinal,4);
619 const auto & f = data(cellPointOrdinal,5);
620 const auto & g = data(cellPointOrdinal,6);
621 const auto & h = data(cellPointOrdinal,7);
622 const auto & i = data(cellPointOrdinal,8);
623 const PointScalar det = det3(a,b,c,d,e,f,g,h,i);
624
625 {
626 const auto val0 = e*i - h*f;
627 const auto val1 = - d*i + g*f;
628 const auto val2 = d*h - g*e;
629
630 invData(cellPointOrdinal,0) = val0/det;
631 invData(cellPointOrdinal,1) = val1/det;
632 invData(cellPointOrdinal,2) = val2/det;
633 }
634 {
635 const auto val0 = h*c - b*i;
636 const auto val1 = a*i - g*c;
637 const auto val2 = - a*h + g*b;
638
639 invData(cellPointOrdinal,3) = val0/det;
640 invData(cellPointOrdinal,4) = val1/det;
641 invData(cellPointOrdinal,5) = val2/det;
642 }
643 {
644 const auto val0 = b*f - e*c;
645 const auto val1 = - a*f + d*c;
646 const auto val2 = a*e - d*b;
647
648 invData(cellPointOrdinal,6) = val0/det;
649 invData(cellPointOrdinal,7) = val1/det;
650 invData(cellPointOrdinal,8) = val2/det;
651 }
652 break;
653 }
654 default:
655 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL in setJacobianInv()");
656 }
657 // handle the remaining, diagonal entries
658 const int offset = blockWidth * blockWidth;
659
660 for (int d=blockWidth; d<spaceDim; d++)
661 {
662 const int index = d-blockWidth+offset;
663 invData(cellPointOrdinal,index) = 1.0 / data(cellPointOrdinal,index);
664 }
665 });
666 }
667 else // neither cell nor point varies
668 {
669 auto data = jacobian.getUnderlyingView1();
670 auto invData = jacobianInv.getUnderlyingView1();
671
672 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,1),
673 KOKKOS_LAMBDA (const int &dummyArg) {
674 const int blockWidth = jacobian.blockPlusDiagonalLastNonDiagonal() + 1;
675 const int numDiagonals = data.extent_int(0) - blockWidth * blockWidth;
676 const int spaceDim = blockWidth + numDiagonals;
677
678 switch (blockWidth)
679 {
680 case 0:
681 break;
682 case 1:
683 {
684 invData(0) = 1.0 / data(0);
685 break;
686 }
687 case 2:
688 {
689 const auto & a = data(0);
690 const auto & b = data(1);
691 const auto & c = data(2);
692 const auto & d = data(3);
693 const PointScalar det = det2(a,b,c,d);
694
695 invData(0) = d/det;
696 invData(1) = - b/det;
697 invData(2) = - c/det;
698 invData(3) = a/det;
699 break;
700 }
701 case 3:
702 {
703 const auto & a = data(0);
704 const auto & b = data(1);
705 const auto & c = data(2);
706 const auto & d = data(3);
707 const auto & e = data(4);
708 const auto & f = data(5);
709 const auto & g = data(6);
710 const auto & h = data(7);
711 const auto & i = data(8);
712 const PointScalar det = det3(a,b,c,d,e,f,g,h,i);
713
714 {
715 const auto val0 = e*i - h*f;
716 const auto val1 = - d*i + g*f;
717 const auto val2 = d*h - g*e;
718
719 invData(0) = val0/det;
720 invData(1) = val1/det;
721 invData(2) = val2/det;
722 }
723 {
724 const auto val0 = h*c - b*i;
725 const auto val1 = a*i - g*c;
726 const auto val2 = - a*h + g*b;
727
728 invData(3) = val0/det;
729 invData(4) = val1/det;
730 invData(5) = val2/det;
731 }
732 {
733 const auto val0 = b*f - e*c;
734 const auto val1 = - a*f + d*c;
735 const auto val2 = a*e - d*b;
736
737 invData(6) = val0/det;
738 invData(7) = val1/det;
739 invData(8) = val2/det;
740 }
741 break;
742 }
743 default:
744 INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(true, std::invalid_argument, "blocks with block width > 3 not supported for BLOCK_PLUS_DIAGONAL in setJacobianInv()");
745 }
746 // handle the remaining, diagonal entries
747 const int offset = blockWidth * blockWidth;
748
749 for (int d=blockWidth; d<spaceDim; d++)
750 {
751 const int index = d-blockWidth+offset;
752 invData(index) = 1.0 / data(index);
753 }
754 });
755 }
756 }
757 else if ( jacobian.getUnderlyingViewRank() == 4 )
758 {
759 auto data = jacobian.getUnderlyingView4();
760 auto invData = jacobianInv.getUnderlyingView4();
761
762 Intrepid2::RealSpaceTools<ExecSpaceType>::inverse(invData, data);
763 }
764 else if ( jacobian.getUnderlyingViewRank() == 3 )
765 {
766 auto data = jacobian.getUnderlyingView3();
767 auto invData = jacobianInv.getUnderlyingView3();
768
769 Intrepid2::RealSpaceTools<ExecSpaceType>::inverse(invData, data);
770 }
771 else if ( jacobian.getUnderlyingViewRank() == 2 ) // Cell, point constant; D1, D2 vary normally
772 {
773 auto data = jacobian.getUnderlyingView2();
774 auto invData = jacobianInv.getUnderlyingView2();
775
776 Kokkos::parallel_for("fill jacobian inverse", Kokkos::RangePolicy<ExecSpaceType>(0,1), KOKKOS_LAMBDA(const int &i)
777 {
778 Intrepid2::Kernels::Serial::inverse(invData,data);
779 });
780 }
781 else
782 {
783 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "jacobian's underlying view must have rank 2,3, or 4");
784 }
785 }
786
787 template<typename DeviceType>
788 template<typename jacobianValueType, class ...jacobianProperties,
789 typename BasisGradientsType,
790 typename WorksetType>
791 void
792 CellTools<DeviceType>::
793 setJacobian( Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian,
794 const WorksetType worksetCell,
795 const BasisGradientsType gradients, const int startCell, const int endCell)
796 {
797 constexpr bool is_accessible =
798 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
799 typename decltype(jacobian)::memory_space>::accessible;
800 static_assert(is_accessible, "CellTools<DeviceType>::setJacobian(..): output view's memory space is not compatible with DeviceType");
801
802 using JacobianViewType = Kokkos::DynRankView<jacobianValueType,jacobianProperties...>;
803 using FunctorType = FunctorCellTools::F_setJacobian<JacobianViewType,WorksetType,BasisGradientsType> ;
804
805 // resolve the -1 default argument for endCell into the true end cell index
806 int endCellResolved = (endCell == -1) ? worksetCell.extent_int(0) : endCell;
807
808 using range_policy_type = Kokkos::Experimental::MDRangePolicy
809 < ExecSpaceType, Kokkos::Experimental::Rank<2>, Kokkos::IndexType<ordinal_type> >;
810 range_policy_type policy( { 0, 0 },
811 { jacobian.extent(0), jacobian.extent(1) } );
812 Kokkos::parallel_for( policy, FunctorType(jacobian, worksetCell, gradients, startCell, endCellResolved) );
813 }
814
815 template<typename DeviceType>
816 template<typename jacobianValueType, class ...jacobianProperties,
817 typename pointValueType, class ...pointProperties,
818 typename WorksetType,
819 typename HGradBasisType>
820 void
821 CellTools<DeviceType>::
822 setJacobian( Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian,
823 const Kokkos::DynRankView<pointValueType,pointProperties...> points,
824 const WorksetType worksetCell,
825 const Teuchos::RCP<HGradBasisType> basis,
826 const int startCell, const int endCell) {
827 constexpr bool are_accessible =
828 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
829 typename decltype(jacobian)::memory_space>::accessible &&
830 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
831 typename decltype(points)::memory_space>::accessible;
832 static_assert(are_accessible, "CellTools<DeviceType>::setJacobian(..): input/output views' memory spaces are not compatible with DeviceType");
833
834#ifdef HAVE_INTREPID2_DEBUG
835 CellTools_setJacobianArgs(jacobian, points, worksetCell, basis->getBaseCellTopology(), startCell, endCell);
836 //static_assert(std::is_same( pointValueType, decltype(basis->getDummyOutputValue()) ));
837#endif
838 const auto cellTopo = basis->getBaseCellTopology();
839 const ordinal_type spaceDim = cellTopo.getDimension();
840 const ordinal_type numCells = jacobian.extent(0);
841
842 //points can be rank-2 (P,D), or rank-3 (C,P,D)
843 const ordinal_type pointRank = points.rank();
844 const ordinal_type numPoints = (pointRank == 2 ? points.extent(0) : points.extent(1));
845 const ordinal_type basisCardinality = basis->getCardinality();
846
847 // the following does not work for RCP; its * operator returns reference not the object
848 //typedef typename decltype(*basis)::output_value_type gradValueType;
849 //typedef Kokkos::DynRankView<decltype(basis->getDummyOutputValue()),DeviceType> gradViewType;
850
851 auto vcprop = Kokkos::common_view_alloc_prop(points);
852 using GradViewType = Kokkos::DynRankView<typename decltype(vcprop)::value_type,DeviceType>;
853
854 GradViewType grads;
855
856 switch (pointRank) {
857 case 2: {
858 // For most FEMs
859 grads = GradViewType(Kokkos::view_alloc("CellTools::setJacobian::grads", vcprop),basisCardinality, numPoints, spaceDim);
860 basis->getValues(grads,
861 points,
862 OPERATOR_GRAD);
863 break;
864 }
865 case 3: {
866 // For CVFEM
867 grads = GradViewType(Kokkos::view_alloc("CellTools::setJacobian::grads", vcprop), numCells, basisCardinality, numPoints, spaceDim);
868 for (ordinal_type cell=0;cell<numCells;++cell)
869 basis->getValues(Kokkos::subview( grads, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() ),
870 Kokkos::subview( points, cell, Kokkos::ALL(), Kokkos::ALL() ),
871 OPERATOR_GRAD);
872 break;
873 }
874 }
875
876 setJacobian(jacobian, worksetCell, grads, startCell, endCell);
877 }
878
879 template<typename DeviceType>
880 template<typename jacobianInvValueType, class ...jacobianInvProperties,
881 typename jacobianValueType, class ...jacobianProperties>
882 void
883 CellTools<DeviceType>::
884 setJacobianInv( Kokkos::DynRankView<jacobianInvValueType,jacobianInvProperties...> jacobianInv,
885 const Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian ) {
886#ifdef HAVE_INTREPID2_DEBUG
887 CellTools_setJacobianInvArgs(jacobianInv, jacobian);
888#endif
889 RealSpaceTools<DeviceType>::inverse(jacobianInv, jacobian);
890 }
891
892 template<typename DeviceType>
893 template<typename jacobianDetValueType, class ...jacobianDetProperties,
894 typename jacobianValueType, class ...jacobianProperties>
895 void
896 CellTools<DeviceType>::
897 setJacobianDet( Kokkos::DynRankView<jacobianDetValueType,jacobianDetProperties...> jacobianDet,
898 const Kokkos::DynRankView<jacobianValueType,jacobianProperties...> jacobian ) {
899#ifdef HAVE_INTREPID2_DEBUG
900 CellTools_setJacobianDetArgs(jacobianDet, jacobian);
901#endif
902 RealSpaceTools<DeviceType>::det(jacobianDet, jacobian);
903 }
904
905}
906
907#endif
Intrepid2::CellTools_setJacobianInvArgs
void CellTools_setJacobianInvArgs(const jacobianInvViewType jacobianInv, const jacobianViewType jacobian)
Validates arguments to Intrepid2::CellTools::setJacobianInv.
Definition: Intrepid2_CellToolsDefValidateArguments.hpp:155
Intrepid2::CellTools_setJacobianDetArgs
void CellTools_setJacobianDetArgs(const jacobianDetViewType jacobianDet, const jacobianViewType jacobian)
Validates arguments to Intrepid2::CellTools::setJacobianDet.
Definition: Intrepid2_CellToolsDefValidateArguments.hpp:187
Intrepid2::CONSTANT
@ CONSTANT
does not vary
Definition: Intrepid2_Data.hpp:33
Intrepid2::BLOCK_PLUS_DIAGONAL
@ BLOCK_PLUS_DIAGONAL
one of two dimensions in a matrix; bottom-right part of matrix is diagonal
Definition: Intrepid2_Data.hpp:35
Intrepid2_Kernels.hpp
Header file for small functions used in Intrepid2.
Intrepid2::getMatchingViewWithLabel
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.
Definition: Intrepid2_Utils.hpp:361
Intrepid2::CellTools::setJacobianDet
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.
Definition: Intrepid2_CellToolsDefJacobian.hpp:897
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:124
Intrepid2::CellTools::setJacobianInv
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.
Definition: Intrepid2_CellToolsDefJacobian.hpp:884
Intrepid2::CellTools::setJacobian
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.
Definition: Intrepid2_CellToolsDefJacobian.hpp:822
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:158
Intrepid2::Data
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
Definition: Intrepid2_Data.hpp:199
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:1764
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:1498
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:1512
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:1505
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:1365
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:1724
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:1491
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:1372
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:1605
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::RealSpaceTools::det
static void det(DeterminantArrayViewType detArray, const MatrixViewType inMats)
Computes determinants of matrices stored in an array of total rank 3 (array of matrices),...
Definition: Intrepid2_RealSpaceToolsDef.hpp:856
Intrepid2::FunctorCellTools::F_setJacobian
Functor for calculation of Jacobian on cell workset see Intrepid2::CellTools for more.
Definition: Intrepid2_CellToolsDefJacobian.hpp:75
Generated by doxygen 1.9.4