Intrepid2
Intrepid2_CellToolsDefPhysToRef.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
43
49#ifndef __INTREPID2_CELLTOOLS_DEF_PHYS_TO_REF_HPP__
50#define __INTREPID2_CELLTOOLS_DEF_PHYS_TO_REF_HPP__
51
52// disable clang warnings
53#if defined (__clang__) && !defined (__INTEL_COMPILER)
54#pragma clang system_header
55#endif
56
57namespace Intrepid2 {
58
59
60 //============================================================================================//
61 // //
62 // Reference-to-physical frame mapping and its inverse //
63 // //
64 //============================================================================================//
65
66 template<typename DeviceType>
67 template<typename refPointValueType, class ...refPointProperties,
68 typename physPointValueType, class ...physPointProperties,
69 typename worksetCellValueType, class ...worksetCellProperties>
70 void
72 mapToReferenceFrame( Kokkos::DynRankView<refPointValueType,refPointProperties...> refPoints,
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;
83
84 static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceFrame(..): input/output views' memory spaces are not compatible with DeviceType");
85
86#ifdef HAVE_INTREPID2_DEBUG
87 CellTools_mapToReferenceFrameArgs(refPoints, physPoints, worksetCell, cellTopo);
88#endif
89 using deviceType = typename decltype(refPoints)::device_type;
90
92 typedef Kokkos::DynRankView<typename ScalarTraits<refPointValueType>::scalar_type,deviceType> refPointViewSpType;
93
94 const auto spaceDim = cellTopo.getDimension();
95 refPointViewSpType
96 cellCenter("CellTools::mapToReferenceFrame::cellCenter", spaceDim);
97 getReferenceCellCenter(cellCenter, cellTopo);
98
99 // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array. Requires (C,P,D) initial guess.
100 const auto numCells = worksetCell.extent(0);
101 const auto numPoints = physPoints.extent(1);
102
103 // init guess is created locally and non fad whatever refpoints type is
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 );
108 //refPointViewSpType initGuess("CellTools::mapToReferenceFrame::initGuess", numCells, numPoints, spaceDim);
109 rst::clone(initGuess, cellCenter);
110
111 mapToReferenceFrameInitGuess(refPoints, initGuess, physPoints, worksetCell, cellTopo);
112 }
113
114
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>
121 void
123 mapToReferenceFrameInitGuess( Kokkos::DynRankView<refPointValueType,refPointProperties...> refPoints,
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());
131
132#endif
133
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;
143
144 static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceFrameInitGuess(..): input/output views' memory spaces are not compatible with DeviceType");
145
146
147 const auto cellTopo = basis->getBaseCellTopology();
148 const auto spaceDim = cellTopo.getDimension();
149
150 // Default: map (C,P,D) array of physical pt. sets to (C,P,D) array.
151 // Requires (C,P,D) temp arrays and (C,P,D,D) Jacobians.
152 const auto numCells = worksetCell.extent(0);
153 const auto numPoints = physPoints.extent(1);
154
155 using rst = RealSpaceTools<DeviceType>;
156 const auto tol = tolerence();
157
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 >;
161
162 // Temp arrays for Newton iterates and Jacobians. Resize according to rank of ref. point array
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);
165
166 // deep copy may not work with FAD but this is right thing to do as it can move data between devices
167 Kokkos::deep_copy(xOld, initGuess);
168
169 // jacobian should select fad dimension between xOld and worksetCell as they are input; no front interface yet
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);
174
175 using errorViewType = Kokkos::DynRankView<typename ScalarTraits<refPointValueType>::scalar_type, DeviceType>;
176 errorViewType
177 xScalarTmp ("CellTools::mapToReferenceFrameInitGuess::xScalarTmp", numCells, numPoints, spaceDim),
178 errorPointwise("CellTools::mapToReferenceFrameInitGuess::errorPointwise", numCells, numPoints),
179 errorCellwise ("CellTools::mapToReferenceFrameInitGuess::errorCellwise", numCells);
180
181 // Newton method to solve the equation F(refPoints) - physPoints = 0:
182 // refPoints = xOld - DF^{-1}(xOld)*(F(xOld) - physPoints) = xOld + DF^{-1}(xOld)*(physPoints - F(xOld))
183 for (ordinal_type iter=0;iter<Parameters::MaxNewton;++iter) {
184
185 // Jacobians at the old iterates and their inverses.
186 setJacobian(jacobian, xOld, worksetCell, basis);
187 setJacobianInv(jacobianInv, jacobian);
188
189 // The Newton step.
190 mapToPhysicalFrame(xTmp, xOld, worksetCell, basis); // xTmp <- F(xOld)
191 rst::subtract(xTmp, physPoints, xTmp); // xTmp <- physPoints - F(xOld)
192 rst::matvec(refPoints, jacobianInv, xTmp); // refPoints <- DF^{-1}( physPoints - F(xOld) )
193 rst::add(refPoints, xOld); // refPoints <- DF^{-1}( physPoints - F(xOld) ) + xOld
194
195 // l2 error (Euclidean distance) between old and new iterates: |xOld - xNew|
196 rst::subtract(xTmp, xOld, refPoints);
197
198 // extract values
199 rst::extractScalarValues(xScalarTmp, xTmp);
200 rst::vectorNorm(errorPointwise, xScalarTmp, NORM_TWO);
201
202 // Average L2 error for a multiple sets of physical points: error is rank-2 (C,P) array
203 rst::vectorNorm(errorCellwise, errorPointwise, NORM_ONE);
204
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);
208
209 // Stopping criterion:
210 if (errorTotal < tol)
211 break;
212
213 // initialize next Newton step ( this is not device friendly )
214 Kokkos::deep_copy(xOld, refPoints);
215 }
216 }
217
218}
219
220#endif
void CellTools_mapToReferenceFrameArgs(const refPointViewType refPoints, const physPointViewType physPoints, const worksetCellViewType worksetCell, const shards::CellTopology cellTopo)
Validates arguments to Intrepid2::CellTools::mapToReferenceFrame with default initial guess.
static void mapToReferenceFrame(Kokkos::DynRankView< refPointValueType, refPointProperties... > refPoints, const Kokkos::DynRankView< physPointValueType, physPointProperties... > physPoints, const Kokkos::DynRankView< worksetCellValueType, worksetCellProperties... > worksetCell, const shards::CellTopology cellTopo)
Computes , the inverse of the reference-to-physical frame map using a default initial guess.
static void mapToReferenceFrameInitGuess(Kokkos::DynRankView< refPointValueType, refPointProperties... > refPoints, const Kokkos::DynRankView< initGuessValueType, initGuessProperties... > initGuess, const Kokkos::DynRankView< physPointValueType, physPointProperties... > physPoints, const Kokkos::DynRankView< worksetCellValueType, worksetCellProperties... > worksetCell, const HGradBasisPtrType basis)
Computation of , the inverse of the reference-to-physical frame map using user-supplied initial guess...
static constexpr ordinal_type MaxNewton
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
Implementation of basic linear algebra functionality in Euclidean space.
layout deduction (temporary meta-function)