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