49 #ifndef __INTREPID2_CELLTOOLS_DEF_PHYS_TO_REF_HPP__
50 #define __INTREPID2_CELLTOOLS_DEF_PHYS_TO_REF_HPP__
53 #if defined (__clang__) && !defined (__INTEL_COMPILER)
54 #pragma clang system_header
66 template<
typename DeviceType>
67 template<
typename refPointValueType,
class ...refPointProperties,
68 typename physPointValueType,
class ...physPointProperties,
69 typename worksetCellValueType,
class ...worksetCellProperties>
73 const Kokkos::DynRankView<physPointValueType,physPointProperties...> physPoints,
74 const Kokkos::DynRankView<worksetCellValueType,worksetCellProperties...> worksetCell,
75 const shards::CellTopology cellTopo ) {
76 constexpr
bool are_accessible =
77 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
78 typename decltype(refPoints)::memory_space>::accessible &&
79 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
80 typename decltype(physPoints)::memory_space>::accessible &&
81 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
82 typename decltype(worksetCell)::memory_space>::accessible;
84 static_assert(are_accessible,
"CellTools<DeviceType>::mapToReferenceFrame(..): input/output views' memory spaces are not compatible with DeviceType");
86 #ifdef HAVE_INTREPID2_DEBUG
87 CellTools_mapToReferenceFrameArgs(refPoints, physPoints, worksetCell, cellTopo);
89 using deviceType =
typename decltype(refPoints)::device_type;
92 typedef Kokkos::DynRankView<typename ScalarTraits<refPointValueType>::scalar_type,deviceType> refPointViewSpType;
94 const auto spaceDim = cellTopo.getDimension();
96 cellCenter(
"CellTools::mapToReferenceFrame::cellCenter", spaceDim);
97 getReferenceCellCenter(cellCenter, cellTopo);
100 const auto numCells = worksetCell.extent(0);
101 const auto numPoints = physPoints.extent(1);
104 using result_layout =
typename DeduceLayout< decltype(refPoints) >::result_layout;
105 auto vcprop = Kokkos::common_view_alloc_prop(refPoints);
106 using common_value_type =
typename decltype(vcprop)::value_type;
107 Kokkos::DynRankView< common_value_type, result_layout, deviceType > initGuess ( Kokkos::view_alloc(
"CellTools::mapToReferenceFrame::initGuess", vcprop), numCells, numPoints, spaceDim );
109 rst::clone(initGuess, cellCenter);
111 mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, worksetCell, cellTopo);
115 template<
typename DeviceType>
116 template<
typename refPointValueType,
class ...refPointProperties,
117 typename initGuessValueType,
class ...initGuessProperties,
118 typename physPointValueType,
class ...physPointProperties,
119 typename worksetCellValueType,
class ...worksetCellProperties,
120 typename HGradBasisPtrType>
124 const Kokkos::DynRankView<initGuessValueType,initGuessProperties...> initGuess,
125 const Kokkos::DynRankView<physPointValueType,physPointProperties...> physPoints,
126 const Kokkos::DynRankView<worksetCellValueType,worksetCellProperties...> worksetCell,
127 const HGradBasisPtrType basis ) {
128 #ifdef HAVE_INTREPID2_DEBUG
129 CellTools_mapToReferenceFrameInitGuessArgs(refPoints, initGuess, physPoints, worksetCell,
130 basis->getBaseCellTopology());
134 constexpr
bool are_accessible =
135 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
136 typename decltype(refPoints)::memory_space>::accessible &&
137 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
138 typename decltype(initGuess)::memory_space>::accessible &&
139 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
140 typename decltype(physPoints)::memory_space>::accessible &&
141 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
142 typename decltype(worksetCell)::memory_space>::accessible;
144 static_assert(are_accessible,
"CellTools<DeviceType>::mapToReferenceFrameInitGuess(..): input/output views' memory spaces are not compatible with DeviceType");
147 const auto cellTopo = basis->getBaseCellTopology();
148 const auto spaceDim = cellTopo.getDimension();
152 const auto numCells = worksetCell.extent(0);
153 const auto numPoints = physPoints.extent(1);
156 const auto tol = tolerence();
158 using result_layout =
typename DeduceLayout< decltype(refPoints) >::result_layout;
159 auto vcprop = Kokkos::common_view_alloc_prop(refPoints);
160 using viewType = Kokkos::DynRankView<typename decltype(vcprop)::value_type, result_layout, DeviceType >;
163 viewType xOld(Kokkos::view_alloc(
"CellTools::mapToReferenceFrameInitGuess::xOld", vcprop), numCells, numPoints, spaceDim);
164 viewType xTmp(Kokkos::view_alloc(
"CellTools::mapToReferenceFrameInitGuess::xTmp", vcprop), numCells, numPoints, spaceDim);
167 Kokkos::deep_copy(xOld, initGuess);
170 auto vcpropJ = Kokkos::common_view_alloc_prop(refPoints, worksetCell);
171 using viewTypeJ = Kokkos::DynRankView<typename decltype(vcpropJ)::value_type, result_layout, DeviceType >;
172 viewTypeJ jacobian(Kokkos::view_alloc(
"CellTools::mapToReferenceFrameInitGuess::jacobian", vcpropJ), numCells, numPoints, spaceDim, spaceDim);
173 viewTypeJ jacobianInv(Kokkos::view_alloc(
"CellTools::mapToReferenceFrameInitGuess::jacobianInv", vcpropJ), numCells, numPoints, spaceDim, spaceDim);
175 using errorViewType = Kokkos::DynRankView<typename ScalarTraits<refPointValueType>::scalar_type, DeviceType>;
177 xScalarTmp (
"CellTools::mapToReferenceFrameInitGuess::xScalarTmp", numCells, numPoints, spaceDim),
178 errorPointwise(
"CellTools::mapToReferenceFrameInitGuess::errorPointwise", numCells, numPoints),
179 errorCellwise (
"CellTools::mapToReferenceFrameInitGuess::errorCellwise", numCells);
186 setJacobian(jacobian, xOld, worksetCell, basis);
187 setJacobianInv(jacobianInv, jacobian);
190 mapToPhysicalFrame(xTmp, xOld, worksetCell, basis);
191 rst::subtract(xTmp, physPoints, xTmp);
192 rst::matvec(refPoints, jacobianInv, xTmp);
193 rst::add(refPoints, xOld);
196 rst::subtract(xTmp, xOld, refPoints);
199 rst::extractScalarValues(xScalarTmp, xTmp);
200 rst::vectorNorm(errorPointwise, xScalarTmp, NORM_TWO);
203 rst::vectorNorm(errorCellwise, errorPointwise, NORM_ONE);
205 auto errorCellwise_h = Kokkos::create_mirror_view(errorCellwise);
206 Kokkos::deep_copy(errorCellwise_h, errorCellwise);
207 const auto errorTotal = rst::Serial::vectorNorm(errorCellwise_h, NORM_ONE);
210 if (errorTotal < tol)
214 Kokkos::deep_copy(xOld, refPoints);
static constexpr ordinal_type MaxNewton
Maximum number of Newton iterations used internally in methods such as computing the action of the in...