Intrepid2
Intrepid2_CellTools_Serial.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
15#ifndef __INTREPID2_CELLTOOLS_SERIAL_HPP__
16#define __INTREPID2_CELLTOOLS_SERIAL_HPP__
17
18#include "Intrepid2_ConfigDefs.hpp"
19
20#include "Shards_CellTopology.hpp"
21
22#include "Intrepid2_Types.hpp"
23#include "Intrepid2_Utils.hpp"
24#include "Intrepid2_Kernels.hpp"
25
26#include "Intrepid2_CellData.hpp"
27
28namespace Intrepid2 {
29
30namespace Impl {
31
36class CellTools {
37public:
38
39 struct Serial {
40
41 // output:
42 // jacobian (physDim,refDim) - jacobian matrix evaluated at a single point
43 // input:
44 // grads (numNodes,refDim) - hgrad basis grad values evaluated at a single point (C1/C2 element only)
45 // nodes (numNodes,physDim) - cell element-to-node connectivity
46 template<typename jacobianViewType,
47 typename basisGradViewType,
48 typename nodeViewType>
49 KOKKOS_INLINE_FUNCTION
50 static void
51 computeJacobian(const jacobianViewType &jacobian, // physDim,refDim
52 const basisGradViewType &grads, // numNodes,refDim
53 const nodeViewType &nodes) { // numNodes,physDim
54 const auto numNodes = nodes.extent(0);
55
56 const auto physDim = jacobian.extent(0);
57 const auto refDim = jacobian.extent(1);
58
59 INTREPID2_TEST_FOR_ABORT( numNodes != grads.extent(0), "grad dimension_0 does not match to cardinality.");
60 INTREPID2_TEST_FOR_ABORT(refDim != grads.extent(1), "grad dimension_1 does not match to space dim.");
61 INTREPID2_TEST_FOR_ABORT( physDim != nodes.extent(1), "node dimension_1 does not match to space dim.");
62
63 Kernels::Serial::gemm_trans_notrans(1.0, nodes, grads, 0.0, jacobian);
64 }
65
66 // output:
67 // point (physDim) - mapped physical point
68 // input:
69 // vals (numNodes) - hgrad basis values evaluated at a single point (C1/C2 element only)
70 // nodes (numNodes,physDim) - cell element-to-node connectivity
71 template<typename PointViewType,
72 typename basisValViewType,
73 typename nodeViewType>
74 KOKKOS_INLINE_FUNCTION
75 static void
76 mapToPhysicalFrame(const PointViewType &point, // physDim
77 const basisValViewType &vals, // numNodes
78 const nodeViewType &nodes) { // numNodes,physDim
79 const auto numNodes = vals.extent(0);
80 const auto physDim = point.extent(0);
81
82 INTREPID2_TEST_FOR_ABORT(numNodes != nodes.extent(0), "nodes dimension_0 does not match to vals dimension_0.");
83 INTREPID2_TEST_FOR_ABORT(physDim != nodes.extent(1), "node dimension_1 does not match to space dim.");
84
85 Kernels::Serial::gemv_trans(1.0, nodes, vals, 0.0, point);
86 }
87
113 template<typename implBasisType,
114 typename refPointViewType,
115 typename phyPointViewType,
116 typename nodeViewType>
117 KOKKOS_INLINE_FUNCTION
118 static bool
119 mapToReferenceFrame(const refPointViewType &refPoint, // refDim
120 const phyPointViewType &physPoint, // physDim
121 const nodeViewType &nodes,
122 const double tol = tolerence()) { // numNodes,physDim
123 const ordinal_type refDim = refPoint.extent(0);
124 const ordinal_type physDim = physPoint.extent(0);
125 const ordinal_type numNodes = nodes.extent(0);
126
127 INTREPID2_TEST_FOR_ABORT(refDim > physDim, "the dimension of the reference cell is greater than physical cell dimension.");
128 INTREPID2_TEST_FOR_ABORT(physDim != static_cast<ordinal_type>(nodes.extent(1)), "physPoint dimension_0 does not match to space dim.");
129 INTREPID2_TEST_FOR_ABORT(numNodes > 27, "function hard-coded to support at most mappings with 27 Dofs");
130
131 typedef typename refPointViewType::non_const_value_type value_type;
132
133 // I want to use view instead of dynrankview
134 // NMAX = 27, MAXDIM = 3
135 value_type buf[27*3 + 27 + 9 + 9 + 9 + 9 + 3 + 3] = {}, *ptr = &buf[0];
136 Kokkos::DynRankView<value_type,
137 Kokkos::AnonymousSpace,
138 Kokkos::MemoryUnmanaged> grads(ptr, numNodes, refDim); ptr += numNodes*refDim;
139
140 Kokkos::DynRankView<value_type,
141 Kokkos::AnonymousSpace,
142 Kokkos::MemoryUnmanaged> vals(ptr, numNodes); ptr += numNodes;
143
144 Kokkos::DynRankView<value_type,
145 Kokkos::AnonymousSpace,
146 Kokkos::MemoryUnmanaged> jac(ptr, physDim, refDim); ptr += physDim*refDim;
147
148 Kokkos::DynRankView<value_type,
149 Kokkos::AnonymousSpace,
150 Kokkos::MemoryUnmanaged> metric(ptr, refDim, refDim); ptr += refDim*refDim;
151
152 Kokkos::DynRankView<value_type,
153 Kokkos::AnonymousSpace,
154 Kokkos::MemoryUnmanaged> invMetric(ptr, refDim, refDim); ptr += refDim*refDim;
155
156 Kokkos::DynRankView<value_type,
157 Kokkos::AnonymousSpace,
158 Kokkos::MemoryUnmanaged> invDf(ptr, refDim, physDim); ptr += refDim*physDim;
159
160 Kokkos::DynRankView<value_type,
161 Kokkos::AnonymousSpace,
162 Kokkos::MemoryUnmanaged> tmpPhysPoint(ptr, physDim); ptr += physDim;
163
164 Kokkos::DynRankView<value_type,
165 Kokkos::AnonymousSpace,
166 Kokkos::MemoryUnmanaged> oldRefPoint(ptr, refDim); ptr += refDim;
167
168 // set initial guess
169 for (ordinal_type j=0;j<refDim;++j) oldRefPoint(j) = 0;
170
171 double xphyNorm = Kernels::Serial::norm(physPoint, NORM_ONE);
172
173 bool converged = false;
174 for (ordinal_type iter=0;iter<Parameters::MaxNewton;++iter) {
175 // xtmp := F(oldRefPoint);
176 implBasisType::template Serial<OPERATOR_VALUE>::getValues(vals, oldRefPoint);
177 mapToPhysicalFrame(tmpPhysPoint, vals, nodes);
178 Kernels::Serial::z_is_axby(tmpPhysPoint, 1.0, physPoint, -1.0, tmpPhysPoint); //residual xtmp := physPoint - F(oldRefPoint);
179 double residualNorm = Kernels::Serial::norm(tmpPhysPoint, NORM_ONE);
180 if (residualNorm < tol*xphyNorm) {
181 converged = true;
182 break;
183 }
184
185 // DF^{-1}
186 implBasisType::template Serial<OPERATOR_GRAD>::getValues(grads, oldRefPoint);
187 CellTools::Serial::computeJacobian(jac, grads, nodes);
188
189 if(physDim == refDim) { //jacobian inverse
190 Kernels::Serial::inverse(invDf, jac);
191 } else { //jacobian is not square, we compute the pseudo-inverse (jac^T jac)^{-1} jac^T
192 Kernels::Serial::gemm_trans_notrans(1.0, jac, jac, 0.0, metric);
193 Kernels::Serial::inverse(invMetric, metric);
194 Kernels::Serial::gemm_notrans_trans(1.0, invMetric, jac, 0.0, invDf);
195 }
196
197 // Newton
198 Kernels::Serial::gemv_notrans(1.0, invDf, tmpPhysPoint, 0.0, refPoint); // refPoint := DF^{-1}( physPoint - F(oldRefPoint))
199 Kernels::Serial::z_is_axby(refPoint, 1.0, oldRefPoint, 1.0, refPoint); // refPoint += oldRefPoint
200
201 // temporarily overwriting oldRefPoint with oldRefPoint-refPoint
202 Kernels::Serial::z_is_axby(oldRefPoint, 1.0, oldRefPoint, -1.0, refPoint);
203
204 double err = Kernels::Serial::norm(oldRefPoint, NORM_ONE);
205 if (err < tol) {
206 converged = true;
207 break;
208 }
209
210 Kernels::Serial::copy(oldRefPoint, refPoint);
211 }
212 return converged;
213 }
214
215 template<typename refEdgeTanViewType, typename ParamViewType>
216 KOKKOS_INLINE_FUNCTION
217 static void
218 getReferenceEdgeTangent(const refEdgeTanViewType refEdgeTangent,
219 const ParamViewType edgeParametrization,
220 const ordinal_type edgeOrdinal ) {
221
222 ordinal_type dim = edgeParametrization.extent(1);
223 for(ordinal_type i = 0; i < dim; ++i) {
224 refEdgeTangent(i) = edgeParametrization(edgeOrdinal, i, 1);
225 }
226 }
227
228 template<typename refFaceTanViewType, typename ParamViewType>
229 KOKKOS_INLINE_FUNCTION
230 static void
231 getReferenceFaceTangent(const refFaceTanViewType refFaceTanU,
232 const refFaceTanViewType refFaceTanV,
233 const ParamViewType faceParametrization,
234 const ordinal_type faceOrdinal) {
235
236 // set refFaceTanU -> C_1(*)
237 // set refFaceTanV -> C_2(*)
238 ordinal_type dim = faceParametrization.extent(1);
239 for(ordinal_type i = 0; i < dim; ++i) {
240 refFaceTanU(i) = faceParametrization(faceOrdinal, i, 1);
241 refFaceTanV(i) = faceParametrization(faceOrdinal, i, 2);
242 }
243 }
244
245 template<typename edgeTangentViewType, typename ParamViewType, typename jacobianViewType>
246 KOKKOS_INLINE_FUNCTION
247 static void
248 getPhysicalEdgeTangent(const edgeTangentViewType edgeTangent, // physDim
249 const ParamViewType edgeParametrization,
250 const jacobianViewType jacobian, // physDim, physDim
251 const ordinal_type edgeOrdinal) {
252 typedef typename ParamViewType::non_const_value_type value_type;
253 const ordinal_type dim = edgeParametrization.extent(1);
254 value_type buf[3];
255 Kokkos::DynRankView<value_type, Kokkos::AnonymousSpace,
256 Kokkos::MemoryUnmanaged> refEdgeTangent(&buf[0], dim);
257
258 getReferenceEdgeTangent(refEdgeTangent, edgeParametrization, edgeOrdinal);
259 Kernels::Serial::matvec_product(edgeTangent, jacobian, refEdgeTangent);
260 }
261
262 template<typename faceTanViewType, typename ParamViewType, typename jacobianViewType>
263 KOKKOS_INLINE_FUNCTION
264 static void
265 getPhysicalFaceTangents( faceTanViewType faceTanU, // physDim
266 faceTanViewType faceTanV, // physDim
267 const ParamViewType faceParametrization,
268 const jacobianViewType jacobian, // physDim, physDim
269 const ordinal_type faceOrdinal) {
270 typedef typename ParamViewType::non_const_value_type value_type;
271 const ordinal_type dim = faceParametrization.extent(1);
272 INTREPID2_TEST_FOR_ABORT(dim != 3,
273 "computing face tangents requires dimension 3.");
274 value_type buf[6];
275 Kokkos::DynRankView<value_type,
276 Kokkos::AnonymousSpace,
277 Kokkos::MemoryUnmanaged> refFaceTanU(&buf[0], dim), refFaceTanV(&buf[3], dim);
278
279 getReferenceFaceTangent(refFaceTanU,
280 refFaceTanV,
281 faceParametrization,
282 faceOrdinal);
283
284 Kernels::Serial::matvec_product_d3(faceTanU, jacobian, refFaceTanU);
285 Kernels::Serial::matvec_product_d3(faceTanV, jacobian, refFaceTanV);
286 }
287
288
289 template<typename faceNormalViewType, typename ParamViewType, typename jacobianViewType>
290 KOKKOS_INLINE_FUNCTION
291 static void
292 getPhysicalFaceNormal(const faceNormalViewType faceNormal, // physDim
293 const ParamViewType faceParametrization,
294 const jacobianViewType jacobian, // physDim, physDim
295 const ordinal_type faceOrdinal) {
296 typedef typename ParamViewType::non_const_value_type value_type;
297 const ordinal_type dim = faceParametrization.extent(1);
298 INTREPID2_TEST_FOR_ABORT(dim != 3,
299 "computing face normal requires dimension 3.");
300 value_type buf[6];
301 Kokkos::DynRankView<value_type,
302 Kokkos::AnonymousSpace,
303 Kokkos::MemoryUnmanaged> faceTanU(&buf[0], dim), faceTanV(&buf[3], dim);
304
305 getPhysicalFaceTangents(faceTanU, faceTanV,
306 faceParametrization,
307 jacobian,
308 faceOrdinal);
309 Kernels::Serial::vector_product_d3(faceNormal, faceTanU, faceTanV);
310 }
311
312 template<typename sideNormalViewType, typename ParamViewType, typename jacobianViewType>
313 KOKKOS_INLINE_FUNCTION
314 static void
315 getPhysicalSideNormal(const sideNormalViewType sideNormal, // physDim
316 const ParamViewType sideParametrization,
317 const jacobianViewType jacobian, // physDim, physDim
318 const ordinal_type sideOrdinal) {
319 const ordinal_type dim = sideParametrization.extent(1);
320 typedef typename ParamViewType::non_const_value_type value_type;
321 switch (dim) {
322 case 2: {
323 value_type buf[3];
324 Kokkos::DynRankView<value_type,
325 Kokkos::AnonymousSpace,
326 Kokkos::MemoryUnmanaged> edgeTangent(&buf[0], dim);
327 getPhysicalEdgeTangent(edgeTangent, sideParametrization, jacobian, sideOrdinal);
328 sideNormal(0) = edgeTangent(1);
329 sideNormal(1) = -edgeTangent(0);
330 break;
331 }
332 case 3: {
333 getPhysicalFaceNormal(sideNormal, sideParametrization, jacobian, sideOrdinal);
334 break;
335 }
336 default: {
337 INTREPID2_TEST_FOR_ABORT(true, "cell dimension is out of range.");
338 break;
339 }
340 }
341 }
342 };
343};
344}
345}
346
347#endif
348
Intrepid2_CellData.hpp
Header file for the classes: Intrepid2::RefSubcellParametrization, Intrepid2::RefCellNodes,...
Intrepid2_Kernels.hpp
Header file for small functions used in Intrepid2.
Intrepid2_Types.hpp
Contains definitions of custom data types in Intrepid2.
Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
Intrepid2::Parameters::MaxNewton
static constexpr ordinal_type MaxNewton
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
Definition Intrepid2_Types.hpp:116
Intrepid2::Impl::CellTools::Serial::mapToReferenceFrame
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.
Definition Intrepid2_CellTools_Serial.hpp:119
Generated by doxygen 1.9.8