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
115 template<typename implBasisType,
116 typename refPointViewType,
117 typename phyPointViewType,
118 typename nodeViewType>
119 KOKKOS_INLINE_FUNCTION
120 static bool
121 mapToReferenceFrame(const refPointViewType &refPoint, // refDim
122 const phyPointViewType &physPoint, // physDim
123 const nodeViewType &nodes,
124 const double tol = tolerance(),
125 const typename phyPointViewType::value_type shellThickness = -1.0) { // numNodes,physDim
126
127 bool isBeamOrShell = (shellThickness > 0);
128
129 //for Beam or Shell elements, refPoint has dimension refDim+1 (e.g., shellQuadrilateral is a 3D topology, but maps according to the QUAD basis functions)
130 const ordinal_type refDim = refPoint.extent(0) - ordinal_type(isBeamOrShell);
131 const ordinal_type physDim = physPoint.extent(0);
132 const ordinal_type numNodes = nodes.extent(0);
133
134 INTREPID2_TEST_FOR_ABORT(refDim > physDim, "the dimension of the reference cell is greater than physical cell dimension.");
135 INTREPID2_TEST_FOR_ABORT(physDim != static_cast<ordinal_type>(nodes.extent(1)), "physPoint dimension_0 does not match to space dim.");
136 INTREPID2_TEST_FOR_ABORT(numNodes > 27, "function hard-coded to support at most mappings with 27 Dofs");
137
138 typedef typename refPointViewType::non_const_value_type value_type;
139
140 // I want to use view instead of dynrankview
141 // NMAX = 27, MAXDIM = 3
142 value_type buf[27*3 + 27 + 9 + 9 + 9 + 9 + 3 + 3] = {}, *ptr = &buf[0];
143 Kokkos::DynRankView<value_type,
144 Kokkos::AnonymousSpace,
145 Kokkos::MemoryUnmanaged> grads(ptr, numNodes, refDim); ptr += numNodes*refDim;
146
147 Kokkos::DynRankView<value_type,
148 Kokkos::AnonymousSpace,
149 Kokkos::MemoryUnmanaged> vals(ptr, numNodes); ptr += numNodes;
150
151 Kokkos::DynRankView<value_type,
152 Kokkos::AnonymousSpace,
153 Kokkos::MemoryUnmanaged> jac(ptr, physDim, refDim); ptr += physDim*refDim;
154
155 Kokkos::DynRankView<value_type,
156 Kokkos::AnonymousSpace,
157 Kokkos::MemoryUnmanaged> metric(ptr, refDim, refDim); ptr += refDim*refDim;
158
159 Kokkos::DynRankView<value_type,
160 Kokkos::AnonymousSpace,
161 Kokkos::MemoryUnmanaged> invMetric(ptr, refDim, refDim); ptr += refDim*refDim;
162
163 Kokkos::DynRankView<value_type,
164 Kokkos::AnonymousSpace,
165 Kokkos::MemoryUnmanaged> invDf(ptr, refDim, physDim); ptr += refDim*physDim;
166
167 Kokkos::DynRankView<value_type,
168 Kokkos::AnonymousSpace,
169 Kokkos::MemoryUnmanaged> tmpPhysPoint(ptr, physDim); ptr += physDim;
170
171 Kokkos::DynRankView<value_type,
172 Kokkos::AnonymousSpace,
173 Kokkos::MemoryUnmanaged> oldRefPoint(ptr, refDim); ptr += refDim;
174
175 double xphyNorm = Kernels::Serial::norm(physPoint, NORM_TWO);
176
177 const auto refDimRange = Kokkos::pair<ordinal_type,ordinal_type>(0, refDim);
178 auto refPt = Kokkos::subview(refPoint,refDimRange); //needed for shell elements
179
180 // set initial guess
181 Kernels::Serial::copy(oldRefPoint, refPt);
182
183 bool converged = false;
184 for (ordinal_type iter=0;iter<Parameters::MaxNewton;++iter) {
185 // xtmp := F(oldRefPoint);
186 implBasisType::template Serial<OPERATOR_VALUE>::getValues(vals, oldRefPoint);
187 mapToPhysicalFrame(tmpPhysPoint, vals, nodes);
188 Kernels::Serial::z_is_axby(tmpPhysPoint, 1.0, physPoint, -1.0, tmpPhysPoint); //residual xtmp := physPoint - F(oldRefPoint);
189 double residualNorm = Kernels::Serial::norm(tmpPhysPoint, NORM_TWO);
190 if (residualNorm < tol*xphyNorm) {
191 converged = true;
192 break;
193 }
194
195 // DF^{-1}
196 implBasisType::template Serial<OPERATOR_GRAD>::getValues(grads, oldRefPoint);
197 CellTools::Serial::computeJacobian(jac, grads, nodes);
198
199 if(physDim == refDim) { //jacobian inverse
200 Kernels::Serial::inverse(invDf, jac);
201 } else { //jacobian is not square, we compute the pseudo-inverse (jac^T jac)^{-1} jac^T
202 Kernels::Serial::gemm_trans_notrans(1.0, jac, jac, 0.0, metric);
203 Kernels::Serial::inverse(invMetric, metric);
204 Kernels::Serial::gemm_notrans_trans(1.0, invMetric, jac, 0.0, invDf);
205 }
206
207 // Newton
208 Kernels::Serial::gemv_notrans(1.0, invDf, tmpPhysPoint, 0.0, refPt); // refPoint := DF^{-1}( physPoint - F(oldRefPoint))
209
210 double err = Kernels::Serial::norm(refPt, NORM_TWO); //error on increment
211
212 Kernels::Serial::z_is_axby(refPt, 1.0, oldRefPoint, 1.0, refPt); // refPoint += oldRefPoint
213
214 if (err < tol) {
215 converged = true;
216 break;
217 }
218
219 Kernels::Serial::copy(oldRefPoint, refPt);
220 }
221
222 if(isBeamOrShell) {
223 implBasisType::template Serial<OPERATOR_GRAD>::getValues(grads, refPt);
224 CellTools::Serial::computeJacobian(jac, grads, nodes);
225 if (physDim==3) {
226 auto t1 = Kokkos::subview( jac, Kokkos::ALL(), 0);
227 auto t2 = Kokkos::subview( jac, Kokkos::ALL(), 1);
228 auto n = Kokkos::subview( invDf, 0, Kokkos::ALL());
229 Kernels::Serial::vector_product_d3(n, t1, t2);
230 refPoint(refDim) = (n(0)*tmpPhysPoint(0) + n(1)*tmpPhysPoint(1) + n(2)*tmpPhysPoint(2))*2.0/shellThickness/Kernels::Serial::norm(n, NORM_TWO);
231 } else { //physDim == 2
232 auto t1 = Kokkos::subview( jac, Kokkos::ALL(), 0);
233 refPoint(refDim) = (-t1(1)*tmpPhysPoint(0) + t1(0)*tmpPhysPoint(1))*2.0/shellThickness/Kernels::Serial::norm(t1, NORM_TWO);
234 }
235 }
236
237 return converged;
238 }
239
240 template<typename refEdgeTanViewType, typename ParamViewType>
241 KOKKOS_INLINE_FUNCTION
242 static void
243 getReferenceEdgeTangent(const refEdgeTanViewType refEdgeTangent,
244 const ParamViewType edgeParametrization,
245 const ordinal_type edgeOrdinal ) {
246
247 ordinal_type dim = edgeParametrization.extent(1);
248 for(ordinal_type i = 0; i < dim; ++i) {
249 refEdgeTangent(i) = edgeParametrization(edgeOrdinal, i, 1);
250 }
251 }
252
253 template<typename refFaceTanViewType, typename ParamViewType>
254 KOKKOS_INLINE_FUNCTION
255 static void
256 getReferenceFaceTangent(const refFaceTanViewType refFaceTanU,
257 const refFaceTanViewType refFaceTanV,
258 const ParamViewType faceParametrization,
259 const ordinal_type faceOrdinal) {
260
261 // set refFaceTanU -> C_1(*)
262 // set refFaceTanV -> C_2(*)
263 ordinal_type dim = faceParametrization.extent(1);
264 for(ordinal_type i = 0; i < dim; ++i) {
265 refFaceTanU(i) = faceParametrization(faceOrdinal, i, 1);
266 refFaceTanV(i) = faceParametrization(faceOrdinal, i, 2);
267 }
268 }
269
270 template<typename edgeTangentViewType, typename ParamViewType, typename jacobianViewType>
271 KOKKOS_INLINE_FUNCTION
272 static void
273 getPhysicalEdgeTangent(const edgeTangentViewType edgeTangent, // physDim
274 const ParamViewType edgeParametrization,
275 const jacobianViewType jacobian, // physDim, physDim
276 const ordinal_type edgeOrdinal) {
277 typedef typename ParamViewType::non_const_value_type value_type;
278 const ordinal_type dim = edgeParametrization.extent(1);
279 value_type buf[3];
280 Kokkos::DynRankView<value_type, Kokkos::AnonymousSpace,
281 Kokkos::MemoryUnmanaged> refEdgeTangent(&buf[0], dim);
282
283 getReferenceEdgeTangent(refEdgeTangent, edgeParametrization, edgeOrdinal);
284 Kernels::Serial::matvec_product(edgeTangent, jacobian, refEdgeTangent);
285 }
286
287 template<typename faceTanViewType, typename ParamViewType, typename jacobianViewType>
288 KOKKOS_INLINE_FUNCTION
289 static void
290 getPhysicalFaceTangents( faceTanViewType faceTanU, // physDim
291 faceTanViewType faceTanV, // physDim
292 const ParamViewType faceParametrization,
293 const jacobianViewType jacobian, // physDim, physDim
294 const ordinal_type faceOrdinal) {
295 typedef typename ParamViewType::non_const_value_type value_type;
296 const ordinal_type dim = faceParametrization.extent(1);
297 INTREPID2_TEST_FOR_ABORT(dim != 3,
298 "computing face tangents requires dimension 3.");
299 value_type buf[6];
300 Kokkos::DynRankView<value_type,
301 Kokkos::AnonymousSpace,
302 Kokkos::MemoryUnmanaged> refFaceTanU(&buf[0], dim), refFaceTanV(&buf[3], dim);
303
304 getReferenceFaceTangent(refFaceTanU,
305 refFaceTanV,
306 faceParametrization,
307 faceOrdinal);
308
309 Kernels::Serial::matvec_product_d3(faceTanU, jacobian, refFaceTanU);
310 Kernels::Serial::matvec_product_d3(faceTanV, jacobian, refFaceTanV);
311 }
312
313
314 template<typename faceNormalViewType, typename ParamViewType, typename jacobianViewType>
315 KOKKOS_INLINE_FUNCTION
316 static void
317 getPhysicalFaceNormal(const faceNormalViewType faceNormal, // physDim
318 const ParamViewType faceParametrization,
319 const jacobianViewType jacobian, // physDim, physDim
320 const ordinal_type faceOrdinal) {
321 typedef typename ParamViewType::non_const_value_type value_type;
322 const ordinal_type dim = faceParametrization.extent(1);
323 INTREPID2_TEST_FOR_ABORT(dim != 3,
324 "computing face normal requires dimension 3.");
325 value_type buf[6];
326 Kokkos::DynRankView<value_type,
327 Kokkos::AnonymousSpace,
328 Kokkos::MemoryUnmanaged> faceTanU(&buf[0], dim), faceTanV(&buf[3], dim);
329
330 getPhysicalFaceTangents(faceTanU, faceTanV,
331 faceParametrization,
332 jacobian,
333 faceOrdinal);
334 Kernels::Serial::vector_product_d3(faceNormal, faceTanU, faceTanV);
335 }
336
337 template<typename sideNormalViewType, typename ParamViewType, typename jacobianViewType>
338 KOKKOS_INLINE_FUNCTION
339 static void
340 getPhysicalSideNormal(const sideNormalViewType sideNormal, // physDim
341 const ParamViewType sideParametrization,
342 const jacobianViewType jacobian, // physDim, physDim
343 const ordinal_type sideOrdinal) {
344 const ordinal_type dim = sideParametrization.extent(1);
345 typedef typename ParamViewType::non_const_value_type value_type;
346 switch (dim) {
347 case 2: {
348 value_type buf[3];
349 Kokkos::DynRankView<value_type,
350 Kokkos::AnonymousSpace,
351 Kokkos::MemoryUnmanaged> edgeTangent(&buf[0], dim);
352 getPhysicalEdgeTangent(edgeTangent, sideParametrization, jacobian, sideOrdinal);
353 sideNormal(0) = edgeTangent(1);
354 sideNormal(1) = -edgeTangent(0);
355 break;
356 }
357 case 3: {
358 getPhysicalFaceNormal(sideNormal, sideParametrization, jacobian, sideOrdinal);
359 break;
360 }
361 default: {
362 INTREPID2_TEST_FOR_ABORT(true, "cell dimension is out of range.");
363 break;
364 }
365 }
366 }
367 };
368};
369}
370}
371
372#endif
373
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=tolerance(), const typename phyPointViewType::value_type shellThickness=-1.0)
Computation of , the inverse of the reference-to-physical frame map.
Definition Intrepid2_CellTools_Serial.hpp:121
Generated by doxygen 1.9.8