Intrepid2
Intrepid2_CellToolsDefControlVolume.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
49#ifndef __INTREPID2_CELLTOOLS_DEF_CONTROL_VOLUME_HPP__
50#define __INTREPID2_CELLTOOLS_DEF_CONTROL_VOLUME_HPP__
51
52// disable clang warnings
53#if defined (__clang__) && !defined (__INTEL_COMPILER)
54#pragma clang system_header
55#endif
56
57namespace Intrepid2 {
58
59 //============================================================================================//
60 // //
61 // Control Volume Coordinates //
62 // //
63 //============================================================================================//
64
65 namespace FunctorCellTools {
66
71 template<typename centerViewType, typename vertViewType>
72 KOKKOS_INLINE_FUNCTION
73 void getBarycenterPolygon2D( centerViewType center,
74 const vertViewType verts) {
75 // the enumeration already assumes the ordering of vertices (circling around the polygon)
76 const ordinal_type nvert = verts.extent(0);
77
78 center(0) = 0;
79 center(1) = 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;
86 area += 0.5*scale;
87 }
88 center(0) /= (6.0*area);
89 center(1) /= (6.0*area);
90 }
91
92 template<typename midPointViewType, 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);
99
100 for (ordinal_type i=0;i<npts;++i) {
101 // first entry is the number of subcell vertices
102 const ordinal_type nvert_per_subcell = map(i, 0);
103 for (ordinal_type j=0;j<dim;++j) {
104 midpts(i,j) = 0;
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;
108 }
109 }
110 }
111
112 template<typename subcvCoordViewType,
113 typename cellCoordViewType,
114 typename mapViewType>
119 subcvCoordViewType _subcvCoords;
120 const cellCoordViewType _cellCoords;
121 const mapViewType _edgeMap;
122
123 KOKKOS_INLINE_FUNCTION
124 F_getSubcvCoords_Polygon2D( subcvCoordViewType subcvCoords_,
125 cellCoordViewType cellCoords_,
126 mapViewType edgeMap_ )
127 : _subcvCoords(subcvCoords_), _cellCoords(cellCoords_), _edgeMap(edgeMap_) {}
128
129 KOKKOS_INLINE_FUNCTION
130 void operator()(const ordinal_type cell) const {
131 // vertices of cell (P,D)
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);
136
137 // control volume coords (N,P,D), here N corresponds to cell vertices
138 auto cvCoords = Kokkos::subdynrankview( _subcvCoords, cell,
139 Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
140
141 // work space for barycenter and midpoints on edges
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);
146
147 getBarycenterPolygon2D(center, verts);
148 getMidPoints(midpts, _edgeMap, verts);
149
150 for (ordinal_type i=0;i<nvert;++i) {
151 for (ordinal_type j=0;j<dim;++j) {
152 // control volume is always quad
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);
157 }
158 }
159 }
160 };
161
162 template<typename subcvCoordViewType,
163 typename cellCoordViewType,
164 typename mapViewType>
169 subcvCoordViewType _subcvCoords;
170 const cellCoordViewType _cellCoords;
171 const mapViewType _edgeMap, _faceMap;
172
173 KOKKOS_INLINE_FUNCTION
174 F_getSubcvCoords_Hexahedron( subcvCoordViewType subcvCoords_,
175 cellCoordViewType cellCoords_,
176 mapViewType edgeMap_,
177 mapViewType faceMap_ )
178 : _subcvCoords(subcvCoords_), _cellCoords(cellCoords_), _edgeMap(edgeMap_), _faceMap(faceMap_) {}
179
180 KOKKOS_INLINE_FUNCTION
181 void operator()(const ordinal_type cell) const {
182 // vertices of cell (P,D)
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);
187
188 // control volume coords (N,P,D), here N corresponds to cell vertices
189 auto cvCoords = Kokkos::subdynrankview( _subcvCoords, cell,
190 Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
191
192 // work space for barycenter and midpoints on edges
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);
198
199 // find barycenter
200 //Warning! I think this assumes the Hexa is affinely mapped from the reference Hexa
201 for (ordinal_type j=0;j<3;++j) {
202 center(j) = 0;
203 for (ordinal_type i=0;i<nvert;++i)
204 center(j) += verts(i,j);
205 center(j) /= nvert;
206 }
207
208 getMidPoints(edge_midpts, _edgeMap, verts);
209 getMidPoints(face_midpts, _faceMap, verts);
210
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) {
214
215 // set first node of bottom hex to primary cell node
216 // and fifth node of upper hex
217 cvCoords(i, 0,j) = verts(i, j);
218 cvCoords(i+4,4,j) = verts(i+4,j);
219
220 // set second node of bottom hex to adjacent edge midpoint
221 // and sixth node of upper hex
222 cvCoords(i, 1,j) = edge_midpts(i, j);
223 cvCoords(i+4,5,j) = edge_midpts(i+4,j);
224
225 // set third node of bottom hex to bottom face midpoint (number 4)
226 // and seventh node of upper hex to top face midpoint
227 cvCoords(i, 2,j) = face_midpts(4,j);
228 cvCoords(i+4,6,j) = face_midpts(5,j);
229
230 // set fourth node of bottom hex to other adjacent edge midpoint
231 // and eight node of upper hex to other adjacent edge midpoint
232 cvCoords(i, 3,j) = edge_midpts(ii, j);
233 cvCoords(i+4,7,j) = edge_midpts(ii+4,j);
234
235 // set fifth node to vertical edge
236 // same as first node of upper hex
237 cvCoords(i, 4,j) = edge_midpts(i+8,j);
238 cvCoords(i+4,0,j) = edge_midpts(i+8,j);
239
240 // set sixth node to adjacent face midpoint
241 // same as second node of upper hex
242 cvCoords(i, 5,j) = face_midpts(i,j);
243 cvCoords(i+4,1,j) = face_midpts(i,j);
244
245 // set seventh node to barycenter
246 // same as third node of upper hex
247 cvCoords(i, 6,j) = center(j);
248 cvCoords(i+4,2,j) = center(j);
249
250 // set eighth node to other adjacent face midpoint
251 // same as fourth node of upper hex
252 cvCoords(i, 7,j) = face_midpts(ii,j);
253 cvCoords(i+4,3,j) = face_midpts(ii,j);
254 }
255 }
256 }
257 };
258
259
260 template<typename subcvCoordViewType,
261 typename cellCoordViewType,
262 typename mapViewType>
267 subcvCoordViewType _subcvCoords;
268 const cellCoordViewType _cellCoords;
269 const mapViewType _edgeMap, _faceMap;
270
271 KOKKOS_INLINE_FUNCTION
272 F_getSubcvCoords_Tetrahedron( subcvCoordViewType subcvCoords_,
273 cellCoordViewType cellCoords_,
274 mapViewType edgeMap_,
275 mapViewType faceMap_ )
276 : _subcvCoords(subcvCoords_), _cellCoords(cellCoords_), _edgeMap(edgeMap_), _faceMap(faceMap_) {}
277
278 KOKKOS_INLINE_FUNCTION
279 void operator()(const ordinal_type cell) const {
280 // ** vertices of cell (P,D)
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);
285
286 // control volume coords (N,P,D), here N corresponds to cell vertices
287 auto cvCoords = Kokkos::subdynrankview( _subcvCoords, cell,
288 Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL() );
289
290 // work space for barycenter and midpoints on edges
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);
296
297 // find barycenter
298 for (ordinal_type j=0;j<3;++j) {
299 center(j) = 0;
300 for (ordinal_type i=0;i<nvert;++i)
301 center(j) += verts(i,j);
302 center(j) /= nvert;
303 }
304
305 getMidPoints(edge_midpts, _edgeMap, verts);
306 getMidPoints(face_midpts, _faceMap, verts);
307
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) {
311 // set first node of bottom hex to primary cell node
312 cvCoords(i,0,j) = verts(i,j);
313
314 // set second node of bottom hex to adjacent edge midpoint
315 cvCoords(i,1,j) = edge_midpts(i,j);
316
317 // set third node of bottom hex to bottom face midpoint (number 3)
318 cvCoords(i,2,j) = face_midpts(3,j);
319
320 // set fourth node of bottom hex to other adjacent edge midpoint
321 cvCoords(i,3,j) = edge_midpts(ii,j);
322
323 // set fifth node to vertical edge
324 cvCoords(i,4,j) = edge_midpts(i+3,j);
325
326 // set sixth node to adjacent face midpoint
327 cvCoords(i,5,j) = face_midpts(i,j);
328
329 // set seventh node to barycenter
330 cvCoords(i,6,j) = center(j);
331
332 // set eighth node to other adjacent face midpoint
333 cvCoords(i,7,j) = face_midpts(ii,j);
334 }
335 }
336
337 for (ordinal_type j=0;j<dim;++j) {
338 // Control volume attached to fourth node
339 // set first node of bottom hex to primary cell node
340 cvCoords(3,0,j) = verts(3,j);
341
342 // set second node of bottom hex to adjacent edge midpoint
343 cvCoords(3,1,j) = edge_midpts(3,j);
344
345 // set third node of bottom hex to bottom face midpoint (number 3)
346 cvCoords(3,2,j) = face_midpts(2,j);
347
348 // set fourth node of bottom hex to other adjacent edge midpoint
349 cvCoords(3,3,j) = edge_midpts(5,j);
350
351 // set fifth node to vertical edge
352 cvCoords(3,4,j) = edge_midpts(4,j);
353
354 // set sixth node to adjacent face midpoint
355 cvCoords(3,5,j) = face_midpts(0,j);
356
357 // set seventh node to barycenter
358 cvCoords(3,6,j) = center(j);
359
360 // set eighth node to other adjacent face midpoint
361 cvCoords(3,7,j) = face_midpts(1,j);
362 }
363 }
364 };
365
366 }
367
368 template<typename DeviceType>
369 template<typename subcvCoordValueType, class ...subcvCoordProperties,
370 typename cellCoordValueType, class ...cellCoordProperties>
371 void
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." );
379
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." );
382
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." );
385#endif
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");
392
393
394 // get array dimensions
395 const ordinal_type numCells = cellCoords.extent(0);
396 //const ordinal_type numVerts = cellCoords.extent(1);
397 const ordinal_type spaceDim = cellCoords.extent(2);
398
399 // construct edge and face map for the cell type
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);
407 }
408
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);
416 }
417
418 Kokkos::deep_copy(edgeMap, edgeMapHost);
419 Kokkos::deep_copy(faceMap, faceMapHost);
420
421 // parallel run
422 using subcvCoordViewType = decltype(subcvCoords);
423 using cellCoordViewType = decltype(cellCoords);
424 using mapViewType = decltype(edgeMap);
425
426 const auto loopSize = numCells;
427 Kokkos::RangePolicy<ExecSpaceType, Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
428
429 switch (primaryCell.getKey()) {
430 case shards::Triangle<3>::key:
431 case shards::Quadrilateral<4>::key: {
432 // 2D polygon
434 Kokkos::parallel_for( policy, FunctorType(subcvCoords, cellCoords, edgeMap) );
435 break;
436 }
437 case shards::Hexahedron<8>::key: {
439 Kokkos::parallel_for( policy, FunctorType(subcvCoords, cellCoords, edgeMap, faceMap) );
440 break;
441 }
442 case shards::Tetrahedron<4>::key: {
444 Kokkos::parallel_for( policy, FunctorType(subcvCoords, cellCoords, edgeMap, faceMap) );
445 break;
446 }
447 default: {
448 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
449 ">>> ERROR (Intrepid2::CellTools::getSubcvCoords: the give cell topology is not supported.");
450 }
451 }
452 }
453}
454
455#endif
KOKKOS_INLINE_FUNCTION void getBarycenterPolygon2D(centerViewType center, const vertViewType verts)
Computes barycenter of polygonal cells.
static void getSubcvCoords(Kokkos::DynRankView< subcvCoordValueType, subcvCoordProperties... > subcvCoords, const Kokkos::DynRankView< cellCoordValueType, cellCoordProperties... > cellCoords, const shards::CellTopology primaryCell)
Computes coordinates of sub-control volumes in each primary cell.
Functor for calculation of sub-control volume coordinates on hexahedra see Intrepid2::CellTools for m...
Functor for calculation of sub-control volume coordinates on polygons see Intrepid2::CellTools for mo...
Functor for calculation of sub-control volume coordinates on tetrahedra see Intrepid2::CellTools for ...