Intrepid2
Intrepid2_ArrayToolsDefContractions.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid 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
49#ifndef __INTREPID2_ARRAYTOOLS_DEF_CONTRACTIONS_HPP__
50#define __INTREPID2_ARRAYTOOLS_DEF_CONTRACTIONS_HPP__
51
52namespace Intrepid2 {
53
54
55 namespace FunctorArrayTools {
59 template < typename outFieldViewType , typename leftFieldViewType , typename rightFieldViewType >
61 outFieldViewType _outputFields;
62 leftFieldViewType _leftFields;
63 rightFieldViewType _rightFields;
64 const bool _sumInto;
65 typedef typename outFieldViewType::value_type value_type;
66
67 KOKKOS_INLINE_FUNCTION
68 F_contractFieldField(outFieldViewType outputFields_,
69 leftFieldViewType leftFields_,
70 rightFieldViewType rightFields_,
71 const bool sumInto_)
72 : _outputFields(outputFields_), _leftFields(leftFields_), _rightFields(rightFields_), _sumInto(sumInto_) {}
73
74 KOKKOS_INLINE_FUNCTION
75 void operator()(const size_type iter) const {
76 size_type cl, lbf, rbf;
77 unrollIndex( cl, lbf, rbf,
78 _outputFields.extent(0),
79 _outputFields.extent(1),
80 _outputFields.extent(2),
81 iter );
82
83 const size_type npts = _leftFields.extent(2);
84 const ordinal_type iend = _leftFields.extent(3);
85 const ordinal_type jend = _leftFields.extent(4);
86
87 value_type tmp(0);
88 for (size_type qp = 0; qp < npts; ++qp)
89 for (ordinal_type i = 0; i < iend; ++i)
90 for (ordinal_type j = 0; j < jend; ++j)
91 tmp += _leftFields(cl, lbf, qp, i, j)*_rightFields(cl, rbf, qp, i, j);
92 if (_sumInto)
93 _outputFields( cl, lbf, rbf ) = _outputFields( cl, lbf, rbf ) + tmp;
94 else
95 _outputFields( cl, lbf, rbf ) = tmp;
96 }
97 };
98 } //end namespace
99
100 template<typename DeviceType>
101 template<typename outputFieldValueType, class ...outputFieldProperties,
102 typename leftFieldValueType, class ...leftFieldProperties,
103 typename rightFieldValueType, class ...rightFieldProperties>
104 void
106 contractFieldField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
107 const Kokkos::DynRankView<leftFieldValueType, leftFieldProperties...> leftFields,
108 const Kokkos::DynRankView<rightFieldValueType, rightFieldProperties...> rightFields,
109 const bool sumInto ) {
110
111 typedef Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outFieldViewType;
112 typedef Kokkos::DynRankView<leftFieldValueType,leftFieldProperties...> leftFieldViewType;
113 typedef Kokkos::DynRankView<rightFieldValueType,rightFieldProperties...> rightFieldViewType;
115
116 const size_type loopSize = leftFields.extent(0)*leftFields.extent(1)*rightFields.extent(1);
117 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
118 Kokkos::parallel_for( policy, FunctorType(outputFields, leftFields, rightFields, sumInto) );
119 }
120
121
122 namespace FunctorArrayTools {
126 template < typename outputFieldsViewType , typename inputDataViewType , typename inputFieldsViewType >
128 outputFieldsViewType _outputFields;
129 inputDataViewType _inputData;
130 inputFieldsViewType _inputFields;
131 const bool _sumInto;
132 typedef typename outputFieldsViewType::value_type value_type;
133
134 KOKKOS_INLINE_FUNCTION
135 F_contractDataField(outputFieldsViewType outputFields_,
136 inputDataViewType inputData_,
137 inputFieldsViewType inputFields_,
138 const bool sumInto_)
139 : _outputFields(outputFields_), _inputData(inputData_), _inputFields(inputFields_), _sumInto(sumInto_) {}
140
141 KOKKOS_DEFAULTED_FUNCTION
142 ~F_contractDataField() = default;
143
144 KOKKOS_INLINE_FUNCTION
145 void operator()(const size_type iter) const {
146 size_type cl, bf;
147 unrollIndex( cl, bf,
148 _inputFields.extent(0),
149 _inputFields.extent(1),
150 iter );
151
152 auto result = Kokkos::subview( _outputFields, cl, bf );
153
154 const auto field = Kokkos::subview( _inputFields, cl, bf, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
155 const auto data = Kokkos::subview( _inputData, cl, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
156
157 const size_type npts = field.extent(0);
158 const ordinal_type iend = field.extent(1);
159 const ordinal_type jend = field.extent(2);
160
161 value_type tmp(0);
162
163 if(_inputData.extent(1) != 1)
164 for (size_type qp = 0; qp < npts; ++qp)
165 for (ordinal_type i = 0; i < iend; ++i)
166 for (ordinal_type j = 0; j < jend; ++j)
167 tmp += field(qp, i, j) * data(qp, i, j);
168 else
169 for (size_type qp = 0; qp < npts; ++qp)
170 for (ordinal_type i = 0; i < iend; ++i)
171 for (ordinal_type j = 0; j < jend; ++j)
172 tmp += field(qp, i, j) * data(0, i, j);
173
174 if (_sumInto)
175 result() = result() + tmp;
176 else
177 result() = tmp;
178 }
179 };
180 } //namespace
181
182 template<typename DeviceType>
183 template<typename outputFieldValueType, class ...outputFieldProperties,
184 typename inputDataValueType, class ...inputDataProperties,
185 typename inputFieldValueType, class ...inputFieldProperties>
186 void
188 contractDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
189 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
190 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
191 const bool sumInto ) {
192
193 typedef Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFieldsViewType;
194 typedef Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputDataViewType;
195 typedef Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFieldsViewType;
197
198 const size_type loopSize = inputFields.extent(0)*inputFields.extent(1);
199 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
200 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields, sumInto) );
201 }
202
203
204 namespace FunctorArrayTools {
208 template < typename outputDataViewType , typename inputDataLeftViewType , typename inputDataRightViewType >
210 outputDataViewType _outputData;
211 inputDataLeftViewType _inputDataLeft;
212 inputDataRightViewType _inputDataRight;
213 const bool _sumInto;
214 typedef typename outputDataViewType::value_type value_type;
215
216 KOKKOS_INLINE_FUNCTION
217 F_contractDataData(outputDataViewType outputData_,
218 inputDataLeftViewType inputDataLeft_,
219 inputDataRightViewType inputDataRight_,
220 const bool sumInto_)
221 : _outputData(outputData_), _inputDataLeft(inputDataLeft_), _inputDataRight(inputDataRight_), _sumInto(sumInto_) {}
222
223 KOKKOS_DEFAULTED_FUNCTION
224 ~F_contractDataData() = default;
225
226 KOKKOS_INLINE_FUNCTION
227 void operator()(const size_type iter) const {
228 const size_type cl = iter;
229
230 auto result = Kokkos::subview( _outputData, cl );
231 const auto left = Kokkos::subview( _inputDataLeft, cl, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
232 const auto right = Kokkos::subview( _inputDataRight, cl, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
233
234 size_type npts = left.extent(0);
235 ordinal_type iend = left.extent(1);
236 ordinal_type jend = left.extent(2);
237
238 value_type tmp(0);
239 for (size_type qp = 0; qp < npts; ++qp)
240 for (ordinal_type i = 0; i < iend; ++i)
241 for (ordinal_type j = 0; j < jend; ++j)
242 tmp += left(qp, i, j)*right(qp, i, j);
243
244 if (_sumInto)
245 result() = result() + tmp;
246 else
247 result() = tmp;
248 }
249 };
250 } //namespace
251
252 template<typename DeviceType>
253 template<typename outputDataValueType, class ...outputDataProperties,
254 typename inputDataLeftValueType, class ...inputDataLeftProperties,
255 typename inputDataRightValueType, class ...inputDataRightProperties>
256 void
258 contractDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
259 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
260 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
261 const bool sumInto ) {
262 typedef Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputDataViewType;
263 typedef Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeftViewType;
264 typedef Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRightViewType;
266
267 const size_type loopSize = inputDataLeft.extent(0);
268 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
269 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight, sumInto) );
270 }
271
272
273
274 template<typename DeviceType>
275 template<typename outputFieldValueType, class ...outputFieldProperties,
276 typename leftFieldValueType, class ...leftFieldProperties,
277 typename rightFieldValueType, class ...rightFieldProperties>
278 void
280 contractFieldFieldScalar( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
281 const Kokkos::DynRankView<leftFieldValueType, leftFieldProperties...> leftFields,
282 const Kokkos::DynRankView<rightFieldValueType, rightFieldProperties...> rightFields,
283 const bool sumInto ) {
284
285#ifdef HAVE_INTREPID2_DEBUG
286 {
287 INTREPID2_TEST_FOR_EXCEPTION( leftFields.rank() != 3, std::invalid_argument,
288 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of the left input argument must equal 3!");
289 INTREPID2_TEST_FOR_EXCEPTION( rightFields.rank() != 3, std::invalid_argument,
290 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of right input argument must equal 3!");
291 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
292 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of output argument must equal 3!");
293 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(0) != rightFields.extent(0), std::invalid_argument,
294 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
295 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(2) != rightFields.extent(2), std::invalid_argument,
296 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
297 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != rightFields.extent(0), std::invalid_argument,
298 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
299 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(1) != leftFields.extent(1), std::invalid_argument,
300 ">>> ERROR (ArrayTools::contractFieldFieldScalar): First dimension of output container and first dimension of left input container must agree!");
301 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != rightFields.extent(1), std::invalid_argument,
302 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimension of output container and first dimension of right input container must agree!");
303
304 }
305#endif
306
308 leftFields,
309 rightFields,
310 sumInto );
311 }
312
313
314 template<typename DeviceType>
315 template<typename outputFieldValueType, class ...outputFieldProperties,
316 typename leftFieldValueType, class ...leftFieldProperties,
317 typename rightFieldValueType, class ...rightFieldProperties>
318 void
320 contractFieldFieldVector( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
321 const Kokkos::DynRankView<leftFieldValueType, leftFieldProperties...> leftFields,
322 const Kokkos::DynRankView<rightFieldValueType, rightFieldProperties...> rightFields,
323 const bool sumInto ) {
324
325#ifdef HAVE_INTREPID2_DEBUG
326 {
327 INTREPID2_TEST_FOR_EXCEPTION( leftFields.rank() != 4, std::invalid_argument,
328 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of the left input argument must equal 4!");
329 INTREPID2_TEST_FOR_EXCEPTION( rightFields.rank() != 4, std::invalid_argument,
330 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of right input argument must equal 4!");
331 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
332 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of output argument must equal 3!");
333 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(0) != rightFields.extent(0), std::invalid_argument,
334 ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
335 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(2) != rightFields.extent(2), std::invalid_argument,
336 ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
337 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(3) != rightFields.extent(3), std::invalid_argument,
338 ">>> ERROR (ArrayTools::contractFieldFieldVector): Third dimensions (numbers of vector components) of the left and right input containers must agree!");
339 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != rightFields.extent(0), std::invalid_argument,
340 ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
341 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(1) != leftFields.extent(1), std::invalid_argument,
342 ">>> ERROR (ArrayTools::contractFieldFieldVector): First dimension of output container and first dimension of left input container must agree!");
343 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != rightFields.extent(1), std::invalid_argument,
344 ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimension of output container and first dimension of right input container must agree!");
345 }
346#endif
347
349 leftFields,
350 rightFields,
351 sumInto );
352 }
353
354
355 template<typename DeviceType>
356 template<typename outputFieldValueType, class ...outputFieldProperties,
357 typename leftFieldValueType, class ...leftFieldProperties,
358 typename rightFieldValueType, class ...rightFieldProperties>
359 void
361 contractFieldFieldTensor( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
362 const Kokkos::DynRankView<leftFieldValueType, leftFieldProperties...> leftFields,
363 const Kokkos::DynRankView<rightFieldValueType, rightFieldProperties...> rightFields,
364 const bool sumInto ) {
365
366#ifdef HAVE_INTREPID2_DEBUG
367 {
368 INTREPID2_TEST_FOR_EXCEPTION( leftFields.rank() != 5, std::invalid_argument,
369 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of the left input argument must equal 5!");
370 INTREPID2_TEST_FOR_EXCEPTION( rightFields.rank() != 5, std::invalid_argument,
371 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of right input argument must equal 5!");
372 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
373 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of output argument must equal 3!");
374 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(0) != rightFields.extent(0), std::invalid_argument,
375 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
376 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(2) != rightFields.extent(2), std::invalid_argument,
377 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimensions (numbers of integration points) of the left and right input containers must agree!");
378 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(3) != rightFields.extent(3), std::invalid_argument,
379 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Third dimensions (first tensor dimensions) of the left and right input containers must agree!");
380 INTREPID2_TEST_FOR_EXCEPTION( leftFields.extent(4) != rightFields.extent(4), std::invalid_argument,
381 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Fourth dimensions (second tensor dimensions) of the left and right input containers must agree!");
382 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != rightFields.extent(0), std::invalid_argument,
383 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
384 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(1) != leftFields.extent(1), std::invalid_argument,
385 ">>> ERROR (ArrayTools::contractFieldFieldTensor): First dimension of output container and first dimension of left input container must agree!");
386 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != rightFields.extent(1), std::invalid_argument,
387 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimension of output container and first dimension of right input container must agree!");
388 }
389#endif
390
392 leftFields,
393 rightFields,
394 sumInto );
395 }
396
397
398 template<typename DeviceType>
399 template<typename outputFieldValueType, class ...outputFieldProperties,
400 typename inputDataValueType, class ...inputDataProperties,
401 typename inputFieldValueType, class ...inputFieldProperties>
402 void
404 contractDataFieldScalar( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
405 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
406 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
407 const bool sumInto ) {
408
409#ifdef HAVE_INTREPID2_DEBUG
410 {
411 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != 3, std::invalid_argument,
412 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the fields input argument must equal 3!");
413 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 2, std::invalid_argument,
414 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the data input argument must equal 2!");
415 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 2, std::invalid_argument,
416 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of output argument must equal 2!");
417 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != inputData.extent(0), std::invalid_argument,
418 ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
419
420 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2) &&
421 inputData.extent(1) != 1, std::invalid_argument,
422 ">>> ERROR (ArrayTools::contractDataFieldScalar): Second dimension of fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
423 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputFields.extent(0), std::invalid_argument,
424 ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
425 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(1) != inputFields.extent(1), std::invalid_argument,
426 ">>> ERROR (ArrayTools::contractDataFieldScalar): First dimensions (number of fields) of the fields input and output containers must agree!");
427 }
428#endif
429
431 inputData,
432 inputFields,
433 sumInto );
434 }
435
436
437 template<typename DeviceType>
438 template<typename outputFieldValueType, class ...outputFieldProperties,
439 typename inputDataValueType, class ...inputDataProperties,
440 typename inputFieldValueType, class ...inputFieldProperties>
441 void
443 contractDataFieldVector( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
444 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
445 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
446 const bool sumInto ) {
447
448#ifdef HAVE_INTREPID2_DEBUG
449 {
450 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != 4, std::invalid_argument,
451 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the fields input argument must equal 4!");
452 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
453 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the data input argument must equal 3!");
454 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 2, std::invalid_argument,
455 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of output argument must equal 2!");
456 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != inputData.extent(0), std::invalid_argument,
457 ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
458 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2) &&
459 inputData.extent(1) != 1, std::invalid_argument,
460 ">>> ERROR (ArrayTools::contractDataFieldVector): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
461 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(3) != inputData.extent(2), std::invalid_argument,
462 ">>> ERROR (ArrayTools::contractDataFieldVector): Third dimension of the fields input container and second dimension of data input container (vector index) must agree!");
463 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputFields.extent(0), std::invalid_argument,
464 ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
465 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(1) != inputFields.extent(1), std::invalid_argument,
466 ">>> ERROR (ArrayTools::contractDataFieldVector): First dimensions of output container and fields input container (number of fields) must agree!");
467 }
468#endif
469
471 inputData,
472 inputFields,
473 sumInto );
474 }
475
476
477
478 template<typename DeviceType>
479 template<typename outputFieldValueType, class ...outputFieldProperties,
480 typename inputDataValueType, class ...inputDataProperties,
481 typename inputFieldValueType, class ...inputFieldProperties>
482 void
484 contractDataFieldTensor( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
485 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
486 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
487 const bool sumInto ) {
488
489#ifdef HAVE_INTREPID2_DEBUG
490 {
491 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != 5, std::invalid_argument,
492 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the fields input argument must equal 5!");
493 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 4, std::invalid_argument,
494 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the data input argument must equal 4!");
495 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 2, std::invalid_argument,
496 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of output argument must equal 2!");
497 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != inputData.extent(0), std::invalid_argument,
498 ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
499 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2) && inputData.extent(1) != 1, std::invalid_argument,
500 ">>> ERROR (ArrayTools::contractDataFieldTensor): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
501 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(3) != inputData.extent(2), std::invalid_argument,
502 ">>> ERROR (ArrayTools::contractDataFieldTensor): Third dimension of the fields input container and second dimension of data input container (first tensor dimension) must agree!");
503 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(4) != inputData.extent(3), std::invalid_argument,
504 ">>> ERROR (ArrayTools::contractDataFieldTensor): Fourth dimension of the fields input container and third dimension of data input container (second tensor dimension) must agree!");
505 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputFields.extent(0), std::invalid_argument,
506 ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!");
507 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(1) != inputFields.extent(1), std::invalid_argument,
508 ">>> ERROR (ArrayTools::contractDataFieldTensor): First dimensions (number of fields) of output container and fields input container must agree!");
509 }
510#endif
511
513 inputData,
514 inputFields,
515 sumInto );
516 }
517
518
519
520 template<typename DeviceType>
521 template<typename outputDataValueType, class ...outputDataProperties,
522 typename inputDataLeftValueType, class ...inputDataLeftProperties,
523 typename inputDataRightValueType, class ...inputDataRightProperties>
524 void
526 contractDataDataScalar( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
527 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
528 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
529 const bool sumInto ) {
530
531#ifdef HAVE_INTREPID2_DEBUG
532 {
533 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 2, std::invalid_argument,
534 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of the left input argument must equal 2!");
535 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != 2, std::invalid_argument,
536 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of right input argument must equal 2!");
537 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 1, std::invalid_argument,
538 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of output argument must equal 1!");
539 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
540 ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
541 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
542 ">>> ERROR (ArrayTools::contractDataDataScalar): First dimensions (numbers of integration points) of the left and right input containers must agree!");
543 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataRight.extent(0), std::invalid_argument,
544 ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
545 }
546#endif
547
549 inputDataLeft,
550 inputDataRight,
551 sumInto );
552 }
553
554
555 template<typename DeviceType>
556 template<typename outputDataValueType, class ...outputDataProperties,
557 typename inputDataLeftValueType, class ...inputDataLeftProperties,
558 typename inputDataRightValueType, class ...inputDataRightProperties>
559 void
561 contractDataDataVector( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
562 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
563 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
564 const bool sumInto ) {
565
566#ifdef HAVE_INTREPID2_DEBUG
567 {
568 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
569 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of the left input argument must equal 3!");
570 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != 3, std::invalid_argument,
571 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of right input argument must equal 3!");
572 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 1, std::invalid_argument,
573 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of output argument must equal 1!");
574 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
575 ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
576 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
577 ">>> ERROR (ArrayTools::contractDataDataVector): First dimensions (numbers of integration points) of the left and right input containers must agree!");
578 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) != inputDataRight.extent(2), std::invalid_argument,
579 ">>> ERROR (ArrayTools::contractDataDataVector): Second dimensions (numbers of vector components) of the left and right input containers must agree!");
580 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataRight.extent(0), std::invalid_argument,
581 ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
582 }
583#endif
584
586 inputDataLeft,
587 inputDataRight,
588 sumInto );
589 }
590
591
592 template<typename DeviceType>
593 template<typename outputDataValueType, class ...outputDataProperties,
594 typename inputDataLeftValueType, class ...inputDataLeftProperties,
595 typename inputDataRightValueType, class ...inputDataRightProperties>
596 void
598 contractDataDataTensor( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
599 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
600 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
601 const bool sumInto ) {
602
603#ifdef HAVE_INTREPID2_DEBUG
604 {
605 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 4, std::invalid_argument,
606 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of the left input argument must equal 4");
607 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != 4, std::invalid_argument,
608 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of right input argument must equal 4!");
609 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 1, std::invalid_argument,
610 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of output argument must equal 1!");
611 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
612 ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!");
613 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
614 ">>> ERROR (ArrayTools::contractDataDataTensor): First dimensions (numbers of integration points) of the left and right input containers must agree!");
615 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) != inputDataRight.extent(2), std::invalid_argument,
616 ">>> ERROR (ArrayTools::contractDataDataTensor): Second dimensions (first tensor dimensions) of the left and right input containers must agree!");
617 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) != inputDataRight.extent(3), std::invalid_argument,
618 ">>> ERROR (ArrayTools::contractDataDataTensor): Third dimensions (second tensor dimensions) of the left and right input containers must agree!");
619 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataRight.extent(0), std::invalid_argument,
620 ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!");
621 }
622#endif
623
625 inputDataLeft,
626 inputDataRight,
627 sumInto );
628 }
629
630}
631#endif
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
static void contractFieldFieldTensor(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< leftFieldValueType, leftFieldProperties... > leftFields, const Kokkos::DynRankView< rightFieldValueType, rightFieldProperties... > rightFields, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1, and D2 of two rank-5 containers with dimensions (...
static void contractDataDataVector(Kokkos::DynRankView< outputDataValueType, outputDataProperties... > outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties... > inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties... > inputDataRight, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D of rank-3 containers with dimensions (C,...
static void contractFieldFieldScalar(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< leftFieldValueType, leftFieldProperties... > leftFields, const Kokkos::DynRankView< rightFieldValueType, rightFieldProperties... > rightFields, const bool sumInto=false)
Contracts the "point" dimension P of two rank-3 containers with dimensions (C,L,P) and (C,...
static void contractDataFieldScalar(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties... > inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties... > inputFields, const bool sumInto=false)
Contracts the "point" dimensions P of a rank-3 containers and a rank-2 container with dimensions (C,...
static void contractDataDataTensor(Kokkos::DynRankView< outputDataValueType, outputDataProperties... > outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties... > inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties... > inputDataRight, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1 and D2 of rank-4 containers with dimensions (C,...
static void contractDataFieldVector(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties... > inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties... > inputFields, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D of a rank-4 container and a rank-3 container wit...
static void contractDataFieldTensor(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties... > inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties... > inputFields, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1 and D2 of a rank-5 container and a rank-4 containe...
static void contractFieldFieldVector(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< leftFieldValueType, leftFieldProperties... > leftFields, const Kokkos::DynRankView< rightFieldValueType, rightFieldProperties... > rightFields, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D1 of two rank-4 containers with dimensions (C,...
static void contractDataDataScalar(Kokkos::DynRankView< outputDataValueType, outputDataProperties... > outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties... > inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties... > inputDataRight, const bool sumInto=false)
Contracts the "point" dimensions P of rank-2 containers with dimensions (C,P), and returns the result...
Functor to contractDataData see Intrepid2::ArrayTools for more.
Functor to contractDataField see Intrepid2::ArrayTools for more.
Functor to contractFieldField see Intrepid2::ArrayTools for more.