Intrepid2
Intrepid2_CubatureTensorDef.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
49#ifndef __INTREPID2_CUBATURE_TENSOR_DEF_HPP__
50#define __INTREPID2_CUBATURE_TENSOR_DEF_HPP__
51
52namespace Intrepid2 {
53
54 template<typename DT, typename PT, typename WT>
55 template<typename cubPointValueType, class ...cubPointProperties,
56 typename cubWeightValueType, class ...cubWeightProperties>
57 void
59 getCubatureImpl( Kokkos::DynRankView<cubPointValueType, cubPointProperties...> cubPoints,
60 Kokkos::DynRankView<cubWeightValueType,cubWeightProperties...> cubWeights ) const {
61#ifdef HAVE_INTREPID2_DEBUG
62 // check size of cubPoints and cubWeights
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.");
67#endif
68 using cubPointViewType = Kokkos::DynRankView<cubPointValueType, DT>;
69 using cubWeightViewType = Kokkos::DynRankView<cubWeightValueType,DT>;
70
71 // mirroring and where the data is problematic... when it becomes a problem, then deal with it.
72 cubPointViewType tmpPoints [Parameters::MaxTensorComponents];
73 cubWeightViewType tmpWeights[Parameters::MaxTensorComponents];
74
75 // this temporary allocation can be member of cubature; for now, let's do this way.
76 // this is cubature setup on the reference cell and called for tensor elements.
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]);
82 }
84 // when the input containers are device space, this is better computed on host and copy to devices
85 // fill tensor cubature
86 {
87 ordinal_type offset[Parameters::MaxTensorComponents+1] = {};
88 for (auto k=0;k<this->numCubatures_;++k) {
89 offset[k+1] = offset[k] + this->cubatures_[k].getDimension();
90 }
91 ordinal_type ii = 0, i[3] = {};
92
93 cubPointViewType cubPoints_device("cubPoints_device", cubPoints.extent(0), cubPoints.extent(1));
94 cubWeightViewType cubWeights_device("cubWeights_device", cubWeights.extent(0));
95
96 auto cubPoints_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cubPoints_device);
97 auto cubWeights_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cubWeights_device);
98
99 Kokkos::deep_copy(cubPoints_host, 0.0);
100 Kokkos::deep_copy(cubWeights_host, 1.0);
101
102 switch (this->numCubatures_) {
103 case 2: {
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() };
106
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]);
113
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]);
117
118 Kokkos::deep_copy(tmpPoints_host[0], tmpPoints[0]);
119 Kokkos::deep_copy(tmpPoints_host[1], tmpPoints[1]);
120
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);
127 }
128 ++ii;
129 }
130 break;
131 }
132 case 3: {
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() };
135
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]);
141
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]);
145
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]);
150
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]);
154
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);
162 }
163 ++ii;
164 }
165 break;
166 }
167 default: {
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.");
170 }
171 }
172 Kokkos::deep_copy(cubPoints_device, cubPoints_host);
173 Kokkos::deep_copy(cubWeights_device, cubWeights_host);
174
175 Kokkos::deep_copy(cubPoints, cubPoints_device);
176 Kokkos::deep_copy(cubWeights, cubWeights_device);
177 }
178 }
179
180} // end namespace Intrepid2
181
182#endif
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...