Intrepid2
Intrepid2_OrientationToolsDefModifyBasis.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
48#ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_BASIS_HPP__
49#define __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_BASIS_HPP__
50
51// disable clang warnings
52#if defined (__clang__) && !defined (__INTEL_COMPILER)
53#pragma clang system_header
54#endif
55
57
58namespace Intrepid2 {
59
60 template<typename DT>
61 template<typename elemOrtValueType, class ...elemOrtProperties,
62 typename elemNodeValueType, class ...elemNodeProperties>
63 void
65 getOrientation( Kokkos::DynRankView<elemOrtValueType,elemOrtProperties...> elemOrts,
66 const Kokkos::DynRankView<elemNodeValueType,elemNodeProperties...> elemNodes,
67 const shards::CellTopology cellTopo,
68 bool isSide) {
69 // small meta data modification and it uses shards; let's do this on host
70 auto elemOrtsHost = Kokkos::create_mirror_view(elemOrts);
71 auto elemNodesHost = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elemNodes);
72
73 const ordinal_type numCells = elemNodes.extent(0);
74 for (auto cell=0;cell<numCells;++cell) {
75 const auto nodes = Kokkos::subview(elemNodesHost, cell, Kokkos::ALL());
76 elemOrtsHost(cell) = Orientation::getOrientation(cellTopo, nodes, isSide);
77 }
78
79 Kokkos::deep_copy(elemOrts, elemOrtsHost);
80 }
81
82 template<typename ortViewType,
83 typename OutputViewType,
84 typename inputViewType,
85 typename o2tViewType,
86 typename t2oViewType,
87 typename dataViewType>
89 ortViewType orts;
90 OutputViewType output;
91 inputViewType input;
92 o2tViewType ordinalToTag;
93 t2oViewType tagToOrdinal;
94
95 const dataViewType matData;
96 const ordinal_type cellDim, numVerts, numEdges, numFaces, numPoints, dimBasis;
97
98 F_modifyBasisByOrientation(ortViewType orts_,
99 OutputViewType output_,
100 inputViewType input_,
101 o2tViewType ordinalToTag_,
102 t2oViewType tagToOrdinal_,
103 const dataViewType matData_,
104 const ordinal_type cellDim_,
105 const ordinal_type numVerts_,
106 const ordinal_type numEdges_,
107 const ordinal_type numFaces_,
108 const ordinal_type numPoints_,
109 const ordinal_type dimBasis_)
110 : orts(orts_),
111 output(output_),
112 input(input_),
113 ordinalToTag(ordinalToTag_),
114 tagToOrdinal(tagToOrdinal_),
115 matData(matData_),
116 cellDim(cellDim_),
117 numVerts(numVerts_),
118 numEdges(numEdges_),
119 numFaces(numFaces_),
120 numPoints(numPoints_),
121 dimBasis(dimBasis_)
122 {}
123
124 KOKKOS_INLINE_FUNCTION
125 void operator()(const ordinal_type cell) const {
126 typedef typename inputViewType::non_const_value_type input_value_type;
127
128 auto out = Kokkos::subview(output, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
129 auto in = (input.rank() == output.rank()) ?
130 Kokkos::subview(input, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL())
131 : Kokkos::subview(input, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
132
133 // edge transformation
134 ordinal_type existEdgeDofs = 0;
135 if (numEdges > 0) {
136 ordinal_type ortEdges[12];
137 orts(cell).getEdgeOrientation(ortEdges, numEdges);
138
139 // apply coeff matrix
140 for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
141 const ordinal_type ordEdge = (1 < tagToOrdinal.extent(0) ? (static_cast<size_type>(edgeId) < tagToOrdinal.extent(1) ? tagToOrdinal(1, edgeId, 0) : -1) : -1);
142
143 if (ordEdge != -1) {
144 existEdgeDofs = 1;
145 const ordinal_type ndofEdge = ordinalToTag(ordEdge, 3);
146 const auto mat = Kokkos::subview(matData,
147 edgeId, ortEdges[edgeId],
148 Kokkos::ALL(), Kokkos::ALL());
149
150 for (ordinal_type j=0;j<numPoints;++j)
151 for (ordinal_type i=0;i<ndofEdge;++i) {
152 const ordinal_type ii = tagToOrdinal(1, edgeId, i);
153
154 for (ordinal_type k=0;k<dimBasis;++k) {
155 input_value_type temp = 0.0;
156 for (ordinal_type l=0;l<ndofEdge;++l) {
157 const ordinal_type ll = tagToOrdinal(1, edgeId, l);
158 temp += mat(i,l)*in(ll, j, k);
159 }
160 out(ii, j, k) = temp;
161 }
162 }
163 }
164 }
165 }
166
167 // face transformation
168 if (numFaces > 0) {
169 ordinal_type ortFaces[12];
170 orts(cell).getFaceOrientation(ortFaces, numFaces);
171
172 // apply coeff matrix
173 for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
174 const ordinal_type ordFace = (2 < tagToOrdinal.extent(0) ? (static_cast<size_type>(faceId) < tagToOrdinal.extent(1) ? tagToOrdinal(2, faceId, 0) : -1) : -1);
175
176 if (ordFace != -1) {
177 const ordinal_type ndofFace = ordinalToTag(ordFace, 3);
178 const auto mat = Kokkos::subview(matData,
179 numEdges*existEdgeDofs+faceId, ortFaces[faceId],
180 Kokkos::ALL(), Kokkos::ALL());
181
182 for (ordinal_type j=0;j<numPoints;++j)
183 for (ordinal_type i=0;i<ndofFace;++i) {
184 const ordinal_type ii = tagToOrdinal(2, faceId, i);
185
186 for (ordinal_type k=0;k<dimBasis;++k) {
187 input_value_type temp = 0.0;
188 for (ordinal_type l=0;l<ndofFace;++l) {
189 const ordinal_type ll = tagToOrdinal(2, faceId, l);
190 temp += mat(i,l)*in(ll, j, k);
191 }
192 out(ii, j, k) = temp;
193 }
194 }
195 }
196 }
197 }
198
199 //side orientations
200 ordinal_type faceOrt(0), edgeOrt(0);
201 if(cellDim == 2) orts(cell).getFaceOrientation(&faceOrt, 1);
202 if (faceOrt != 0) {
203 const ordinal_type ordFace = (2 < tagToOrdinal.extent(0) ? (static_cast<size_type>(0) < tagToOrdinal.extent(1) ? tagToOrdinal(2, 0, 0) : -1) : -1);
204
205 if (ordFace != -1) {
206 const ordinal_type ndofFace = ordinalToTag(ordFace, 3);
207 const auto mat = Kokkos::subview(matData,
208 numEdges*existEdgeDofs, faceOrt,
209 Kokkos::ALL(), Kokkos::ALL());
210
211 for (ordinal_type j=0;j<numPoints;++j)
212 for (ordinal_type i=0;i<ndofFace;++i) {
213 const ordinal_type ii = tagToOrdinal(2, 0, i);
214
215 for (ordinal_type k=0;k<dimBasis;++k) {
216 input_value_type temp = 0.0;
217 for (ordinal_type l=0;l<ndofFace;++l) {
218 const ordinal_type ll = tagToOrdinal(2, 0, l);
219 temp += mat(i,l)*in(ll, j, k);
220 }
221 out(ii, j, k) = temp;
222 }
223 }
224 }
225 }
226
227 if(cellDim == 1) orts(cell).getEdgeOrientation(&edgeOrt, 1);
228 if (edgeOrt != 0) {
229 const ordinal_type ordEdge = (1 < tagToOrdinal.extent(0) ? (static_cast<size_type>(0) < tagToOrdinal.extent(1) ? tagToOrdinal(1, 0, 0) : -1) : -1);
230
231 if (ordEdge != -1) {
232 const ordinal_type ndofEdge = ordinalToTag(ordEdge, 3);
233 const auto mat = Kokkos::subview(matData,
234 0, edgeOrt,
235 Kokkos::ALL(), Kokkos::ALL());
236
237 for (ordinal_type j=0;j<numPoints;++j)
238 for (ordinal_type i=0;i<ndofEdge;++i) {
239 const ordinal_type ii = tagToOrdinal(1, 0, i);
240
241 for (ordinal_type k=0;k<dimBasis;++k) {
242 input_value_type temp = 0.0;
243 for (ordinal_type l=0;l<ndofEdge;++l) {
244 const ordinal_type ll = tagToOrdinal(1, 0, l);
245 temp += mat(i,l)*in(ll, j, k);
246 }
247 out(ii, j, k) = temp;
248 }
249 }
250 }
251 }
252 }
253 };
254
255 template<typename DT>
256 template<typename outputValueType, class ...outputProperties,
257 typename inputValueType, class ...inputProperties,
258 typename OrientationViewType,
259 typename BasisType>
260 void
262 modifyBasisByOrientation( Kokkos::DynRankView<outputValueType,outputProperties...> output,
263 const Kokkos::DynRankView<inputValueType, inputProperties...> input,
264 const OrientationViewType orts,
265 const BasisType* basis ) {
266#ifdef HAVE_INTREPID2_DEBUG
267 {
268 if (input.rank() == output.rank())
269 {
270 for (size_type i=0;i<input.rank();++i)
271 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i), std::invalid_argument,
272 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input and output dimension does not match.");
273 }
274 else if (input.rank() == output.rank() - 1)
275 {
276 for (size_type i=0;i<input.rank();++i)
277 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i+1), std::invalid_argument,
278 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input dimensions must match output dimensions exactly, or else match all but the first dimension (in the case that input does not have a 'cell' dimension).");
279 }
280 else
281 {
282 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument,
283 ">>> ERROR (OrientationTools::modifyBasisByOrientation): input and output ranks must either match, or input rank must be one less than that of output.")
284 }
285
286 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(1)) != basis->getCardinality(), std::invalid_argument,
287 ">>> ERROR (OrientationTools::modifyBasisByOrientation): Field dimension of input/output does not match to basis cardinality.");
288 }
289#endif
290
291 const shards::CellTopology cellTopo = basis->getBaseCellTopology();
292 const ordinal_type cellDim = cellTopo.getDimension();
293
294 //Initialize output with values from input
295 if(input.rank() == output.rank())
296 Kokkos::deep_copy(output, input);
297 else
298 RealSpaceTools<DT>::clone(output, input);
299
300 if ((cellDim < 3) || basis->requireOrientation()) {
301 auto ordinalToTag = Kokkos::create_mirror_view_and_copy(typename DT::memory_space(), basis->getAllDofTags());
302 auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(typename DT::memory_space(), basis->getAllDofOrdinal());
303
304 const ordinal_type
305 numCells = output.extent(0),
306 //numBasis = output.extent(1),
307 numPoints = output.extent(2),
308 dimBasis = output.extent(3); //returns 1 when output.rank() < 4;
309
310 const CoeffMatrixDataViewType matData = createCoeffMatrix(basis);
311
312 ordinal_type numVerts(0), numEdges(0), numFaces(0);
313
314 if (basis->requireOrientation()) {
315 numVerts = cellTopo.getVertexCount()*ordinal_type(basis->getDofCount(0, 0) > 0);
316 numEdges = cellTopo.getEdgeCount()*ordinal_type(basis->getDofCount(1, 0) > 0);
317 numFaces = cellTopo.getFaceCount()*ordinal_type(basis->getDofCount(2, 0) > 0);
318 }
319
320 const Kokkos::RangePolicy<typename DT::execution_space> policy(0, numCells);
322 <decltype(orts),
323 decltype(output),decltype(input),
324 decltype(ordinalToTag),decltype(tagToOrdinal),
325 decltype(matData)> FunctorType;
326 Kokkos::parallel_for
327 (policy,
328 FunctorType(orts,
329 output, input,
330 ordinalToTag, tagToOrdinal,
331 matData,
332 cellDim, numVerts, numEdges, numFaces,
333 numPoints, dimBasis));
334 }
335 }
336}
337
338#endif
Header file for the Intrepid2::Orientation class.
Kokkos::View< double ****, DeviceType > CoeffMatrixDataViewType
subcell ordinal, orientation, matrix m x n
static void modifyBasisByOrientation(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input, const OrientationViewType orts, const BasisType *basis)
Modify basis due to orientation.
static void getOrientation(Kokkos::DynRankView< elemOrtValueType, elemOrtProperties... > elemOrts, const Kokkos::DynRankView< elemNodeValueType, elemNodeProperties... > elemNodes, const shards::CellTopology cellTopo, bool isSide=false)
Compute orientations of cells in a workset.
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.