Intrepid2
Intrepid2_CellTools_Serial.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
48#ifndef __INTREPID2_CELLTOOLS_SERIAL_HPP__
49#define __INTREPID2_CELLTOOLS_SERIAL_HPP__
50
51#include "Intrepid2_ConfigDefs.hpp"
52
53#include "Shards_CellTopology.hpp"
54
55#include "Intrepid2_Types.hpp"
56#include "Intrepid2_Utils.hpp"
57#include "Intrepid2_Kernels.hpp"
58
60
61namespace Intrepid2 {
62
63namespace Impl {
64
69class CellTools {
70public:
71
72 struct Serial {
73
74 // output:
75 // jacobian (physDim,refDim) - jacobian matrix evaluated at a single point
76 // input:
77 // grads (numNodes,refDim) - hgrad basis grad values evaluated at a single point (C1/C2 element only)
78 // nodes (numNodes,physDim) - cell element-to-node connectivity
79 template<typename jacobianViewType,
80 typename basisGradViewType,
81 typename nodeViewType>
82 KOKKOS_INLINE_FUNCTION
83 static void
84 computeJacobian(const jacobianViewType &jacobian, // physDim,refDim
85 const basisGradViewType &grads, // numNodes,refDim
86 const nodeViewType &nodes) { // numNodes,physDim
87 const auto numNodes = nodes.extent(0);
88
89 const auto physDim = jacobian.extent(0);
90 const auto refDim = jacobian.extent(1);
91
92 INTREPID2_TEST_FOR_ABORT( numNodes != grads.extent(0), "grad dimension_0 does not match to cardinality.");
93 INTREPID2_TEST_FOR_ABORT(refDim != grads.extent(1), "grad dimension_1 does not match to space dim.");
94 INTREPID2_TEST_FOR_ABORT( physDim != nodes.extent(1), "node dimension_1 does not match to space dim.");
95
96 Kernels::Serial::gemm_trans_notrans(1.0, nodes, grads, 0.0, jacobian);
97 }
98
99 // output:
100 // point (physDim) - mapped physical point
101 // input:
102 // vals (numNodes) - hgrad basis values evaluated at a single point (C1/C2 element only)
103 // nodes (numNodes,physDim) - cell element-to-node connectivity
104 template<typename PointViewType,
105 typename basisValViewType,
106 typename nodeViewType>
107 KOKKOS_INLINE_FUNCTION
108 static void
109 mapToPhysicalFrame(const PointViewType &point, // physDim
110 const basisValViewType &vals, // numNodes
111 const nodeViewType &nodes) { // numNodes,physDim
112 const auto numNodes = vals.extent(0);
113 const auto physDim = point.extent(0);
114
115 INTREPID2_TEST_FOR_ABORT(numNodes != nodes.extent(0), "nodes dimension_0 does not match to vals dimension_0.");
116 INTREPID2_TEST_FOR_ABORT(physDim != nodes.extent(1), "node dimension_1 does not match to space dim.");
117
118 Kernels::Serial::gemv_trans(1.0, nodes, vals, 0.0, point);
119 }
120
146 template<typename implBasisType,
147 typename refPointViewType,
148 typename phyPointViewType,
149 typename nodeViewType>
150 KOKKOS_INLINE_FUNCTION
151 static bool
152 mapToReferenceFrame(const refPointViewType &refPoint, // refDim
153 const phyPointViewType &physPoint, // physDim
154 const nodeViewType &nodes,
155 const double tol = tolerence()) { // numNodes,physDim
156 const ordinal_type refDim = refPoint.extent(0);
157 const ordinal_type physDim = physPoint.extent(0);
158 const ordinal_type numNodes = nodes.extent(0);
159
160 INTREPID2_TEST_FOR_ABORT(refDim > physDim, "the dimension of the reference cell is greater than physical cell dimension.");
161 INTREPID2_TEST_FOR_ABORT(physDim != static_cast<ordinal_type>(nodes.extent(1)), "physPoint dimension_0 does not match to space dim.");
162 INTREPID2_TEST_FOR_ABORT(numNodes > 27, "function hard-coded to support at most mappings with 27 Dofs");
163
164 typedef typename refPointViewType::non_const_value_type value_type;
165
166 // I want to use view instead of dynrankview
167 // NMAX = 27, MAXDIM = 3
168 value_type buf[27*3 + 27 + 9 + 9 + 9 + 9 + 3 + 3] = {}, *ptr = &buf[0];
169 Kokkos::DynRankView<value_type,
170 Kokkos::AnonymousSpace,
171 Kokkos::MemoryUnmanaged> grads(ptr, numNodes, refDim); ptr += numNodes*refDim;
172
173 Kokkos::DynRankView<value_type,
174 Kokkos::AnonymousSpace,
175 Kokkos::MemoryUnmanaged> vals(ptr, numNodes); ptr += numNodes;
176
177 Kokkos::DynRankView<value_type,
178 Kokkos::AnonymousSpace,
179 Kokkos::MemoryUnmanaged> jac(ptr, physDim, refDim); ptr += physDim*refDim;
180
181 Kokkos::DynRankView<value_type,
182 Kokkos::AnonymousSpace,
183 Kokkos::MemoryUnmanaged> metric(ptr, refDim, refDim); ptr += refDim*refDim;
184
185 Kokkos::DynRankView<value_type,
186 Kokkos::AnonymousSpace,
187 Kokkos::MemoryUnmanaged> invMetric(ptr, refDim, refDim); ptr += refDim*refDim;
188
189 Kokkos::DynRankView<value_type,
190 Kokkos::AnonymousSpace,
191 Kokkos::MemoryUnmanaged> invDf(ptr, refDim, physDim); ptr += refDim*physDim;
192
193 Kokkos::DynRankView<value_type,
194 Kokkos::AnonymousSpace,
195 Kokkos::MemoryUnmanaged> tmpPhysPoint(ptr, physDim); ptr += physDim;
196
197 Kokkos::DynRankView<value_type,
198 Kokkos::AnonymousSpace,
199 Kokkos::MemoryUnmanaged> oldRefPoint(ptr, refDim); ptr += refDim;
200
201 // set initial guess
202 for (ordinal_type j=0;j<refDim;++j) oldRefPoint(j) = 0;
203
204 double xphyNorm = Kernels::Serial::norm(physPoint, NORM_ONE);
205
206 bool converged = false;
207 for (ordinal_type iter=0;iter<Parameters::MaxNewton;++iter) {
208 // xtmp := F(oldRefPoint);
209 implBasisType::template Serial<OPERATOR_VALUE>::getValues(vals, oldRefPoint);
210 mapToPhysicalFrame(tmpPhysPoint, vals, nodes);
211 Kernels::Serial::z_is_axby(tmpPhysPoint, 1.0, physPoint, -1.0, tmpPhysPoint); //residual xtmp := physPoint - F(oldRefPoint);
212 double residualNorm = Kernels::Serial::norm(tmpPhysPoint, NORM_ONE);
213 if (residualNorm < tol*xphyNorm) {
214 converged = true;
215 break;
216 }
217
218 // DF^{-1}
219 implBasisType::template Serial<OPERATOR_GRAD>::getValues(grads, oldRefPoint);
220 CellTools::Serial::computeJacobian(jac, grads, nodes);
221
222 if(physDim == refDim) { //jacobian inverse
223 Kernels::Serial::inverse(invDf, jac);
224 } else { //jacobian is not square, we compute the pseudo-inverse (jac^T jac)^{-1} jac^T
225 Kernels::Serial::gemm_trans_notrans(1.0, jac, jac, 0.0, metric);
226 Kernels::Serial::inverse(invMetric, metric);
227 Kernels::Serial::gemm_notrans_trans(1.0, invMetric, jac, 0.0, invDf);
228 }
229
230 // Newton
231 Kernels::Serial::gemv_notrans(1.0, invDf, tmpPhysPoint, 0.0, refPoint); // refPoint := DF^{-1}( physPoint - F(oldRefPoint))
232 Kernels::Serial::z_is_axby(refPoint, 1.0, oldRefPoint, 1.0, refPoint); // refPoint += oldRefPoint
233
234 // temporarily overwriting oldRefPoint with oldRefPoint-refPoint
235 Kernels::Serial::z_is_axby(oldRefPoint, 1.0, oldRefPoint, -1.0, refPoint);
236
237 double err = Kernels::Serial::norm(oldRefPoint, NORM_ONE);
238 if (err < tol) {
239 converged = true;
240 break;
241 }
242
243 Kernels::Serial::copy(oldRefPoint, refPoint);
244 }
245 return converged;
246 }
247
248 template<typename refEdgeTanViewType, typename ParamViewType>
249 KOKKOS_INLINE_FUNCTION
250 static void
251 getReferenceEdgeTangent(const refEdgeTanViewType refEdgeTangent,
252 const ParamViewType edgeParametrization,
253 const ordinal_type edgeOrdinal ) {
254
255 ordinal_type dim = edgeParametrization.extent(1);
256 for(ordinal_type i = 0; i < dim; ++i) {
257 refEdgeTangent(i) = edgeParametrization(edgeOrdinal, i, 1);
258 }
259 }
260
261 template<typename refFaceTanViewType, typename ParamViewType>
262 KOKKOS_INLINE_FUNCTION
263 static void
264 getReferenceFaceTangent(const refFaceTanViewType refFaceTanU,
265 const refFaceTanViewType refFaceTanV,
266 const ParamViewType faceParametrization,
267 const ordinal_type faceOrdinal) {
268
269 // set refFaceTanU -> C_1(*)
270 // set refFaceTanV -> C_2(*)
271 ordinal_type dim = faceParametrization.extent(1);
272 for(ordinal_type i = 0; i < dim; ++i) {
273 refFaceTanU(i) = faceParametrization(faceOrdinal, i, 1);
274 refFaceTanV(i) = faceParametrization(faceOrdinal, i, 2);
275 }
276 }
277
278 template<typename edgeTangentViewType, typename ParamViewType, typename jacobianViewType>
279 KOKKOS_INLINE_FUNCTION
280 static void
281 getPhysicalEdgeTangent(const edgeTangentViewType edgeTangent, // physDim
282 const ParamViewType edgeParametrization,
283 const jacobianViewType jacobian, // physDim, physDim
284 const ordinal_type edgeOrdinal) {
285 typedef typename ParamViewType::non_const_value_type value_type;
286 const ordinal_type dim = edgeParametrization.extent(1);
287 value_type buf[3];
288 Kokkos::DynRankView<value_type, Kokkos::AnonymousSpace,
289 Kokkos::MemoryUnmanaged> refEdgeTangent(&buf[0], dim);
290
291 getReferenceEdgeTangent(refEdgeTangent, edgeParametrization, edgeOrdinal);
292 Kernels::Serial::matvec_product(edgeTangent, jacobian, refEdgeTangent);
293 }
294
295 template<typename faceTanViewType, typename ParamViewType, typename jacobianViewType>
296 KOKKOS_INLINE_FUNCTION
297 static void
298 getPhysicalFaceTangents( faceTanViewType faceTanU, // physDim
299 faceTanViewType faceTanV, // physDim
300 const ParamViewType faceParametrization,
301 const jacobianViewType jacobian, // physDim, physDim
302 const ordinal_type faceOrdinal) {
303 typedef typename ParamViewType::non_const_value_type value_type;
304 const ordinal_type dim = faceParametrization.extent(1);
305 INTREPID2_TEST_FOR_ABORT(dim != 3,
306 "computing face tangents requires dimension 3.");
307 value_type buf[6];
308 Kokkos::DynRankView<value_type,
309 Kokkos::AnonymousSpace,
310 Kokkos::MemoryUnmanaged> refFaceTanU(&buf[0], dim), refFaceTanV(&buf[3], dim);
311
312 getReferenceFaceTangent(refFaceTanU,
313 refFaceTanV,
314 faceParametrization,
315 faceOrdinal);
316
317 Kernels::Serial::matvec_product_d3(faceTanU, jacobian, refFaceTanU);
318 Kernels::Serial::matvec_product_d3(faceTanV, jacobian, refFaceTanV);
319 }
320
321
322 template<typename faceNormalViewType, typename ParamViewType, typename jacobianViewType>
323 KOKKOS_INLINE_FUNCTION
324 static void
325 getPhysicalFaceNormal(const faceNormalViewType faceNormal, // physDim
326 const ParamViewType faceParametrization,
327 const jacobianViewType jacobian, // physDim, physDim
328 const ordinal_type faceOrdinal) {
329 typedef typename ParamViewType::non_const_value_type value_type;
330 const ordinal_type dim = faceParametrization.extent(1);
331 INTREPID2_TEST_FOR_ABORT(dim != 3,
332 "computing face normal requires dimension 3.");
333 value_type buf[6];
334 Kokkos::DynRankView<value_type,
335 Kokkos::AnonymousSpace,
336 Kokkos::MemoryUnmanaged> faceTanU(&buf[0], dim), faceTanV(&buf[3], dim);
337
338 getPhysicalFaceTangents(faceTanU, faceTanV,
339 faceParametrization,
340 jacobian,
341 faceOrdinal);
342 Kernels::Serial::vector_product_d3(faceNormal, faceTanU, faceTanV);
343 }
344
345 template<typename sideNormalViewType, typename ParamViewType, typename jacobianViewType>
346 KOKKOS_INLINE_FUNCTION
347 static void
348 getPhysicalSideNormal(const sideNormalViewType sideNormal, // physDim
349 const ParamViewType sideParametrization,
350 const jacobianViewType jacobian, // physDim, physDim
351 const ordinal_type sideOrdinal) {
352 const ordinal_type dim = sideParametrization.extent(1);
353 typedef typename ParamViewType::non_const_value_type value_type;
354 switch (dim) {
355 case 2: {
356 value_type buf[3];
357 Kokkos::DynRankView<value_type,
358 Kokkos::AnonymousSpace,
359 Kokkos::MemoryUnmanaged> edgeTangent(&buf[0], dim);
360 getPhysicalEdgeTangent(edgeTangent, sideParametrization, jacobian, sideOrdinal);
361 sideNormal(0) = edgeTangent(1);
362 sideNormal(1) = -edgeTangent(0);
363 break;
364 }
365 case 3: {
366 getPhysicalFaceNormal(sideNormal, sideParametrization, jacobian, sideOrdinal);
367 break;
368 }
369 default: {
370 INTREPID2_TEST_FOR_ABORT(true, "cell dimension is out of range.");
371 break;
372 }
373 }
374 }
375 };
376};
377}
378}
379
380#endif
381
Header file for the classes: Intrepid2::RefSubcellParametrization, Intrepid2::RefCellNodes,...
Header file for small functions used in Intrepid2.
Contains definitions of custom data types in Intrepid2.
Header function for Intrepid2::Util class and other utility functions.
static constexpr ordinal_type MaxNewton
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
static KOKKOS_INLINE_FUNCTION bool mapToReferenceFrame(const refPointViewType &refPoint, const phyPointViewType &physPoint, const nodeViewType &nodes, const double tol=tolerence())
Computation of , the inverse of the reference-to-physical frame map.