49 #ifndef __INTREPID2_CUBATURE_TENSOR_DEF_HPP__
50 #define __INTREPID2_CUBATURE_TENSOR_DEF_HPP__
54 template<
typename DT,
typename PT,
typename WT>
55 template<
typename cubPointValueType,
class ...cubPointProperties,
56 typename cubWeightValueType,
class ...cubWeightProperties>
59 getCubatureImpl( Kokkos::DynRankView<cubPointValueType, cubPointProperties...> cubPoints,
60 Kokkos::DynRankView<cubWeightValueType,cubWeightProperties...> cubWeights )
const {
61 #ifdef HAVE_INTREPID2_DEBUG
63 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(cubPoints.extent(0)) < this->getNumPoints() ||
64 static_cast<ordinal_type
>(cubPoints.extent(1)) < this->getDimension() ||
65 static_cast<ordinal_type
>(cubWeights.extent(0)) < this->getNumPoints(), std::out_of_range,
66 ">>> ERROR (CubatureTensor): Insufficient space allocated for cubature points or weights.");
68 using cubPointViewType = Kokkos::DynRankView<cubPointValueType, DT>;
69 using cubWeightViewType = Kokkos::DynRankView<cubWeightValueType,DT>;
77 for (
auto k=0;k<this->numCubatures_;++k) {
78 const auto &cub = this->cubatures_[k];
79 tmpPoints [k] = cubPointViewType (
"CubatureTensor::getCubature::tmpPoints", cub.getNumPoints(), cub.getDimension());
80 tmpWeights[k] = cubWeightViewType(
"CubatureTensor::getCubature::tmpWeights", cub.getNumPoints());
81 cub.getCubature(tmpPoints[k], tmpWeights[k]);
88 for (
auto k=0;k<this->numCubatures_;++k) {
89 offset[k+1] = offset[k] + this->cubatures_[k].getDimension();
91 ordinal_type ii = 0, i[3] = {};
93 cubPointViewType cubPoints_device(
"cubPoints_device", cubPoints.extent(0), cubPoints.extent(1));
94 cubWeightViewType cubWeights_device(
"cubWeights_device", cubWeights.extent(0));
96 auto cubPoints_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cubPoints_device);
97 auto cubWeights_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cubWeights_device);
99 Kokkos::deep_copy(cubPoints_host, 0.0);
100 Kokkos::deep_copy(cubWeights_host, 1.0);
102 switch (this->numCubatures_) {
104 const ordinal_type npts[] = { this->cubatures_[0].getNumPoints(), this->cubatures_[1].getNumPoints() };
105 const ordinal_type dim [] = { this->cubatures_[0].getDimension(), this->cubatures_[1].getDimension() };
108 typename cubWeightViewType::HostMirror tmpWeights_host[2];
109 tmpWeights_host[0] = Kokkos::create_mirror_view(tmpWeights[0]);
110 tmpWeights_host[1] = Kokkos::create_mirror_view(tmpWeights[1]);
111 Kokkos::deep_copy(tmpWeights_host[0], tmpWeights[0]);
112 Kokkos::deep_copy(tmpWeights_host[1], tmpWeights[1]);
114 typename cubPointViewType::HostMirror tmpPoints_host[2];
115 tmpPoints_host[0] = Kokkos::create_mirror_view(tmpPoints[0]);
116 tmpPoints_host[1] = Kokkos::create_mirror_view(tmpPoints[1]);
118 Kokkos::deep_copy(tmpPoints_host[0], tmpPoints[0]);
119 Kokkos::deep_copy(tmpPoints_host[1], tmpPoints[1]);
121 for (i[1]=0;i[1]<npts[1];++i[1])
122 for (i[0]=0;i[0]<npts[0];++i[0]) {
123 for (
auto nc=0;nc<2;++nc) {
124 cubWeights_host(ii) *= tmpWeights_host[nc](i[nc]);
125 for (ordinal_type j=0;j<dim[nc];++j)
126 cubPoints_host(ii, offset[nc]+j) = tmpPoints_host[nc](i[nc], j);
133 const ordinal_type npts[] = { this->cubatures_[0].getNumPoints(), this->cubatures_[1].getNumPoints(), this->cubatures_[2].getNumPoints() };
134 const ordinal_type dim [] = { this->cubatures_[0].getDimension(), this->cubatures_[1].getDimension(), this->cubatures_[2].getDimension() };
137 typename cubWeightViewType::HostMirror tmpWeights_host[3];
138 tmpWeights_host[0] = Kokkos::create_mirror_view(tmpWeights[0]);
139 tmpWeights_host[1] = Kokkos::create_mirror_view(tmpWeights[1]);
140 tmpWeights_host[2] = Kokkos::create_mirror_view(tmpWeights[2]);
142 Kokkos::deep_copy(tmpWeights_host[0], tmpWeights[0]);
143 Kokkos::deep_copy(tmpWeights_host[1], tmpWeights[1]);
144 Kokkos::deep_copy(tmpWeights_host[2], tmpWeights[2]);
146 typename cubPointViewType::HostMirror tmpPoints_host[3];
147 tmpPoints_host[0] = Kokkos::create_mirror_view(tmpPoints[0]);
148 tmpPoints_host[1] = Kokkos::create_mirror_view(tmpPoints[1]);
149 tmpPoints_host[2] = Kokkos::create_mirror_view(tmpPoints[2]);
151 Kokkos::deep_copy(tmpPoints_host[0], tmpPoints[0]);
152 Kokkos::deep_copy(tmpPoints_host[1], tmpPoints[1]);
153 Kokkos::deep_copy(tmpPoints_host[2], tmpPoints[2]);
155 for (i[2]=0;i[2]<npts[2];++i[2])
156 for (i[1]=0;i[1]<npts[1];++i[1])
157 for (i[0]=0;i[0]<npts[0];++i[0]) {
158 for (
auto nc=0;nc<3;++nc) {
159 cubWeights_host(ii) *= tmpWeights_host[nc](i[nc]);
160 for (ordinal_type j=0;j<dim[nc];++j)
161 cubPoints_host(ii, offset[nc]+j) = tmpPoints_host[nc](i[nc], j);
168 INTREPID2_TEST_FOR_EXCEPTION( !(this->numCubatures_ == 2 || this->numCubatures_ == 3), std::runtime_error,
169 ">>> ERROR (CubatureTensor::getCubature): CubatureTensor supports only 2 or 3 component direct cubatures.");
172 Kokkos::deep_copy(cubPoints_device, cubPoints_host);
173 Kokkos::deep_copy(cubWeights_device, cubWeights_host);
175 Kokkos::deep_copy(cubPoints, cubPoints_device);
176 Kokkos::deep_copy(cubWeights, cubWeights_device);
void getCubatureImpl(Kokkos::DynRankView< cubPointValueType, cubPointProperties...> cubPoints, Kokkos::DynRankView< cubWeightValueType, cubWeightProperties...> cubWeights) const
static constexpr ordinal_type MaxTensorComponents
Maximum number of tensor/Cartesian products that can be taken: this allows hypercube basis in 7D to b...