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 
74 namespace Intrepid2 {
75 
76 namespace Impl {
77 
78 namespace Debug {
79 
80 #ifdef HAVE_INTREPID2_DEBUG
81 template<typename subcellBasisType,
82 typename cellBasisType>
83 inline
84 void
85 check_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 must be greater 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 
169 template<typename OutputViewType,
170 typename subcellBasisHostType,
171 typename cellBasisHostType>
172 inline
173 void
175 getCoeffMatrix_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  auto subcellParam = Intrepid2::RefSubcellParametrization<host_device_type>::get(subcellDim, cellTopo.getKey());
259 
260  // refPtsCell = F_s (\eta_o (refPtsSubcell))
261  Kokkos::DynRankView<value_type,host_device_type> refPtsCell("refPtsCell", ndofSubcell, cellDim);
262  mapSubcellCoordsToRefCell(refPtsCell,refPtsSubcell, subcellParam, subcellBaseKey, subcellId, subcellOrt);
263 
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> trJacobianF("trJacobianF", subcellDim, cellDim );
267  OrientationTools::getRefSubcellTangents(trJacobianF, subcellParam, subcellBaseKey, subcellId, subcellOrt);
268 
269 
270 
271  // cellBasisValues = \psi_k(F_s (\eta_o (\xi_j)))
272  Kokkos::DynRankView<value_type,host_device_type> cellBasisValues("cellBasisValues", numCellBasis, ndofSubcell, cellDim);
273  cellBasis.getValues(cellBasisValues, refPtsCell, OPERATOR_VALUE);
274 
275  //
276  // Compute Psi_jk = \psi_k(F_s (\eta_o (\xi_j))) \cdot (J_F J_\eta t_j)
277  // and Phi_ji = \phi_i (\xi_j) \ctot t_j, and solve
278  // Psi A^T = Phi
279  //
280 
281  // construct Psi and Phi matrices. LAPACK wants left layout
282  Kokkos::DynRankView<value_type,Kokkos::LayoutLeft,host_device_type> // left layout for lapack
283  PsiMat("PsiMat", ndofSubcell, ndofSubcell),
284  PhiMat("PhiMat", ndofSubcell, ndofSubcell);
285 
286  auto cellTagToOrdinal = cellBasis.getAllDofOrdinal();
287  auto subcellTagToOrdinal = subcellBasis.getAllDofOrdinal();
288 
289  for (ordinal_type i=0;i<ndofSubcell;++i) {
290  const ordinal_type ic = cellTagToOrdinal(subcellDim, subcellId, i);
291  for (ordinal_type j=0;j<ndofSubcell;++j) {
292  const ordinal_type isc = subcellTagToOrdinal(subcellDim, 0, i);
293  value_type refEntry = 0, ortEntry =0;
294  for (ordinal_type k=0;k<subcellDim;++k) {
295  ortEntry += subCellValues(isc,j,k)*refSubcellTangents(j,k);
296  for (ordinal_type d=0; d<cellDim; ++d)
297  refEntry += cellBasisValues(ic,j,d)*trJacobianF(k,d)*refSubcellTangents(j,k);
298  }
299  PsiMat(j,i) = refEntry;
300  PhiMat(j,i) = ortEntry;
301  }
302  }
303 
304  // Solve the system using Lapack
305  {
306  Teuchos::LAPACK<ordinal_type,value_type> lapack;
307  ordinal_type info = 0;
308 
309 
310  /*
311  Kokkos::View<value_type**,Kokkos::LayoutLeft,host_space_type> work("pivVec", 2*ndofSubcell, 1);
312  lapack.GELS('N', ndofSubcell*card, ndofSubcell, ndofSubcell,
313  PsiMat.data(),
314  PsiMat.stride_1(),
315  PhiMat.data(),
316  PhiMat.stride_1(),
317  work.data(), work.extent(0),
318  &info);
319 
320  */
321  Kokkos::DynRankView<ordinal_type,host_device_type> pivVec("pivVec", ndofSubcell);
322  lapack.GESV(ndofSubcell, ndofSubcell,
323  PsiMat.data(),
324  PhiMat.stride_1(),
325  pivVec.data(),
326  PhiMat.data(),
327  PhiMat.stride_1(),
328  &info);
329  //*/
330  if (info) {
331  std::stringstream ss;
332  ss << ">>> ERROR (Intrepid::OrientationTools::getCoeffMatrix_HCURL): "
333  << "LAPACK return with error code: "
334  << info;
335  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
336  }
337 
338  //After solving the system w/ LAPACK, Phi contains A^T
339 
340  // transpose and clean up numerical noise (for permutation matrices)
341  const double eps = tolerence();
342  for (ordinal_type i=0;i<ndofSubcell;++i) {
343  auto intmatii = std::round(PhiMat(i,i));
344  PhiMat(i,i) = (std::abs(PhiMat(i,i) - intmatii) < eps) ? intmatii : PhiMat(i,i);
345  for (ordinal_type j=i+1;j<ndofSubcell;++j) {
346  auto matij = PhiMat(i,j);
347 
348  auto intmatji = std::round(PhiMat(j,i));
349  PhiMat(i,j) = (std::abs(PhiMat(j,i) - intmatji) < eps) ? intmatji : PhiMat(j,i);
350 
351  auto intmatij = std::round(matij);
352  PhiMat(j,i) = (std::abs(matij - intmatij) < eps) ? intmatij : matij;
353  }
354  }
355 
356 
357 
358  // Print A matrix
359  /*
360  {
361  std::cout << "|";
362  for (ordinal_type i=0;i<ndofSubcell;++i) {
363  for (ordinal_type j=0;j<ndofSubcell;++j) {
364  std::cout << PhiMat(i,j) << " ";
365  }
366  std::cout << "| ";
367  }
368  std::cout <<std::endl;
369  }
370  */
371  }
372 
373  {
374  // move the data to original device memory
375  const Kokkos::pair<ordinal_type,ordinal_type> range(0, ndofSubcell);
376  auto suboutput = Kokkos::subview(output, range, range);
377  auto tmp = Kokkos::create_mirror_view_and_copy(typename OutputViewType::device_type::memory_space(), PhiMat);
378  Kokkos::deep_copy(suboutput, tmp);
379  }
380 }
381 }
382 
383 }
384 #endif
Implementation of the default H(curl)-compatible FEM basis on Quadrilateral cell. ...
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...
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 getRefSubcellTangents(TanViewType tangents, const ParamViewType subCellParametrization, const unsigned subcellTopoKey, const ordinal_type subCellOrd, const ordinal_type ort)
Computes the (oriented) subCell tangents.
Header file for the Intrepid2::Basis_HVOL_LINE_Cn_FEM class.
Header file for the Intrepid2::Basis_HCURL_QUAD_In_FEM class.
Implementation of the locally HVOL-compatible FEM basis of variable order on the [-1,1] reference line cell, using Lagrange polynomials.
Implementation of the default H(curl)-compatible Nedelec (first kind) basis of arbitrary degree on Tr...
Header file for the Intrepid2::Basis_HCURL_TRI_In_FEM class.