49 #ifndef __INTREPID2_CELLTOOLS_DEF_CONTROL_VOLUME_HPP__
50 #define __INTREPID2_CELLTOOLS_DEF_CONTROL_VOLUME_HPP__
53 #if defined (__clang__) && !defined (__INTEL_COMPILER)
54 #pragma clang system_header
65 namespace FunctorCellTools {
71 template<
typename centerViewType,
typename vertViewType>
72 KOKKOS_INLINE_FUNCTION
73 void getBarycenterPolygon2D( centerViewType center,
74 const vertViewType verts) {
76 const ordinal_type nvert = verts.extent(0);
80 typename centerViewType::value_type area = 0;
81 for (ordinal_type i=0;i<nvert;++i) {
82 const ordinal_type j = (i + 1)%nvert;
83 const auto scale = verts(i,0)*verts(j,1) - verts(j,0)*verts(i,1);
84 center(0) += (verts(i,0) + verts(j,0))*scale;
85 center(1) += (verts(i,1) + verts(j,1))*scale;
88 center(0) /= (6.0*area);
89 center(1) /= (6.0*area);
92 template<
typename m
idPo
intViewType,
typename nodeMapViewType,
typename vertViewType>
93 KOKKOS_INLINE_FUNCTION
94 void getMidPoints( midPointViewType midpts,
95 const nodeMapViewType map,
96 const vertViewType verts) {
97 const ordinal_type npts = map.extent(0);
98 const ordinal_type dim = verts.extent(1);
100 for (ordinal_type i=0;i<npts;++i) {
102 const ordinal_type nvert_per_subcell = map(i, 0);
103 for (ordinal_type j=0;j<dim;++j) {
105 for (ordinal_type k=1;k<=nvert_per_subcell;++k)
106 midpts(i,j) += verts(map(i,k),j);
107 midpts(i,j) /= nvert_per_subcell;
112 template<
typename subcvCoordViewType,
113 typename cellCoordViewType,
114 typename mapViewType>
119 subcvCoordViewType _subcvCoords;
120 const cellCoordViewType _cellCoords;
121 const mapViewType _edgeMap;
123 KOKKOS_INLINE_FUNCTION
125 cellCoordViewType cellCoords_,
126 mapViewType edgeMap_ )
127 : _subcvCoords(subcvCoords_), _cellCoords(cellCoords_), _edgeMap(edgeMap_) {}
129 KOKKOS_INLINE_FUNCTION
130 void operator()(
const ordinal_type cell)
const {
132 const auto verts = Kokkos::subdynrankview( _cellCoords, cell,
133 Kokkos::ALL(), Kokkos::ALL() );
134 const ordinal_type nvert = verts.extent(0);
135 const ordinal_type dim = verts.extent(1);
138 auto cvCoords = Kokkos::subdynrankview( _subcvCoords, cell,
139 Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
142 typedef typename subcvCoordViewType::value_type value_type;
143 value_type buf_center[2], buf_midpts[4*2];
144 Kokkos::View<value_type*,Kokkos::AnonymousSpace> center(buf_center, 2);
145 Kokkos::View<value_type**,Kokkos::AnonymousSpace> midpts(buf_midpts, 4, 2);
147 getBarycenterPolygon2D(center, verts);
148 getMidPoints(midpts, _edgeMap, verts);
150 for (ordinal_type i=0;i<nvert;++i) {
151 for (ordinal_type j=0;j<dim;++j) {
153 cvCoords(i, 0, j) = verts(i, j);
154 cvCoords(i, 1, j) = midpts(i, j);
155 cvCoords(i, 2, j) = center(j);
156 cvCoords(i, 3, j) = midpts((i+nvert-1)%nvert, j);
162 template<
typename subcvCoordViewType,
163 typename cellCoordViewType,
164 typename mapViewType>
169 subcvCoordViewType _subcvCoords;
170 const cellCoordViewType _cellCoords;
171 const mapViewType _edgeMap, _faceMap;
173 KOKKOS_INLINE_FUNCTION
175 cellCoordViewType cellCoords_,
176 mapViewType edgeMap_,
177 mapViewType faceMap_ )
178 : _subcvCoords(subcvCoords_), _cellCoords(cellCoords_), _edgeMap(edgeMap_), _faceMap(faceMap_) {}
180 KOKKOS_INLINE_FUNCTION
181 void operator()(
const ordinal_type cell)
const {
183 const auto verts = Kokkos::subdynrankview( _cellCoords, cell,
184 Kokkos::ALL(), Kokkos::ALL() );
185 const ordinal_type nvert = verts.extent(0);
186 const ordinal_type dim = verts.extent(1);
189 auto cvCoords = Kokkos::subdynrankview( _subcvCoords, cell,
190 Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
193 typedef typename subcvCoordViewType::value_type value_type;
194 value_type buf_center[3], buf_edge_midpts[12*3], buf_face_midpts[6*3];
195 Kokkos::View<value_type*,Kokkos::AnonymousSpace> center(buf_center, 3);
196 Kokkos::View<value_type**,Kokkos::AnonymousSpace> edge_midpts(buf_edge_midpts, 12, 3);
197 Kokkos::View<value_type**,Kokkos::AnonymousSpace> face_midpts(buf_face_midpts, 6, 3);
201 for (ordinal_type j=0;j<3;++j) {
203 for (ordinal_type i=0;i<nvert;++i)
204 center(j) += verts(i,j);
208 getMidPoints(edge_midpts, _edgeMap, verts);
209 getMidPoints(face_midpts, _faceMap, verts);
211 for (ordinal_type i=0;i<4;++i) {
212 const ordinal_type ii = (i+4-1)%4;
213 for (ordinal_type j=0;j<dim;++j) {
217 cvCoords(i, 0,j) = verts(i, j);
218 cvCoords(i+4,4,j) = verts(i+4,j);
222 cvCoords(i, 1,j) = edge_midpts(i, j);
223 cvCoords(i+4,5,j) = edge_midpts(i+4,j);
227 cvCoords(i, 2,j) = face_midpts(4,j);
228 cvCoords(i+4,6,j) = face_midpts(5,j);
232 cvCoords(i, 3,j) = edge_midpts(ii, j);
233 cvCoords(i+4,7,j) = edge_midpts(ii+4,j);
237 cvCoords(i, 4,j) = edge_midpts(i+8,j);
238 cvCoords(i+4,0,j) = edge_midpts(i+8,j);
242 cvCoords(i, 5,j) = face_midpts(i,j);
243 cvCoords(i+4,1,j) = face_midpts(i,j);
247 cvCoords(i, 6,j) = center(j);
248 cvCoords(i+4,2,j) = center(j);
252 cvCoords(i, 7,j) = face_midpts(ii,j);
253 cvCoords(i+4,3,j) = face_midpts(ii,j);
260 template<
typename subcvCoordViewType,
261 typename cellCoordViewType,
262 typename mapViewType>
267 subcvCoordViewType _subcvCoords;
268 const cellCoordViewType _cellCoords;
269 const mapViewType _edgeMap, _faceMap;
271 KOKKOS_INLINE_FUNCTION
273 cellCoordViewType cellCoords_,
274 mapViewType edgeMap_,
275 mapViewType faceMap_ )
276 : _subcvCoords(subcvCoords_), _cellCoords(cellCoords_), _edgeMap(edgeMap_), _faceMap(faceMap_) {}
278 KOKKOS_INLINE_FUNCTION
279 void operator()(
const ordinal_type cell)
const {
281 const auto verts = Kokkos::subdynrankview( _cellCoords, cell,
282 Kokkos::ALL(), Kokkos::ALL() );
283 const ordinal_type nvert = verts.extent(0);
284 const ordinal_type dim = verts.extent(1);
287 auto cvCoords = Kokkos::subdynrankview( _subcvCoords, cell,
288 Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
291 typedef typename subcvCoordViewType::value_type value_type;
292 value_type buf_center[3], buf_edge_midpts[6*3], buf_face_midpts[4*3];
293 Kokkos::View<value_type*,Kokkos::AnonymousSpace> center(buf_center, 3);
294 Kokkos::View<value_type**,Kokkos::AnonymousSpace> edge_midpts(buf_edge_midpts, 6, 3);
295 Kokkos::View<value_type**,Kokkos::AnonymousSpace> face_midpts(buf_face_midpts, 4, 3);
298 for (ordinal_type j=0;j<3;++j) {
300 for (ordinal_type i=0;i<nvert;++i)
301 center(j) += verts(i,j);
305 getMidPoints(edge_midpts, _edgeMap, verts);
306 getMidPoints(face_midpts, _faceMap, verts);
308 for (ordinal_type i=0;i<3;++i) {
309 const ordinal_type ii = (i+3-1)%3;
310 for (ordinal_type j=0;j<dim;++j) {
312 cvCoords(i,0,j) = verts(i,j);
315 cvCoords(i,1,j) = edge_midpts(i,j);
318 cvCoords(i,2,j) = face_midpts(3,j);
321 cvCoords(i,3,j) = edge_midpts(ii,j);
324 cvCoords(i,4,j) = edge_midpts(i+3,j);
327 cvCoords(i,5,j) = face_midpts(i,j);
330 cvCoords(i,6,j) = center(j);
333 cvCoords(i,7,j) = face_midpts(ii,j);
337 for (ordinal_type j=0;j<dim;++j) {
340 cvCoords(3,0,j) = verts(3,j);
343 cvCoords(3,1,j) = edge_midpts(3,j);
346 cvCoords(3,2,j) = face_midpts(2,j);
349 cvCoords(3,3,j) = edge_midpts(5,j);
352 cvCoords(3,4,j) = edge_midpts(4,j);
355 cvCoords(3,5,j) = face_midpts(0,j);
358 cvCoords(3,6,j) = center(j);
361 cvCoords(3,7,j) = face_midpts(1,j);
368 template<
typename DeviceType>
369 template<
typename subcvCoordValueType,
class ...subcvCoordProperties,
370 typename cellCoordValueType,
class ...cellCoordProperties>
373 getSubcvCoords( Kokkos::DynRankView<subcvCoordValueType,subcvCoordProperties...> subcvCoords,
374 const Kokkos::DynRankView<cellCoordValueType,cellCoordProperties...> cellCoords,
375 const shards::CellTopology primaryCell ) {
376 #ifdef HAVE_INTREPID2_DEBUG
377 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(primaryCell), std::invalid_argument,
378 ">>> ERROR (Intrepid2::CellTools::getSubcvCoords): the primary cell must have a reference cell." );
380 INTREPID2_TEST_FOR_EXCEPTION( cellCoords.extent(1) != primaryCell.getVertexCount(), std::invalid_argument,
381 ">>> ERROR (Intrepid2::CellTools::getSubcvCoords): cell coords dimension(1) does not match to # of vertices of the cell." );
383 INTREPID2_TEST_FOR_EXCEPTION( cellCoords.extent(2) != primaryCell.getDimension(), std::invalid_argument,
384 ">>> ERROR (Intrepid2::CellTools::getSubcvCoords): cell coords dimension(2) does not match to the dimension of the cell." );
386 constexpr
bool are_accessible =
387 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
388 typename decltype(subcvCoords)::memory_space>::accessible &&
389 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
390 typename decltype(cellCoords)::memory_space>::accessible;
391 static_assert(are_accessible,
"CellTools<DeviceType>::getSubcvCoords(..): input/output views' memory spaces are not compatible with DeviceType");
395 const ordinal_type numCells = cellCoords.extent(0);
397 const ordinal_type spaceDim = cellCoords.extent(2);
400 const ordinal_type numEdge = primaryCell.getSubcellCount(1);
401 Kokkos::View<ordinal_type**,DeviceType> edgeMap(
"CellTools::getSubcvCoords::edgeMap", numEdge, 3);
402 auto edgeMapHost = Kokkos::create_mirror_view(edgeMap);
403 for (ordinal_type i=0;i<numEdge;++i) {
404 edgeMapHost(i,0) = primaryCell.getNodeCount(1, i);
405 for (ordinal_type j=0;j<edgeMapHost(i,0);++j)
406 edgeMapHost(i,j+1) = primaryCell.getNodeMap(1, i, j);
409 const ordinal_type numFace = (spaceDim > 2 ? primaryCell.getSubcellCount(2) : 0);
410 Kokkos::View<ordinal_type**,DeviceType> faceMap(
"CellTools::getSubcvCoords::faceMap", numFace, 5);
411 auto faceMapHost = Kokkos::create_mirror_view(faceMap);
412 for (ordinal_type i=0;i<numFace;++i) {
413 faceMapHost(i,0) = primaryCell.getNodeCount(2, i);
414 for (ordinal_type j=0;j<faceMapHost(i,0);++j)
415 faceMapHost(i,j+1) = primaryCell.getNodeMap(2, i, j);
418 Kokkos::deep_copy(edgeMap, edgeMapHost);
419 Kokkos::deep_copy(faceMap, faceMapHost);
422 using subcvCoordViewType = decltype(subcvCoords);
423 using cellCoordViewType = decltype(cellCoords);
424 using mapViewType = decltype(edgeMap);
426 const auto loopSize = numCells;
427 Kokkos::RangePolicy<ExecSpaceType, Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
429 switch (primaryCell.getKey()) {
430 case shards::Triangle<3>::key:
431 case shards::Quadrilateral<4>::key: {
434 Kokkos::parallel_for( policy, FunctorType(subcvCoords, cellCoords, edgeMap) );
437 case shards::Hexahedron<8>::key: {
439 Kokkos::parallel_for( policy, FunctorType(subcvCoords, cellCoords, edgeMap, faceMap) );
442 case shards::Tetrahedron<4>::key: {
444 Kokkos::parallel_for( policy, FunctorType(subcvCoords, cellCoords, edgeMap, faceMap) );
448 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
449 ">>> ERROR (Intrepid2::CellTools::getSubcvCoords: the give cell topology is not supported.");