48 #ifndef __INTREPID2_CELLTOOLS_SERIAL_HPP__
49 #define __INTREPID2_CELLTOOLS_SERIAL_HPP__
51 #include "Intrepid2_ConfigDefs.hpp"
53 #include "Shards_CellTopology.hpp"
79 template<
typename jacobianViewType,
80 typename basisGradViewType,
81 typename nodeViewType>
82 KOKKOS_INLINE_FUNCTION
84 computeJacobian(
const jacobianViewType &jacobian,
85 const basisGradViewType &grads,
86 const nodeViewType &nodes) {
87 const auto numNodes = nodes.extent(0);
89 const auto physDim = jacobian.extent(0);
90 const auto refDim = jacobian.extent(1);
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.");
96 Kernels::Serial::gemm_trans_notrans(1.0, nodes, grads, 0.0, jacobian);
104 template<
typename PointViewType,
105 typename basisValViewType,
106 typename nodeViewType>
107 KOKKOS_INLINE_FUNCTION
109 mapToPhysicalFrame(
const PointViewType &point,
110 const basisValViewType &vals,
111 const nodeViewType &nodes) {
112 const auto numNodes = vals.extent(0);
113 const auto physDim = point.extent(0);
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.");
118 Kernels::Serial::gemv_trans(1.0, nodes, vals, 0.0, point);
146 template<
typename implBasisType,
147 typename refPointViewType,
148 typename phyPointViewType,
149 typename nodeViewType>
150 KOKKOS_INLINE_FUNCTION
153 const phyPointViewType &physPoint,
154 const nodeViewType &nodes,
155 const double tol = tolerence()) {
156 const ordinal_type refDim = refPoint.extent(0);
157 const ordinal_type physDim = physPoint.extent(0);
158 const ordinal_type numNodes = nodes.extent(0);
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");
164 typedef typename refPointViewType::non_const_value_type value_type;
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;
173 Kokkos::DynRankView<value_type,
174 Kokkos::AnonymousSpace,
175 Kokkos::MemoryUnmanaged> vals(ptr, numNodes); ptr += numNodes;
177 Kokkos::DynRankView<value_type,
178 Kokkos::AnonymousSpace,
179 Kokkos::MemoryUnmanaged> jac(ptr, physDim, refDim); ptr += physDim*refDim;
181 Kokkos::DynRankView<value_type,
182 Kokkos::AnonymousSpace,
183 Kokkos::MemoryUnmanaged> metric(ptr, refDim, refDim); ptr += refDim*refDim;
185 Kokkos::DynRankView<value_type,
186 Kokkos::AnonymousSpace,
187 Kokkos::MemoryUnmanaged> invMetric(ptr, refDim, refDim); ptr += refDim*refDim;
189 Kokkos::DynRankView<value_type,
190 Kokkos::AnonymousSpace,
191 Kokkos::MemoryUnmanaged> invDf(ptr, refDim, physDim); ptr += refDim*physDim;
193 Kokkos::DynRankView<value_type,
194 Kokkos::AnonymousSpace,
195 Kokkos::MemoryUnmanaged> tmpPhysPoint(ptr, physDim); ptr += physDim;
197 Kokkos::DynRankView<value_type,
198 Kokkos::AnonymousSpace,
199 Kokkos::MemoryUnmanaged> oldRefPoint(ptr, refDim); ptr += refDim;
202 for (ordinal_type j=0;j<refDim;++j) oldRefPoint(j) = 0;
204 double xphyNorm = Kernels::Serial::norm(physPoint, NORM_ONE);
206 bool converged =
false;
210 mapToPhysicalFrame(tmpPhysPoint, vals, nodes);
211 Kernels::Serial::z_is_axby(tmpPhysPoint, 1.0, physPoint, -1.0, tmpPhysPoint);
212 double residualNorm = Kernels::Serial::norm(tmpPhysPoint, NORM_ONE);
213 if (residualNorm < tol*xphyNorm) {
220 CellTools::Serial::computeJacobian(jac, grads, nodes);
222 if(physDim == refDim) {
223 Kernels::Serial::inverse(invDf, jac);
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);
231 Kernels::Serial::gemv_notrans(1.0, invDf, tmpPhysPoint, 0.0, refPoint);
232 Kernels::Serial::z_is_axby(refPoint, 1.0, oldRefPoint, 1.0, refPoint);
235 Kernels::Serial::z_is_axby(oldRefPoint, 1.0, oldRefPoint, -1.0, refPoint);
237 double err = Kernels::Serial::norm(oldRefPoint, NORM_ONE);
243 Kernels::Serial::copy(oldRefPoint, refPoint);
248 template<
typename refEdgeTanViewType,
typename ParamViewType>
249 KOKKOS_INLINE_FUNCTION
251 getReferenceEdgeTangent(
const refEdgeTanViewType refEdgeTangent,
252 const ParamViewType edgeParametrization,
253 const ordinal_type edgeOrdinal ) {
255 ordinal_type dim = edgeParametrization.extent(1);
256 for(ordinal_type i = 0; i < dim; ++i) {
257 refEdgeTangent(i) = edgeParametrization(edgeOrdinal, i, 1);
261 template<
typename refFaceTanViewType,
typename ParamViewType>
262 KOKKOS_INLINE_FUNCTION
264 getReferenceFaceTangent(
const refFaceTanViewType refFaceTanU,
265 const refFaceTanViewType refFaceTanV,
266 const ParamViewType faceParametrization,
267 const ordinal_type faceOrdinal) {
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);
278 template<
typename edgeTangentViewType,
typename ParamViewType,
typename jacobianViewType>
279 KOKKOS_INLINE_FUNCTION
281 getPhysicalEdgeTangent(
const edgeTangentViewType edgeTangent,
282 const ParamViewType edgeParametrization,
283 const jacobianViewType jacobian,
284 const ordinal_type edgeOrdinal) {
285 typedef typename ParamViewType::non_const_value_type value_type;
286 const ordinal_type dim = edgeParametrization.extent(1);
288 Kokkos::DynRankView<value_type, Kokkos::AnonymousSpace,
289 Kokkos::MemoryUnmanaged> refEdgeTangent(&buf[0], dim);
291 getReferenceEdgeTangent(refEdgeTangent, edgeParametrization, edgeOrdinal);
292 Kernels::Serial::matvec_product(edgeTangent, jacobian, refEdgeTangent);
295 template<
typename faceTanViewType,
typename ParamViewType,
typename jacobianViewType>
296 KOKKOS_INLINE_FUNCTION
298 getPhysicalFaceTangents( faceTanViewType faceTanU,
299 faceTanViewType faceTanV,
300 const ParamViewType faceParametrization,
301 const jacobianViewType jacobian,
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.");
308 Kokkos::DynRankView<value_type,
309 Kokkos::AnonymousSpace,
310 Kokkos::MemoryUnmanaged> refFaceTanU(&buf[0], dim), refFaceTanV(&buf[3], dim);
312 getReferenceFaceTangent(refFaceTanU,
317 Kernels::Serial::matvec_product_d3(faceTanU, jacobian, refFaceTanU);
318 Kernels::Serial::matvec_product_d3(faceTanV, jacobian, refFaceTanV);
322 template<
typename faceNormalViewType,
typename ParamViewType,
typename jacobianViewType>
323 KOKKOS_INLINE_FUNCTION
325 getPhysicalFaceNormal(
const faceNormalViewType faceNormal,
326 const ParamViewType faceParametrization,
327 const jacobianViewType jacobian,
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.");
334 Kokkos::DynRankView<value_type,
335 Kokkos::AnonymousSpace,
336 Kokkos::MemoryUnmanaged> faceTanU(&buf[0], dim), faceTanV(&buf[3], dim);
338 getPhysicalFaceTangents(faceTanU, faceTanV,
342 Kernels::Serial::vector_product_d3(faceNormal, faceTanU, faceTanV);
345 template<
typename s
ideNormalViewType,
typename ParamViewType,
typename jacobianViewType>
346 KOKKOS_INLINE_FUNCTION
348 getPhysicalSideNormal(
const sideNormalViewType sideNormal,
349 const ParamViewType sideParametrization,
350 const jacobianViewType jacobian,
351 const ordinal_type sideOrdinal) {
352 const ordinal_type dim = sideParametrization.extent(1);
353 typedef typename ParamViewType::non_const_value_type value_type;
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);
366 getPhysicalFaceNormal(sideNormal, sideParametrization, jacobian, sideOrdinal);
370 INTREPID2_TEST_FOR_ABORT(
true,
"cell dimension is out of range.");
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...
Header file for the classes: Intrepid2::RefSubcellParametrization, Intrepid2::RefCellNodes, Intrepid2::RefCellCenter.
Contains definitions of custom data types in Intrepid2.
Header file for small functions used in Intrepid2.