Intrepid2
Intrepid2_CubatureTensorDef.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_CUBATURE_TENSOR_DEF_HPP__
17#define __INTREPID2_CUBATURE_TENSOR_DEF_HPP__
18
19namespace Intrepid2 {
20
21 template<typename DT, typename PT, typename WT>
22 template<typename cubPointValueType, class ...cubPointProperties,
23 typename cubWeightValueType, class ...cubWeightProperties>
24 void
25 CubatureTensor<DT,PT,WT>::
26 getCubatureImpl( Kokkos::DynRankView<cubPointValueType, cubPointProperties...> cubPoints,
27 Kokkos::DynRankView<cubWeightValueType,cubWeightProperties...> cubWeights ) const {
28#ifdef HAVE_INTREPID2_DEBUG
29 // check size of cubPoints and cubWeights
30 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(cubPoints.extent(0)) < this->getNumPoints() ||
31 static_cast<ordinal_type>(cubPoints.extent(1)) < this->getDimension() ||
32 static_cast<ordinal_type>(cubWeights.extent(0)) < this->getNumPoints(), std::out_of_range,
33 ">>> ERROR (CubatureTensor): Insufficient space allocated for cubature points or weights.");
34#endif
35 using cubPointViewType = Kokkos::DynRankView<cubPointValueType, DT>;
36 using cubWeightViewType = Kokkos::DynRankView<cubWeightValueType,DT>;
37
38 // mirroring and where the data is problematic... when it becomes a problem, then deal with it.
39 cubPointViewType tmpPoints [Parameters::MaxTensorComponents];
40 cubWeightViewType tmpWeights[Parameters::MaxTensorComponents];
41
42 // this temporary allocation can be member of cubature; for now, let's do this way.
43 // this is cubature setup on the reference cell and called for tensor elements.
44 for (auto k=0;k<this->numCubatures_;++k) {
45 const auto &cub = this->cubatures_[k];
46 tmpPoints [k] = cubPointViewType ("CubatureTensor::getCubature::tmpPoints", cub.getNumPoints(), cub.getDimension());
47 tmpWeights[k] = cubWeightViewType("CubatureTensor::getCubature::tmpWeights", cub.getNumPoints());
48 cub.getCubature(tmpPoints[k], tmpWeights[k]);
49 }
50
51 // when the input containers are device space, this is better computed on host and copy to devices
52 // fill tensor cubature
53 {
54 ordinal_type offset[Parameters::MaxTensorComponents+1] = {};
55 for (auto k=0;k<this->numCubatures_;++k) {
56 offset[k+1] = offset[k] + this->cubatures_[k].getDimension();
57 }
58 ordinal_type ii = 0, i[3] = {};
59
60 cubPointViewType cubPoints_device("cubPoints_device", cubPoints.extent(0), cubPoints.extent(1));
61 cubWeightViewType cubWeights_device("cubWeights_device", cubWeights.extent(0));
62
63 auto cubPoints_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cubPoints_device);
64 auto cubWeights_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cubWeights_device);
65
66 Kokkos::deep_copy(cubPoints_host, 0.0);
67 Kokkos::deep_copy(cubWeights_host, 1.0);
68
69 switch (this->numCubatures_) {
70 case 2: {
71 const ordinal_type npts[] = { this->cubatures_[0].getNumPoints(), this->cubatures_[1].getNumPoints() };
72 const ordinal_type dim [] = { this->cubatures_[0].getDimension(), this->cubatures_[1].getDimension() };
73
75 typename cubWeightViewType::host_mirror_type tmpWeights_host[2];
76 tmpWeights_host[0] = Kokkos::create_mirror_view(tmpWeights[0]);
77 tmpWeights_host[1] = Kokkos::create_mirror_view(tmpWeights[1]);
78 Kokkos::deep_copy(tmpWeights_host[0], tmpWeights[0]);
79 Kokkos::deep_copy(tmpWeights_host[1], tmpWeights[1]);
80
81 typename cubPointViewType::host_mirror_type tmpPoints_host[2];
82 tmpPoints_host[0] = Kokkos::create_mirror_view(tmpPoints[0]);
83 tmpPoints_host[1] = Kokkos::create_mirror_view(tmpPoints[1]);
84
85 Kokkos::deep_copy(tmpPoints_host[0], tmpPoints[0]);
86 Kokkos::deep_copy(tmpPoints_host[1], tmpPoints[1]);
87
88 for (i[1]=0;i[1]<npts[1];++i[1])
89 for (i[0]=0;i[0]<npts[0];++i[0]) {
90 for (auto nc=0;nc<2;++nc) {
91 cubWeights_host(ii) *= tmpWeights_host[nc](i[nc]);
92 for (ordinal_type j=0;j<dim[nc];++j)
93 cubPoints_host(ii, offset[nc]+j) = tmpPoints_host[nc](i[nc], j);
94 }
95 ++ii;
96 }
97 break;
98 }
99 case 3: {
100 const ordinal_type npts[] = { this->cubatures_[0].getNumPoints(), this->cubatures_[1].getNumPoints(), this->cubatures_[2].getNumPoints() };
101 const ordinal_type dim [] = { this->cubatures_[0].getDimension(), this->cubatures_[1].getDimension(), this->cubatures_[2].getDimension() };
102
104 typename cubWeightViewType::host_mirror_type tmpWeights_host[3];
105 tmpWeights_host[0] = Kokkos::create_mirror_view(tmpWeights[0]);
106 tmpWeights_host[1] = Kokkos::create_mirror_view(tmpWeights[1]);
107 tmpWeights_host[2] = Kokkos::create_mirror_view(tmpWeights[2]);
108
109 Kokkos::deep_copy(tmpWeights_host[0], tmpWeights[0]);
110 Kokkos::deep_copy(tmpWeights_host[1], tmpWeights[1]);
111 Kokkos::deep_copy(tmpWeights_host[2], tmpWeights[2]);
112
113 typename cubPointViewType::host_mirror_type tmpPoints_host[3];
114 tmpPoints_host[0] = Kokkos::create_mirror_view(tmpPoints[0]);
115 tmpPoints_host[1] = Kokkos::create_mirror_view(tmpPoints[1]);
116 tmpPoints_host[2] = Kokkos::create_mirror_view(tmpPoints[2]);
117
118 Kokkos::deep_copy(tmpPoints_host[0], tmpPoints[0]);
119 Kokkos::deep_copy(tmpPoints_host[1], tmpPoints[1]);
120 Kokkos::deep_copy(tmpPoints_host[2], tmpPoints[2]);
121
122 for (i[2]=0;i[2]<npts[2];++i[2])
123 for (i[1]=0;i[1]<npts[1];++i[1])
124 for (i[0]=0;i[0]<npts[0];++i[0]) {
125 for (auto nc=0;nc<3;++nc) {
126 cubWeights_host(ii) *= tmpWeights_host[nc](i[nc]);
127 for (ordinal_type j=0;j<dim[nc];++j)
128 cubPoints_host(ii, offset[nc]+j) = tmpPoints_host[nc](i[nc], j);
129 }
130 ++ii;
131 }
132 break;
133 }
134 default: {
135 INTREPID2_TEST_FOR_EXCEPTION( !(this->numCubatures_ == 2 || this->numCubatures_ == 3), std::runtime_error,
136 ">>> ERROR (CubatureTensor::getCubature): CubatureTensor supports only 2 or 3 component direct cubatures.");
137 }
138 }
139 Kokkos::deep_copy(cubPoints_device, cubPoints_host);
140 Kokkos::deep_copy(cubWeights_device, cubWeights_host);
141
142 Kokkos::deep_copy(cubPoints, cubPoints_device);
143 Kokkos::deep_copy(cubWeights, cubWeights_device);
144 }
145 }
146
147} // end namespace Intrepid2
148
149#endif
Intrepid2::CubatureTensor::getCubatureImpl
void getCubatureImpl(Kokkos::DynRankView< cubPointValueType, cubPointProperties... > cubPoints, Kokkos::DynRankView< cubWeightValueType, cubWeightProperties... > cubWeights) const
Definition Intrepid2_CubatureTensorDef.hpp:26
Intrepid2::Parameters::MaxTensorComponents
static constexpr ordinal_type MaxTensorComponents
Maximum number of tensor/Cartesian products that can be taken: this allows hypercube basis in 7D to b...
Definition Intrepid2_Types.hpp:120
Generated by doxygen 1.9.8