Intrepid2
Intrepid2_HCURL_TET_In_FEMDef.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Intrepid2 Package
4 //
5 // Copyright 2007 NTESS and the Intrepid2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
16 #ifndef __INTREPID2_HCURL_TET_IN_FEM_DEF_HPP__
17 #define __INTREPID2_HCURL_TET_IN_FEM_DEF_HPP__
18 
19 #include "Intrepid2_HGRAD_TET_Cn_FEM_ORTH.hpp"
20 #include "Intrepid2_CubatureDirectTetDefault.hpp"
21 #include "Teuchos_SerialDenseMatrix.hpp"
22 
23 namespace Intrepid2 {
24 
25 // -------------------------------------------------------------------------------------
26 
27 namespace Impl {
28 
29 template<EOperator OpType>
30 template<typename OutputViewType,
31 typename InputViewType,
32 typename WorkViewType,
33 typename VinvViewType>
34 KOKKOS_INLINE_FUNCTION
35 void
36 Basis_HCURL_TET_In_FEM::Serial<OpType>::
37 getValues( OutputViewType output,
38  const InputViewType input,
39  WorkViewType work,
40  const VinvViewType coeffs ) {
41 
42  constexpr ordinal_type spaceDim = 3;
43  const ordinal_type
44  cardPn = coeffs.extent(0)/spaceDim,
45  card = coeffs.extent(1),
46  npts = input.extent(0);
47 
48  // compute order
49  ordinal_type order = 0;
50  for (ordinal_type p=0;p<=Parameters::MaxOrder;++p) {
51  if (card == CardinalityHCurlTet(p)) {
52  order = p;
53  break;
54  }
55  }
56 
57  typedef typename Kokkos::DynRankView<typename InputViewType::value_type, typename WorkViewType::memory_space> ViewType;
58  auto vcprop = Kokkos::common_view_alloc_prop(input);
59  auto ptr = work.data();
60 
61  switch (OpType) {
62  case OPERATOR_VALUE: {
63  const ViewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts);
64  ViewType dummyView;
65 
66  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
67  Serial<OpType>::getValues(phis, input, dummyView, order);
68 
69  for (ordinal_type i=0;i<card;++i)
70  for (ordinal_type j=0;j<npts;++j)
71  for (ordinal_type d=0;d<spaceDim;++d) {
72  output.access(i,j,d) = 0.0;
73  for (ordinal_type k=0;k<cardPn;++k)
74  output.access(i,j,d) += coeffs(k+d*cardPn,i) * phis(k,j);
75  }
76  break;
77  }
78  case OPERATOR_CURL: {
79  const ViewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts, spaceDim);
80  ptr += card*npts*spaceDim*get_dimension_scalar(input);
81  const ViewType workView(Kokkos::view_wrap(ptr, vcprop), card, npts, spaceDim+1);
82 
83  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
84  Serial<OPERATOR_GRAD>::getValues(phis, input, workView, order);
85 
86  for (ordinal_type i=0;i<card;++i) {
87  for (ordinal_type j=0;j<npts;++j) {
88  for (ordinal_type d=0; d< spaceDim; ++d) {
89  output.access(i,j,d) = 0.0;
90  ordinal_type d1 = (d+1) % spaceDim, d2 = (d+2) % spaceDim;
91  for (ordinal_type k=0; k<cardPn; ++k) //\sum_k (coeffs_k, coeffs_{k+cardPn}, coeffs_{k+2 cardPn}) \times phis_kj (cross product)
92  output.access(i,j,d) += coeffs(k+d2*cardPn,i)*phis(k,j,d1)
93  -coeffs(k+d1*cardPn,i)*phis(k,j,d2);
94  }
95  }
96  }
97  break;
98  }
99  default: {
100  INTREPID2_TEST_FOR_ABORT( true,
101  ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");
102  }
103  }
104 }
105 
106 template<typename DT, ordinal_type numPtsPerEval,
107 typename outputValueValueType, class ...outputValueProperties,
108 typename inputPointValueType, class ...inputPointProperties,
109 typename vinvValueType, class ...vinvProperties>
110 void
111 Basis_HCURL_TET_In_FEM::
112 getValues(
113  const typename DT::execution_space& space,
114  Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
115  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
116  const Kokkos::DynRankView<vinvValueType, vinvProperties...> coeffs,
117  const EOperator operatorType) {
118  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
119  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
120  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
121  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
122 
123  // loopSize corresponds to cardinality
124  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
125  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
126  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
127  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(space, 0, loopSize);
128 
129  typedef typename inputPointViewType::value_type inputPointType;
130 
131  const ordinal_type cardinality = outputValues.extent(0);
132  const ordinal_type spaceDim = 3;
133 
134  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
135  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
136 
137  switch (operatorType) {
138  case OPERATOR_VALUE: {
139  workViewType work(Kokkos::view_alloc(space, "Basis_HCURL_TET_In_FEM::getValues::work", vcprop), cardinality, inputPoints.extent(0));
140  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
141  OPERATOR_VALUE,numPtsPerEval> FunctorType;
142  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, coeffs, work) );
143  break;
144  }
145  case OPERATOR_CURL: {
146  workViewType work(Kokkos::view_alloc(space, "Basis_HCURL_TET_In_FEM::getValues::work", vcprop), cardinality*(2*spaceDim+1), inputPoints.extent(0));
147  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
148  OPERATOR_CURL,numPtsPerEval> FunctorType;
149  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, coeffs, work) );
150  break;
151  }
152  default: {
153  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
154  ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented" );
155  }
156  }
157 }
158 }
159 
160 // -------------------------------------------------------------------------------------
161 template<typename DT, typename OT, typename PT>
162 Basis_HCURL_TET_In_FEM<DT,OT,PT>::
163 Basis_HCURL_TET_In_FEM( const ordinal_type order,
164  const EPointType pointType ) {
165 
166  constexpr ordinal_type spaceDim = 3;
167  this->basisCardinality_ = CardinalityHCurlTet(order);
168  this->basisDegree_ = order; // small n
169  this->basisCellTopologyKey_ = shards::Tetrahedron<4>::key;
170  this->basisType_ = BASIS_FEM_LAGRANGIAN;
171  this->basisCoordinates_ = COORDINATES_CARTESIAN;
172  this->functionSpace_ = FUNCTION_SPACE_HCURL;
173  pointType_ = pointType;
174  const ordinal_type card = this->basisCardinality_;
175 
176  const ordinal_type cardPn = Intrepid2::getPnCardinality<spaceDim>(order); // dim of (P_{n}) -- smaller space
177  const ordinal_type cardPnm1 = Intrepid2::getPnCardinality<spaceDim>(order-1); // dim of (P_{n-1}) -- smaller space
178  const ordinal_type cardPnm2 = Intrepid2::getPnCardinality<spaceDim>(order-2); // dim of (P_{n-2}) -- smaller space
179  const ordinal_type cardVecPn = spaceDim*cardPn; // dim of (P_{n})^2 -- larger space
180  const ordinal_type cardVecPnm1 = spaceDim*cardPnm1; // dim of (P_{n-1})^2 -- smaller space
181  const ordinal_type cardPnm1H = cardPnm1-cardPnm2; //Homogeneous polynomial of order (n-1)
182 
183  // Note: the only reason why equispaced can't support higher order than Parameters::MaxOrder appears to be the fact that the tags below get stored into a fixed-length array.
184  // TODO: relax the maximum order requirement by setting up tags in a different container, perhaps directly into an OrdinalTypeArray1DHost (tagView, below). (As of this writing (1/25/22), looks like other nodal bases do this in a similar way -- those should be fixed at the same time; maybe search for Parameters::MaxOrder.)
185  INTREPID2_TEST_FOR_EXCEPTION( order > Parameters::MaxOrder, std::invalid_argument, "polynomial order exceeds the max supported by this class");
186 
187  // Basis-dependent initializations
188  constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
189  constexpr ordinal_type maxCard = CardinalityHCurlTet(Parameters::MaxOrder);
190  ordinal_type tags[maxCard][tagSize];
191 
192  // points are computed in the host and will be copied
193  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
194  dofCoords("Hcurl::Tet::In::dofCoords", card, spaceDim);
195 
196  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
197  coeffs("Hcurl::Tet::In::coeffs", cardVecPn, card);
198 
199  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
200  dofCoeffs("Hcurl::Tet::In::dofCoeffs", card, spaceDim);
201 
202  // first, need to project the basis for RT space onto the
203  // orthogonal basis of degree n
204  // get coefficients of PkHx
205 
206  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace> //use LayoutLeft for Lapack
207  V1("Hcurl::Tet::In::V1", cardVecPn, cardVecPnm1 + spaceDim*cardPnm1H);
208 
209 
210  // these two loops get the first three sets of basis functions
211  for (ordinal_type i=0;i<cardPnm1;i++)
212  for (ordinal_type d=0;d<spaceDim;d++)
213  V1(i+d*cardPn,i+d*cardPnm1) = 1.0;
214 
215 
216  // now I need to integrate { (x,y) \times phi } against the big basis
217  // first, get a cubature rule.
218  CubatureDirectTetDefault<Kokkos::HostSpace::execution_space,scalarType,scalarType> myCub( 2 * order );
219  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubPoints("Hcurl::Tet::In::cubPoints", myCub.getNumPoints() , spaceDim );
220  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubWeights("Hcurl::Tet::In::cubWeights", myCub.getNumPoints() );
221  myCub.getCubature( cubPoints , cubWeights );
222 
223  // tabulate the scalar orthonormal basis at cubature points
224  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> phisAtCubPoints("Hcurl::Tet::In::phisAtCubPoints", cardPn , myCub.getNumPoints() );
225  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(typename Kokkos::HostSpace::execution_space{},
226  phisAtCubPoints,
227  cubPoints,
228  order,
229  OPERATOR_VALUE);
230 
231  // Integrate (x psi_j, y psi_j, z psi_j) \times (phi_i, phi_{i+cardPn}, phi_{i+2 cardPn}) cross product. psi are homogeneous polynomials of order (n-1)
232  for (ordinal_type i=0;i<cardPn;i++) {
233  for (ordinal_type j=0;j<cardPnm1H;j++) { // loop over homogeneous polynomials
234  for (ordinal_type d=0; d< spaceDim; ++d) {
235  scalarType integral(0);
236  for (ordinal_type k=0;k<myCub.getNumPoints();k++)
237  integral += cubWeights(k) * cubPoints(k,d)
238  * phisAtCubPoints(cardPnm2+j,k)
239  * phisAtCubPoints(i,k);
240  ordinal_type d1 = (d+1) % spaceDim, d2 = (d+2) % spaceDim;
241  V1(i+d2*cardPn,cardVecPnm1+d1*cardPnm1H + j) = -integral;
242  V1(i+d1*cardPn,cardVecPnm1+d2*cardPnm1H + j) = integral;
243  }
244  }
245  }
246 
247 
248 
249 
250 
251  // now I need to set up an SVD to get a basis for the space
252  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
253  S("Hcurl::Tet::In::S", cardVecPn,1),
254  U("Hcurl::Tet::In::U", cardVecPn, cardVecPn),
255  Vt("Hcurl::Tet::In::Vt", cardVecPn, cardVecPn),
256  work("Hcurl::Tet::In::work", 5*cardVecPn,1),
257  rWork("Hcurl::Tet::In::rW", 1,1);
258 
259 
260 
261  ordinal_type info = 0;
262  Teuchos::LAPACK<ordinal_type,scalarType> lapack;
263 
264 
265  lapack.GESVD( 'A',
266  'N',
267  V1.extent(0) ,
268  V1.extent(1) ,
269  V1.data() ,
270  V1.stride_1() ,
271  S.data() ,
272  U.data() ,
273  U.stride_1() ,
274  Vt.data() ,
275  Vt.stride_1() ,
276  work.data() ,
277  5*cardVecPn ,
278  rWork.data() ,
279  &info );
280 
281 
282 #ifdef HAVE_INTREPID2_DEBUG
283  ordinal_type num_nonzero_sv = 0;
284  for (int i=0;i<cardVecPn;i++)
285  num_nonzero_sv += (S(i,0) > 10*tolerence());
286 
287  INTREPID2_TEST_FOR_EXCEPTION( num_nonzero_sv != card, std::invalid_argument,
288  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM( order, pointType), Matrix V1 should have rank equal to the cardinality of HCURL space");
289 #endif
290 
291  // next, apply the RT nodes (rows) to the basis for (P_n)^2 (columns)
292  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
293  V2("Hcurl::Tet::In::V2", card ,cardVecPn);
294 
295  shards::CellTopology cellTopo(shards::getCellTopologyData<shards::Tetrahedron<4> >());
296  const ordinal_type numEdges = cellTopo.getEdgeCount();
297  const ordinal_type numFaces = cellTopo.getFaceCount();
298 
299  // first numEdges * degree nodes are normals at each edge
300  // get the points on the line
301 
302  shards::CellTopology edgeTopo(shards::getCellTopologyData<shards::Line<2> >() );
303  shards::CellTopology faceTopo(shards::getCellTopologyData<shards::Triangle<3> >() );
304 
305  const int numPtsPerEdge = PointTools::getLatticeSize( edgeTopo ,
306  order+1 ,
307  1 );
308 
309  const int numPtsPerFace = PointTools::getLatticeSize( faceTopo ,
310  order+1 ,
311  1 );
312 
313  const int numPtsPerCell = PointTools::getLatticeSize( cellTopo ,
314  order+1 ,
315  1 );
316 
317  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> linePts("Hcurl::Tet::In::linePts", numPtsPerEdge , 1 );
318  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> triPts("Hcurl::Tet::In::triPts", numPtsPerFace , 2 );
319 
320  // construct lattice
321  const ordinal_type offset = 1;
322 
323 
324 
325  PointTools::getLattice( linePts,
326  edgeTopo,
327  order+1, offset,
328  pointType );
329 
330  PointTools::getLattice( triPts,
331  faceTopo,
332  order+1, offset,
333  pointType );
334 
335  // holds the image of the line points
336  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> edgePts("Hcurl::Tet::In::edgePts", numPtsPerEdge , spaceDim );
337  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> facePts("Hcurl::Tet::In::facePts", numPtsPerFace , spaceDim );
338  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> phisAtEdgePoints("Hcurl::Tet::In::phisAtEdgePoints", cardPn , numPtsPerEdge );
339  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> phisAtFacePoints("Hcurl::Tet::In::phisAtFacePoints", cardPn , numPtsPerFace);
340 
341  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> edgeTan("Hcurl::Tet::In::edgeTan", spaceDim );
342 
343  // these are tangents scaled by the appropriate edge lengths.
344  for (ordinal_type i=0;i<numEdges;i++) { // loop over edges
345  CellTools<Kokkos::HostSpace>::getReferenceEdgeTangent( edgeTan ,
346  i ,
347  cellTopo );
348 
349  CellTools<Kokkos::HostSpace>::mapToReferenceSubcell( edgePts ,
350  linePts ,
351  1 ,
352  i ,
353  cellTopo );
354 
355  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(typename Kokkos::HostSpace::execution_space{},
356  phisAtEdgePoints,
357  edgePts,
358  order,
359  OPERATOR_VALUE);
360 
361  // loop over points (rows of V2)
362  for (ordinal_type j=0;j<numPtsPerEdge;j++) {
363 
364  const ordinal_type i_card = numPtsPerEdge*i+j;
365 
366  // loop over orthonormal basis functions (columns of V2)
367  for (ordinal_type k=0;k<cardPn;k++)
368  for (ordinal_type d=0;d<spaceDim;d++)
369  V2(i_card,k+d*cardPn) = edgeTan(d) * phisAtEdgePoints(k,j);
370 
371  //save dof coordinates and coefficients
372  for(ordinal_type k=0; k<spaceDim; ++k) {
373  dofCoords(i_card,k) = edgePts(j,k);
374  dofCoeffs(i_card,k) = edgeTan(k);
375  }
376 
377  tags[i_card][0] = 1; // edge dof
378  tags[i_card][1] = i; // edge id
379  tags[i_card][2] = j; // local dof id
380  tags[i_card][3] = numPtsPerEdge; // total vert dof
381 
382  }
383  }
384 
385  if(numPtsPerFace >0) {//handle faces if needed (order >1)
386  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> faceTan1("Hcurl::Tet::In::edgeTan", spaceDim );
387  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> faceTan2("Hcurl::Tet::In::edgeTan", spaceDim );
388 
389  for (ordinal_type i=0;i<numFaces;i++) { // loop over faces
390  CellTools<Kokkos::HostSpace>::getReferenceFaceTangents( faceTan1 ,
391  faceTan2,
392  i ,
393  cellTopo );
394 
395  CellTools<Kokkos::HostSpace>::mapToReferenceSubcell( facePts ,
396  triPts ,
397  2 ,
398  i ,
399  cellTopo );
400 
401  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(typename Kokkos::HostSpace::execution_space{},
402  phisAtFacePoints,
403  facePts,
404  order,
405  OPERATOR_VALUE);
406 
407  // loop over points (rows of V2)
408  for (ordinal_type j=0;j<numPtsPerFace;j++) {
409 
410  const ordinal_type i_card = numEdges*numPtsPerEdge+2*numPtsPerFace*i+2*j;
411  const ordinal_type i_card_p1 = i_card+1; // creating a temp otherwise nvcc gets confused
412 
413  // loop over orthonormal basis functions (columns of V2)
414  for (ordinal_type k=0;k<cardPn;k++)
415  for (ordinal_type d=0;d<spaceDim;d++) {
416  V2(i_card,k+d*cardPn) = faceTan1(d) * phisAtFacePoints(k,j);
417  V2(i_card_p1,k+d*cardPn) = faceTan2(d) * phisAtFacePoints(k,j);
418  }
419 
420  //save dof coordinates
421  for(ordinal_type k=0; k<spaceDim; ++k) {
422  dofCoords(i_card,k) = facePts(j,k);
423  dofCoords(i_card_p1,k) = facePts(j,k);
424  dofCoeffs(i_card,k) = faceTan1(k);
425  dofCoeffs(i_card_p1,k) = faceTan2(k);
426  }
427 
428  tags[i_card][0] = 2; // face dof
429  tags[i_card][1] = i; // face id
430  tags[i_card][2] = 2*j; // local face id
431  tags[i_card][3] = 2*numPtsPerFace; // total face dof
432 
433  tags[i_card_p1][0] = 2; // face dof
434  tags[i_card_p1][1] = i; // face id
435  tags[i_card_p1][2] = 2*j+1; // local face id
436  tags[i_card_p1][3] = 2*numPtsPerFace; // total face dof
437 
438  }
439  }
440  }
441 
442 
443  // internal dof, if needed
444  if (numPtsPerCell > 0) {
445  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
446  cellPoints( "Hcurl::Tet::In::cellPoints", numPtsPerCell , spaceDim );
447  PointTools::getLattice( cellPoints ,
448  cellTopo ,
449  order + 1 ,
450  1 ,
451  pointType );
452 
453  Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
454  phisAtCellPoints("Hcurl::Tet::In::phisAtCellPoints", cardPn , numPtsPerCell );
455  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(typename Kokkos::HostSpace::execution_space{},
456  phisAtCellPoints,
457  cellPoints,
458  order,
459  OPERATOR_VALUE);
460 
461  // copy values into right positions of V2
462  for (ordinal_type j=0;j<numPtsPerCell;j++) {
463 
464  const ordinal_type i_card = numEdges*numPtsPerEdge+2*numFaces*numPtsPerFace+spaceDim*j;
465 
466  for (ordinal_type k=0;k<cardPn;k++)
467  for (ordinal_type d=0;d<spaceDim;d++)
468  V2(i_card+d,d*cardPn+k) = phisAtCellPoints(k,j);
469 
470 
471  //save dof coordinates
472  for(ordinal_type d=0; d<spaceDim; ++d) {
473  for(ordinal_type dim=0; dim<spaceDim; ++dim) {
474  dofCoords(i_card+d,dim) = cellPoints(j,dim);
475  dofCoeffs(i_card+d,dim) = (d==dim);
476  }
477 
478  tags[i_card+d][0] = spaceDim; // elem dof
479  tags[i_card+d][1] = 0; // elem id
480  tags[i_card+d][2] = spaceDim*j+d; // local dof id
481  tags[i_card+d][3] = spaceDim*numPtsPerCell; // total vert dof
482  }
483  }
484  }
485 
486  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
487  // so we transpose on copy below.
488  const ordinal_type lwork = card*card;
489  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
490  vmat("Hcurl::Tet::In::vmat", card, card),
491  work1("Hcurl::Tet::In::work", lwork),
492  ipiv("Hcurl::Tet::In::ipiv", card);
493 
494  //vmat = V2*U;
495  for(ordinal_type i=0; i< card; ++i) {
496  for(ordinal_type j=0; j< card; ++j) {
497  scalarType s=0;
498  for(ordinal_type k=0; k< cardVecPn; ++k)
499  s += V2(i,k)*U(k,j);
500  vmat(i,j) = s;
501  }
502  }
503 
504  info = 0;
505 
506  lapack.GETRF(card, card,
507  vmat.data(), vmat.stride_1(),
508  (ordinal_type*)ipiv.data(),
509  &info);
510 
511  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
512  std::runtime_error ,
513  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRF returns nonzero info." );
514 
515  lapack.GETRI(card,
516  vmat.data(), vmat.stride_1(),
517  (ordinal_type*)ipiv.data(),
518  work1.data(), lwork,
519  &info);
520 
521  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
522  std::runtime_error ,
523  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRI returns nonzero info." );
524 
525  for (ordinal_type i=0;i<cardVecPn;++i) {
526  for (ordinal_type j=0;j<card;++j){
527  scalarType s=0;
528  for(ordinal_type k=0; k< card; ++k)
529  s += U(i,k)*vmat(k,j);
530  coeffs(i,j) = s;
531  }
532  }
533 
534  this->coeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), coeffs);
535  Kokkos::deep_copy(this->coeffs_ , coeffs);
536 
537  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
538  Kokkos::deep_copy(this->dofCoords_, dofCoords);
539 
540  this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffs);
541  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
542 
543 
544  // set tags
545  {
546  // Basis-dependent initializations
547  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
548  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
549  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
550 
551  OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
552 
553  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
554  // tags are constructed on host
555  this->setOrdinalTagData(this->tagToOrdinal_,
556  this->ordinalToTag_,
557  tagView,
558  this->basisCardinality_,
559  tagSize,
560  posScDim,
561  posScOrd,
562  posDfOrd);
563  }
564 }
565 
566 template<typename DT, typename OT, typename PT>
567 void
568 Basis_HCURL_TET_In_FEM<DT,OT,PT>::getScratchSpaceSize(
569  ordinal_type& perTeamSpaceSize,
570  ordinal_type& perThreadSpaceSize,
571  const PointViewType inputPoints,
572  const EOperator operatorType) const {
573  perTeamSpaceSize = 0;
574  ordinal_type scalarWorkViewExtent = (operatorType == OPERATOR_VALUE) ? this->basisCardinality_ : 7*this->basisCardinality_;
575  perThreadSpaceSize = scalarWorkViewExtent*get_dimension_scalar(inputPoints)*sizeof(typename BasisBase::scalarType);
576 }
577 
578 template<typename DT, typename OT, typename PT>
579 KOKKOS_INLINE_FUNCTION
580 void
581 Basis_HCURL_TET_In_FEM<DT,OT,PT>::getValues(
582  OutputViewType outputValues,
583  const PointViewType inputPoints,
584  const EOperator operatorType,
585  const typename Kokkos::TeamPolicy<typename DT::execution_space>::member_type& team_member,
586  const typename DT::execution_space::scratch_memory_space & scratchStorage,
587  const ordinal_type subcellDim,
588  const ordinal_type subcellOrdinal) const {
589 
590  INTREPID2_TEST_FOR_ABORT( !((subcellDim == -1) && (subcellOrdinal == -1)),
591  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM::getValues), The capability of selecting subsets of basis functions has not been implemented yet.");
592 
593  const int numPoints = inputPoints.extent(0);
594  using ScalarType = typename ScalarTraits<typename PointViewType::value_type>::scalar_type;
595  using WorkViewType = Kokkos::DynRankView< ScalarType,typename DT::execution_space::scratch_memory_space,Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
596  ordinal_type scalarSizePerPoint = (operatorType == OPERATOR_VALUE) ? this->basisCardinality_ : 7*this->basisCardinality_;
597  ordinal_type sizePerPoint = scalarSizePerPoint*get_dimension_scalar(inputPoints);
598  WorkViewType workView(scratchStorage, sizePerPoint*team_member.team_size());
599  using range_type = Kokkos::pair<ordinal_type,ordinal_type>;
600 
601  switch(operatorType) {
602  case OPERATOR_VALUE:
603  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
604  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type (pt,pt+1), Kokkos::ALL() );
605  const auto input = Kokkos::subview( inputPoints, range_type(pt, pt+1), Kokkos::ALL() );
606  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
607  Impl::Basis_HCURL_TET_In_FEM::Serial<OPERATOR_VALUE>::getValues( output, input, work, this->coeffs_ );
608  });
609  break;
610  case OPERATOR_CURL:
611  Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=] (ordinal_type& pt) {
612  auto output = Kokkos::subview( outputValues, Kokkos::ALL(), range_type(pt,pt+1), Kokkos::ALL() );
613  const auto input = Kokkos::subview( inputPoints, range_type(pt,pt+1), Kokkos::ALL() );
614  WorkViewType work(workView.data() + sizePerPoint*team_member.team_rank(), sizePerPoint);
615  Impl::Basis_HCURL_TET_In_FEM::Serial<OPERATOR_CURL>::getValues( output, input, work, this->coeffs_ );
616  });
617  break;
618  default: {
619  INTREPID2_TEST_FOR_ABORT( true,
620  ">>> ERROR (Basis_HCURL_TET_In_FEM): getValues not implemented for this operator");
621  }
622  }
623 }
624 } // namespace Intrepid2
625 #endif
Intrepid2::Basis< DeviceType, outputValueType, pointValueType >::scalarType
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
Definition: Intrepid2_Basis.hpp:157
Intrepid2::CubatureDirect::getCubature
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights) const override
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Definition: Intrepid2_CubatureDirect.hpp:155
Intrepid2::Basis_HCURL_TET_In_FEM::Basis_HCURL_TET_In_FEM
Basis_HCURL_TET_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition: Intrepid2_HCURL_TET_In_FEMDef.hpp:163
Intrepid2::PointTools::getLattice
static void getLattice(Kokkos::DynRankView< pointValueType, pointProperties...> points, const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference simplex, quadrilateral or hexahedron (cu...
Definition: Intrepid2_PointToolsDef.hpp:86
Intrepid2::Basis_HCURL_TET_In_FEM::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition: Intrepid2_Basis.hpp:338
Intrepid2::CubatureDirect::getNumPoints
virtual ordinal_type getNumPoints() const override
Returns the number of cubature points.
Definition: Intrepid2_CubatureDirect.hpp:164
Intrepid2::CellTools::getReferenceEdgeTangent
static void getReferenceEdgeTangent(RefEdgeTangentViewType refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:211
Intrepid2_CubatureDirectTetDefault.hpp
Header file for the Intrepid2::CubatureDirectTetDefault class.
Intrepid2::Basis_HCURL_TET_In_FEM::getScratchSpaceSize
virtual void getScratchSpaceSize(ordinal_type &perTeamSpaceSize, ordinal_type &perThreadSpaceSize, const PointViewType inputPointsconst, const EOperator operatorType=OPERATOR_VALUE) const override
Return the size of the scratch space, in bytes, needed for using the team-level implementation of get...
Definition: Intrepid2_HCURL_TET_In_FEMDef.hpp:568
Intrepid2::CellTools::mapToReferenceSubcell
static void mapToReferenceSubcell(refSubcellViewType refSubcellPoints, const paramPointViewType paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Intrepid2::Basis_HCURL_TET_In_FEM::getValues
virtual void getValues(const ExecutionSpace &space, OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
Definition: Intrepid2_HCURL_TET_In_FEM.hpp:199
Intrepid2::CubatureDirectTetDefault
Defines direct integration rules on a tetrahedron.
Definition: Intrepid2_CubatureDirectTetDefault.hpp:30
Intrepid2::Basis_HCURL_TET_In_FEM::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition: Intrepid2_Basis.hpp:125
Intrepid2::Impl::Basis_HCURL_TET_In_FEM::Serial
See Intrepid2::Basis_HCURL_TET_In_FEM.
Definition: Intrepid2_HCURL_TET_In_FEM.hpp:74
Intrepid2_HGRAD_TET_Cn_FEM_ORTH.hpp
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition: Intrepid2_Types.hpp:102
Intrepid2::PointTools::getLatticeSize
static ordinal_type getLatticeSize(const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...
Definition: Intrepid2_PointToolsDef.hpp:35
Intrepid2::CellTools::getReferenceFaceTangents
static void getReferenceFaceTangents(RefFaceTanViewType refFaceTanU, RefFaceTanViewType refFaceTanV, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:249
Generated by   doxygen 1.8.5