Intrepid2
Intrepid2_OrientationToolsDefCoeffMatrix_HCURL.hpp
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
62#ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HCURL_HPP__
63#define __INTREPID2_ORIENTATIONTOOLS_DEF_COEFF_MATRIX_HCURL_HPP__
64
68
69// disable clang warnings
70#if defined (__clang__) && !defined (__INTEL_COMPILER)
71#pragma clang system_header
72#endif
73
74namespace Intrepid2 {
75
76namespace Impl {
77
78namespace Debug {
79
80#ifdef HAVE_INTREPID2_DEBUG
81template<typename subcellBasisType,
82typename cellBasisType>
83inline
84void
85check_getCoeffMatrix_HCURL(const subcellBasisType& subcellBasis,
86 const cellBasisType& cellBasis,
87 const ordinal_type subcellId,
88 const ordinal_type subcellOrt) {
89 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
90 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
91
92 const ordinal_type cellDim = cellTopo.getDimension();
93 const ordinal_type subcellDim = subcellTopo.getDimension();
94
95
96
97 INTREPID2_TEST_FOR_EXCEPTION( subcellDim > cellDim,
98 std::logic_error,
99 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
100 "cellDim cannot be smaller than subcellDim.");
101
102 const auto subcellBaseKey = subcellTopo.getBaseKey();
103 const auto cellBaseKey = cellTopo.getBaseKey();
104
105 INTREPID2_TEST_FOR_EXCEPTION( subcellBaseKey != shards::Line<>::key &&
106 subcellBaseKey != shards::Quadrilateral<>::key &&
107 subcellBaseKey != shards::Triangle<>::key,
108 std::logic_error,
109 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
110 "subcellBasis must have line, quad, or triangle topology.");
111
112 INTREPID2_TEST_FOR_EXCEPTION( cellBaseKey != shards::Quadrilateral<>::key &&
113 cellBaseKey != shards::Triangle<>::key &&
114 cellBaseKey != shards::Hexahedron<>::key &&
115 cellBaseKey != shards::Wedge<>::key &&
116 cellBaseKey != shards::Tetrahedron<>::key,
117 std::logic_error,
118 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
119 "cellBasis must have quad, triangle, hexhedron or tetrahedron topology.");
120
121 //
122 // Function space
123 //
124 {
125 const bool cellBasisIsHCURL = cellBasis.getFunctionSpace() == FUNCTION_SPACE_HCURL;
126
127
128 if (cellBasisIsHCURL) {
129 const bool subcellBasisIsHGRAD = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HGRAD;
130 const bool subcellBasisIsHVOL = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HVOL;
131 const bool subcellBasisIsHCURL = subcellBasis.getFunctionSpace() == FUNCTION_SPACE_HCURL;
132 const bool cellIsTri = cellBaseKey == shards::Triangle<>::key;
133 const bool cellIsTet = cellBaseKey == shards::Tetrahedron<>::key;
134 const bool cellIsHex = cellBaseKey == shards::Hexahedron<>::key;
135 const bool cellIsQuad = cellBaseKey == shards::Quadrilateral<>::key;
136
137
138 // edge hcurl is hgrad with gauss legendre points
139 switch (subcellDim) {
140 case 1: {
141 //TODO: Hex, QUAD, TET and TRI element should have the same 1d basis
142 if (cellIsHex || cellIsQuad) {
143 INTREPID2_TEST_FOR_EXCEPTION( !(subcellBasisIsHGRAD||subcellBasisIsHVOL),
144 std::logic_error,
145 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): "
146 "subcellBasis function space (1d) is not consistent to cellBasis.");
147 } else if (cellIsTet || cellIsTri) {
148 INTREPID2_TEST_FOR_EXCEPTION( !subcellBasisIsHVOL,
149 std::logic_error,
150 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): "
151 "subcellBasis function space (1d) is not consistent to cellBasis.");
152 }
153 break;
154 }
155 case 2: {
156 INTREPID2_TEST_FOR_EXCEPTION( !subcellBasisIsHCURL,
157 std::logic_error,
158 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): "
159 "subcellBasis function space (2d) is not consistent to cellBasis.");
160 break;
161 }
162 }
163 }
164 }
165}
166#endif
167}
168
169template<typename OutputViewType,
170typename subcellBasisHostType,
171typename cellBasisHostType>
172inline
173void
175getCoeffMatrix_HCURL(OutputViewType &output,
176 const subcellBasisHostType& subcellBasis,
177 const cellBasisHostType& cellBasis,
178 const ordinal_type subcellId,
179 const ordinal_type subcellOrt) {
180
181#ifdef HAVE_INTREPID2_DEBUG
182 Debug::check_getCoeffMatrix_HCURL(subcellBasis,cellBasis,subcellId,subcellOrt);
183#endif
184
185 using value_type = typename OutputViewType::non_const_value_type;
186 using host_device_type = typename Kokkos::HostSpace::device_type;
187
188 //
189 // Topology
190 //
191 // populate points on a subcell and map to subcell
192 const shards::CellTopology cellTopo = cellBasis.getBaseCellTopology();
193 const shards::CellTopology subcellTopo = subcellBasis.getBaseCellTopology();
194 const ordinal_type cellDim = cellTopo.getDimension();
195 const ordinal_type subcellDim = subcellTopo.getDimension();
196 const auto subcellBaseKey = subcellTopo.getBaseKey();
197 const ordinal_type numCellBasis = cellBasis.getCardinality();
198 const ordinal_type numSubcellBasis = subcellBasis.getCardinality();
199 const ordinal_type ndofSubcell = cellBasis.getDofCount(subcellDim,subcellId);
200
201 // Compute reference points \xi_j and tangents t_j on the subcell
202 // To do so we use the DoF coordinates and DoF coefficients of a Lagrangian HCURL basis
203 // on the subcell spanning the same space as the bases \phi_j
204
205 // Tangents t_j
206 Kokkos::DynRankView<value_type,host_device_type> subcellTangents("subcellTangents", numSubcellBasis, subcellDim);
207 auto degree = subcellBasis.getDegree();
208 BasisPtr<host_device_type, value_type, value_type> basisPtr;
209 if(subcellBaseKey == shards::Line<>::key) {
211 basisPtr->getDofCoeffs(Kokkos::subview(subcellTangents, Kokkos::ALL(),0));
212 } else if (subcellBaseKey == shards::Triangle<>::key) {
214 basisPtr->getDofCoeffs(subcellTangents);
215 } else if (subcellBaseKey == shards::Quadrilateral<>::key) {
217 basisPtr->getDofCoeffs(subcellTangents);
218 }
219
220 // coordinates \xi_j
221 Kokkos::DynRankView<value_type,host_device_type> subcellDofCoords("subcellDofCoords", basisPtr->getCardinality(), subcellDim);
222 basisPtr->getDofCoords(subcellDofCoords);
223 INTREPID2_TEST_FOR_EXCEPTION( basisPtr->getDofCount(subcellDim,0) != ndofSubcell,
224 std::logic_error,
225 ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): " \
226 "the number of basisPtr internal DoFs should equate those of the subcell");
227
228 // restrict \xi_j (and corresponding t_j) to points internal to the HCURL basis
229 Kokkos::DynRankView<value_type,host_device_type> refPtsSubcell("refPtsSubcell", ndofSubcell, subcellDim);
230 Kokkos::DynRankView<value_type,host_device_type> refSubcellTangents("subcellTangents", ndofSubcell, subcellDim);
231 auto tagToOrdinal = basisPtr->getAllDofOrdinal();
232 for (ordinal_type i=0;i<ndofSubcell;++i) {
233 ordinal_type isc = tagToOrdinal(subcellDim, 0, i);
234 for(ordinal_type d=0; d <subcellDim; ++d){
235 refPtsSubcell(i,d) = subcellDofCoords(isc,d);
236 for(ordinal_type k=0; k <subcellDim; ++k)
237 refSubcellTangents(i,d) = subcellTangents(isc,d);
238 }
239 }
240
241 //
242 // Bases evaluation on the reference points
243 //
244
245 // subcellBasisValues = \phi_i (\xi_j)
246 Kokkos::DynRankView<value_type,host_device_type> subCellValues("subCellValues", numSubcellBasis, ndofSubcell, subcellDim);
247 if(subcellDim==1) {
248 auto lineValues = Kokkos::subview(subCellValues, Kokkos::ALL(), Kokkos::ALL(), 0);
249 subcellBasis.getValues(lineValues, refPtsSubcell, OPERATOR_VALUE);
250 } else {
251 subcellBasis.getValues(subCellValues, refPtsSubcell, OPERATOR_VALUE);
252 }
253
254 //
255 // Basis evaluation on the reference points
256 //
257
258 Kokkos::DynRankView<value_type,host_device_type> refPtsCell("refPtsCell", ndofSubcell, cellDim);
259 Kokkos::DynRankView<value_type,host_device_type> trJacobianF("trJacobianF", subcellDim, cellDim );
260
261 if(cellDim == subcellDim) {
262 // refPtsCell = \eta_o (refPtsSubcell)
263 mapToModifiedReference(refPtsCell,refPtsSubcell,subcellBaseKey,subcellOrt);
264
265 //mapping tangents t_j into parent cell, i.e. computing J_F J_\eta t_j
266 Kokkos::DynRankView<value_type,host_device_type> jac("data", subcellDim, subcellDim);
267 getJacobianOfOrientationMap(jac,subcellBaseKey,subcellOrt);
268 for(ordinal_type d=0; d<subcellDim; ++d)
269 for(ordinal_type j=0; j<cellDim; ++j) {
270 trJacobianF(j,d) = jac(d,j);
271 }
272 }
273 else {
274 auto subcellParam = Intrepid2::RefSubcellParametrization<host_device_type>::get(subcellDim, cellTopo.getKey());
275 // refPtsCell = F_s (\eta_o (refPtsSubcell))
276 mapSubcellCoordsToRefCell(refPtsCell,refPtsSubcell, subcellParam, subcellBaseKey, subcellId, subcellOrt);
277
278 //mapping tangents t_j into parent cell, i.e. computing J_F J_\eta t_j
279 OrientationTools::getRefSubcellTangents(trJacobianF, subcellParam, subcellBaseKey, subcellId, subcellOrt);
280 }
281
282 // cellBasisValues = \psi_k(F_s (\eta_o (\xi_j)))
283 Kokkos::DynRankView<value_type,host_device_type> cellBasisValues("cellBasisValues", numCellBasis, ndofSubcell, cellDim);
284 cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
285
286 //
287 // Compute Psi_jk = \psi_k(F_s (\eta_o (\xi_j))) \cdot (J_F J_\eta t_j)
288 // and Phi_ji = \phi_i (\xi_j) \ctot t_j, and solve
289 // Psi A^T = Phi
290 //
291
292 // construct Psi and Phi matrices. LAPACK wants left layout
293 Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type> // left layout for lapack
294 PsiMat("PsiMat", ndofSubcell, ndofSubcell),
295 PhiMat("PhiMat", ndofSubcell, ndofSubcell);
296
297 auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
298 auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
299
300 for (ordinal_type i=0;i<ndofSubcell;++i) {
301 const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
302 for (ordinal_type j=0;j<ndofSubcell;++j) {
303 const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
304 value_type refEntry = 0, ortEntry =0;
305 for (ordinal_type k=0;k<subcellDim;++k) {
306 ortEntry += subCellValues(isc,j,k)*refSubcellTangents(j,k);
307 for (ordinal_type d=0; d<cellDim; ++d)
308 refEntry += cellBasisValues(ic,j,d)*trJacobianF(k,d)*refSubcellTangents(j,k);
309 }
310 PsiMat(j,i) = refEntry;
311 PhiMat(j,i) = ortEntry;
312 }
313 }
314
315 // Solve the system using Lapack
316 {
317 Teuchos::LAPACK<ordinal_type,value_type> lapack;
318 ordinal_type info = 0;
319
320
321 /*
322 Kokkos::View<value_type**,Kokkos::LayoutLeft,host_space_type> work("pivVec", 2*ndofSubcell, 1);
323 lapack.GELS('N', ndofSubcell*card, ndofSubcell, ndofSubcell,
324 PsiMat.data(),
325 PsiMat.stride_1(),
326 PhiMat.data(),
327 PhiMat.stride_1(),
328 work.data(), work.extent(0),
329 &info);
330
331 */
332 Kokkos::DynRankView<ordinal_type,host_device_type> pivVec("pivVec", ndofSubcell);
333 lapack.GESV(ndofSubcell, ndofSubcell,
334 PsiMat.data(),
335 PhiMat.stride_1(),
336 pivVec.data(),
337 PhiMat.data(),
338 PhiMat.stride_1(),
339 &info);
340 //*/
341 if (info) {
342 std::stringstream ss;
343 ss << ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): "
344 << "LAPACK return with error code: "
345 << info;
346 INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
347 }
348
349 //After solving the system w/ LAPACK, Phi contains A^T
350
351 // transpose and clean up numerical noise (for permutation matrices)
352 const double eps = tolerence();
353 for (ordinal_type i=0;i<ndofSubcell;++i) {
354 auto intmatii = std::round(PhiMat(i,i));
355 PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
356 for (ordinal_type j=i+1;j<ndofSubcell;++j) {
357 auto matij = PhiMat(i,j);
358
359 auto intmatji = std::round(PhiMat(j,i));
360 PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
361
362 auto intmatij = std::round(matij);
363 PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
364 }
365 }
366
367
368
369 // Print A matrix
370 /*
371 {
372 std::cout << "|";
373 for (ordinal_type i=0;i<ndofSubcell;++i) {
374 for (ordinal_type j=0;j<ndofSubcell;++j) {
375 std::cout << PhiMat(i,j) << " ";
376 }
377 std::cout << "| ";
378 }
379 std::cout <<std::endl;
380 }
381 */
382 }
383
384 {
385 // move the data to original device memory
386 const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
387 auto suboutput = Kokkos::subview(output, range, range);
388 auto tmp = Kokkos::create_mirror_view_and_copy(typename OutputViewType::device_type::memory_space(), PhiMat);
389 Kokkos::deep_copy(suboutput, tmp);
390 }
391}
392}
393
394}
395#endif
Header file for the Intrepid2::Basis_HCURL_QUAD_In_FEM class.
Header file for the Intrepid2::Basis_HCURL_TRI_In_FEM class.
Header file for the Intrepid2::Basis_HVOL_LINE_Cn_FEM class.
Implementation of the default H(curl)-compatible FEM basis on Quadrilateral cell.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Tr...
Implementation of the locally HVOL-compatible FEM basis of variable order on the [-1,...
static void getJacobianOfOrientationMap(JacobianViewType jacobian, const shards::CellTopology cellTopo, const ordinal_type cellOrt)
Computes jacobian of the parameterization maps of 1- and 2-subcells with orientation.
static KOKKOS_INLINE_FUNCTION void getRefSubcellTangents(TanViewType tangents, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) subCell tangents.
static void mapToModifiedReference(outPointViewType outPoints, const refPointViewType refPoints, const shards::CellTopology cellTopo, const ordinal_type cellOrt=0)
Computes modified parameterization maps of 1- and 2-subcells with orientation.
static void getCoeffMatrix_HCURL(OutputViewType &output, const subcellBasisHostType &subcellBasis, const cellBasisHostType &cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute orientation matrix for HCURL basis for a given subcell and its reference basis.
static KOKKOS_INLINE_FUNCTION void mapSubcellCoordsToRefCell(coordsViewType cellCoords, const subcellCoordsViewType subCellCoords, const ParamViewType subcellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Maps points defined on the subCell manifold into the parent Cell.
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...