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