Intrepid2
Intrepid2_HCURL_TET_In_FEMDef.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 
49 #ifndef __INTREPID2_HCURL_TET_IN_FEM_DEF_HPP__
50 #define __INTREPID2_HCURL_TET_IN_FEM_DEF_HPP__
51 
54 #include "Teuchos_SerialDenseMatrix.hpp"
55 
56 namespace Intrepid2 {
57 
58 // -------------------------------------------------------------------------------------
59 
60 namespace Impl {
61 
62 template<EOperator opType>
63 template<typename OutputViewType,
64 typename inputViewType,
65 typename workViewType,
66 typename vinvViewType>
67 KOKKOS_INLINE_FUNCTION
68 void
69 Basis_HCURL_TET_In_FEM::Serial<opType>::
70 getValues( OutputViewType output,
71  const inputViewType input,
72  workViewType work,
73  const vinvViewType coeffs ) {
74 
75  constexpr ordinal_type spaceDim = 3;
76  const ordinal_type
77  cardPn = coeffs.extent(0)/spaceDim,
78  card = coeffs.extent(1),
79  npts = input.extent(0);
80 
81  // compute order
82  ordinal_type order = 0;
83  for (ordinal_type p=0;p<=Parameters::MaxOrder;++p) {
84  if (card == CardinalityHCurlTet(p)) {
85  order = p;
86  break;
87  }
88  }
89 
90  typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
91  auto vcprop = Kokkos::common_view_alloc_prop(work);
92  auto ptr = work.data();
93 
94  switch (opType) {
95  case OPERATOR_VALUE: {
96  const viewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts);
97  workViewType dummyView;
98 
99  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
100  Serial<opType>::getValues(phis, input, dummyView, order);
101 
102  for (ordinal_type i=0;i<card;++i)
103  for (ordinal_type j=0;j<npts;++j)
104  for (ordinal_type d=0;d<spaceDim;++d) {
105  output.access(i,j,d) = 0.0;
106  for (ordinal_type k=0;k<cardPn;++k)
107  output.access(i,j,d) += coeffs(k+d*cardPn,i) * phis(k,j);
108  }
109  break;
110  }
111  case OPERATOR_CURL: {
112  const viewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts, spaceDim);
113  ptr += card*npts*spaceDim*get_dimension_scalar(work);
114  const viewType workView(Kokkos::view_wrap(ptr, vcprop), card, npts, spaceDim+1);
115 
116  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
117  Serial<OPERATOR_GRAD>::getValues(phis, input, workView, order);
118 
119  for (ordinal_type i=0;i<card;++i) {
120  for (ordinal_type j=0;j<npts;++j) {
121  for (ordinal_type d=0; d< spaceDim; ++d) {
122  output.access(i,j,d) = 0.0;
123  ordinal_type d1 = (d+1) % spaceDim, d2 = (d+2) % spaceDim;
124  for (ordinal_type k=0; k<cardPn; ++k) //\sum_k (coeffs_k, coeffs_{k+cardPn}, coeffs_{k+2 cardPn}) \times phis_kj (cross product)
125  output.access(i,j,d) += coeffs(k+d2*cardPn,i)*phis(k,j,d1)
126  -coeffs(k+d1*cardPn,i)*phis(k,j,d2);
127  }
128  }
129  }
130  break;
131  }
132  default: {
133  INTREPID2_TEST_FOR_ABORT( true,
134  ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented");
135  }
136  }
137 }
138 
139 template<typename SpT, ordinal_type numPtsPerEval,
140 typename outputValueValueType, class ...outputValueProperties,
141 typename inputPointValueType, class ...inputPointProperties,
142 typename vinvValueType, class ...vinvProperties>
143 void
144 Basis_HCURL_TET_In_FEM::
145 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
146  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
147  const Kokkos::DynRankView<vinvValueType, vinvProperties...> coeffs,
148  const EOperator operatorType) {
149  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
150  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
151  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
152  typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
153 
154  // loopSize corresponds to cardinality
155  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
156  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
157  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
158  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
159 
160  typedef typename inputPointViewType::value_type inputPointType;
161 
162  const ordinal_type cardinality = outputValues.extent(0);
163  const ordinal_type spaceDim = 3;
164 
165  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
166  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
167 
168  switch (operatorType) {
169  case OPERATOR_VALUE: {
170  workViewType work(Kokkos::view_alloc("Basis_HCURL_TET_In_FEM::getValues::work", vcprop), cardinality, inputPoints.extent(0));
171  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
172  OPERATOR_VALUE,numPtsPerEval> FunctorType;
173  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, coeffs, work) );
174  break;
175  }
176  case OPERATOR_CURL: {
177  workViewType work(Kokkos::view_alloc("Basis_HCURL_TET_In_FEM::getValues::work", vcprop), cardinality*(2*spaceDim+1), inputPoints.extent(0));
178  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
179  OPERATOR_CURL,numPtsPerEval> FunctorType;
180  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, coeffs, work) );
181  break;
182  }
183  default: {
184  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
185  ">>> ERROR (Basis_HCURL_TET_In_FEM): Operator type not implemented" );
186  }
187  }
188 }
189 }
190 
191 // -------------------------------------------------------------------------------------
192 template<typename SpT, typename OT, typename PT>
194 Basis_HCURL_TET_In_FEM( const ordinal_type order,
195  const EPointType pointType ) {
196 
197  constexpr ordinal_type spaceDim = 3;
198  this->basisCardinality_ = CardinalityHCurlTet(order);
199  this->basisDegree_ = order; // small n
200  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
201  this->basisType_ = BASIS_FEM_FIAT;
202  this->basisCoordinates_ = COORDINATES_CARTESIAN;
203  this->functionSpace_ = FUNCTION_SPACE_HCURL;
204 
205  const ordinal_type card = this->basisCardinality_;
206 
207  const ordinal_type cardPn = Intrepid2::getPnCardinality<spaceDim>(order); // dim of (P_{n}) -- smaller space
208  const ordinal_type cardPnm1 = Intrepid2::getPnCardinality<spaceDim>(order-1); // dim of (P_{n-1}) -- smaller space
209  const ordinal_type cardPnm2 = Intrepid2::getPnCardinality<spaceDim>(order-2); // dim of (P_{n-2}) -- smaller space
210  const ordinal_type cardVecPn = spaceDim*cardPn; // dim of (P_{n})^2 -- larger space
211  const ordinal_type cardVecPnm1 = spaceDim*cardPnm1; // dim of (P_{n-1})^2 -- smaller space
212  const ordinal_type cardPnm1H = cardPnm1-cardPnm2; //Homogeneous polynomial of order (n-1)
213 
214 
215 
216  // Basis-dependent initializations
217  constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
218  constexpr ordinal_type maxCard = CardinalityHCurlTet(Parameters::MaxOrder);
219  ordinal_type tags[maxCard][tagSize];
220 
221  // points are computed in the host and will be copied
222  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
223  dofCoords("Hcurl::Tet::In::dofCoords", card, spaceDim);
224 
225  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
226  coeffs("Hcurl::Tet::In::coeffs", cardVecPn, card);
227 
228  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
229  dofCoeffs("Hcurl::Tet::In::dofCoeffs", card, spaceDim);
230 
231  // first, need to project the basis for RT space onto the
232  // orthogonal basis of degree n
233  // get coefficients of PkHx
234 
235  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace> //use LayoutLeft for Lapack
236  V1("Hcurl::Tet::In::V1", cardVecPn, cardVecPnm1 + spaceDim*cardPnm1H);
237 
238 
239  // these two loops get the first three sets of basis functions
240  for (ordinal_type i=0;i<cardPnm1;i++)
241  for (ordinal_type d=0;d<spaceDim;d++)
242  V1(i+d*cardPn,i+d*cardPnm1) = 1.0;
243 
244 
245  // now I need to integrate { (x,y) \times phi } against the big basis
246  // first, get a cubature rule.
248  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> cubPoints("Hcurl::Tet::In::cubPoints", myCub.getNumPoints() , spaceDim );
249  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> cubWeights("Hcurl::Tet::In::cubWeights", myCub.getNumPoints() );
250  myCub.getCubature( cubPoints , cubWeights );
251 
252  // tabulate the scalar orthonormal basis at cubature points
253  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> phisAtCubPoints("Hcurl::Tet::In::phisAtCubPoints", cardPn , myCub.getNumPoints() );
254  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtCubPoints, cubPoints, order, OPERATOR_VALUE);
255 
256  // 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)
257  for (ordinal_type i=0;i<cardPn;i++) {
258  for (ordinal_type j=0;j<cardPnm1H;j++) { // loop over homogeneous polynomials
259  for (ordinal_type d=0; d< spaceDim; ++d) {
260  scalarType integral(0);
261  for (ordinal_type k=0;k<myCub.getNumPoints();k++)
262  integral += cubWeights(k) * cubPoints(k,d)
263  * phisAtCubPoints(cardPnm2+j,k)
264  * phisAtCubPoints(i,k);
265  ordinal_type d1 = (d+1) % spaceDim, d2 = (d+2) % spaceDim;
266  V1(i+d2*cardPn,cardVecPnm1+d1*cardPnm1H + j) = -integral;
267  V1(i+d1*cardPn,cardVecPnm1+d2*cardPnm1H + j) = integral;
268  }
269  }
270  }
271 
272 
273 
274 
275 
276  // now I need to set up an SVD to get a basis for the space
277  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
278  S("Hcurl::Tet::In::S", cardVecPn,1),
279  U("Hcurl::Tet::In::U", cardVecPn, cardVecPn),
280  Vt("Hcurl::Tet::In::Vt", cardVecPn, cardVecPn),
281  work("Hcurl::Tet::In::work", 5*cardVecPn,1),
282  rWork("Hcurl::Tet::In::rW", 1,1);
283 
284 
285 
286  ordinal_type info = 0;
287  Teuchos::LAPACK<ordinal_type,scalarType> lapack;
288 
289 
290  lapack.GESVD( 'A',
291  'N',
292  V1.extent(0) ,
293  V1.extent(1) ,
294  V1.data() ,
295  V1.stride_1() ,
296  S.data() ,
297  U.data() ,
298  U.stride_1() ,
299  Vt.data() ,
300  Vt.stride_1() ,
301  work.data() ,
302  5*cardVecPn ,
303  rWork.data() ,
304  &info );
305 
306 
307 #ifdef HAVE_INTREPID2_DEBUG
308  ordinal_type num_nonzero_sv = 0;
309  for (int i=0;i<cardVecPn;i++)
310  num_nonzero_sv += (S(i,0) > tolerence());
311 
312  INTREPID2_TEST_FOR_EXCEPTION( num_nonzero_sv != card, std::invalid_argument,
313  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM( order, pointType), Matrix V1 should have rank equal to the cardinality of HCURL space");
314 #endif
315 
316  // next, apply the RT nodes (rows) to the basis for (P_n)^2 (columns)
317  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
318  V2("Hcurl::Tet::In::V2", card ,cardVecPn);
319 
320  const ordinal_type numEdges = this->basisCellTopology_.getEdgeCount();
321  const ordinal_type numFaces = this->basisCellTopology_.getFaceCount();
322 
323  // first numEdges * degree nodes are normals at each edge
324  // get the points on the line
325 
326  shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
327  shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
328 
329  const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
330  order+1 ,
331  1 );
332 
333  const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
334  order+1 ,
335  1 );
336 
337  const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
338  order+1 ,
339  1 );
340 
341  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> linePts("Hcurl::Tet::In::linePts", numPtsPerEdge , 1 );
342  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> triPts("Hcurl::Tet::In::triPts", numPtsPerFace , 2 );
343 
344  // construct lattice
345  const ordinal_type offset = 1;
346 
347 
348 
349  PointTools::getLattice( linePts,
350  edgeTop,
351  order+1, offset,
352  pointType );
353 
354  PointTools::getLattice( triPts,
355  faceTop,
356  order+1, offset,
357  pointType );
358 
359  // holds the image of the line points
360  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> edgePts("Hcurl::Tet::In::edgePts", numPtsPerEdge , spaceDim );
361  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> facePts("Hcurl::Tet::In::facePts", numPtsPerFace , spaceDim );
362  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> phisAtEdgePoints("Hcurl::Tet::In::phisAtEdgePoints", cardPn , numPtsPerEdge );
363  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> phisAtFacePoints("Hcurl::Tet::In::phisAtFacePoints", cardPn , numPtsPerFace);
364 
365  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> edgeTan("Hcurl::Tet::In::edgeTan", spaceDim );
366 
367  // these are tangents scaled by the appropriate edge lengths.
368  for (ordinal_type i=0;i<numEdges;i++) { // loop over edges
370  i ,
371  this->basisCellTopology_ );
372  /* multiply by measure of reference edge so that magnitude of the edgeTan is equal to the edge measure */
373  const scalarType refEdgeMeasure = 2.0;
374  for (ordinal_type j=0;j<spaceDim;j++)
375  edgeTan(j) *= refEdgeMeasure;
376 
378  linePts ,
379  1 ,
380  i ,
381  this->basisCellTopology_ );
382 
383  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtEdgePoints , edgePts, order, OPERATOR_VALUE);
384 
385  // loop over points (rows of V2)
386  for (ordinal_type j=0;j<numPtsPerEdge;j++) {
387 
388  const ordinal_type i_card = numPtsPerEdge*i+j;
389 
390  // loop over orthonormal basis functions (columns of V2)
391  for (ordinal_type k=0;k<cardPn;k++)
392  for (ordinal_type d=0;d<spaceDim;d++)
393  V2(i_card,k+d*cardPn) = edgeTan(d) * phisAtEdgePoints(k,j);
394 
395  //save dof coordinates and coefficients
396  for(ordinal_type k=0; k<spaceDim; ++k) {
397  dofCoords(i_card,k) = edgePts(j,k);
398  dofCoeffs(i_card,k) = edgeTan(k);
399  }
400 
401  tags[i_card][0] = 1; // edge dof
402  tags[i_card][1] = i; // edge id
403  tags[i_card][2] = j; // local dof id
404  tags[i_card][3] = numPtsPerEdge; // total vert dof
405 
406  }
407  }
408 
409  if(numPtsPerFace >0) {//handle faces if needed (order >1)
410  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> faceTan1("Hcurl::Tet::In::edgeTan", spaceDim );
411  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> faceTan2("Hcurl::Tet::In::edgeTan", spaceDim );
412 
413  for (ordinal_type i=0;i<numFaces;i++) { // loop over faces
415  faceTan2,
416  i ,
417  this->basisCellTopology_ );
418 
420  triPts ,
421  2 ,
422  i ,
423  this->basisCellTopology_ );
424 
425  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtFacePoints , facePts, order, OPERATOR_VALUE);
426 
427  // loop over points (rows of V2)
428  for (ordinal_type j=0;j<numPtsPerFace;j++) {
429 
430  const ordinal_type i_card = numEdges*numPtsPerEdge+2*numPtsPerFace*i+2*j;
431  const ordinal_type i_card_p1 = i_card+1; // creating a temp otherwise nvcc gets confused
432 
433  // loop over orthonormal basis functions (columns of V2)
434  for (ordinal_type k=0;k<cardPn;k++)
435  for (ordinal_type d=0;d<spaceDim;d++) {
436  V2(i_card,k+d*cardPn) = faceTan1(d) * phisAtFacePoints(k,j);
437  V2(i_card_p1,k+d*cardPn) = faceTan2(d) * phisAtFacePoints(k,j);
438  }
439 
440  //save dof coordinates
441  for(ordinal_type k=0; k<spaceDim; ++k) {
442  dofCoords(i_card,k) = facePts(j,k);
443  dofCoords(i_card_p1,k) = facePts(j,k);
444  dofCoeffs(i_card,k) = faceTan1(k);
445  dofCoeffs(i_card_p1,k) = faceTan2(k);
446  }
447 
448  tags[i_card][0] = 2; // face dof
449  tags[i_card][1] = i; // face id
450  tags[i_card][2] = 2*j; // local face id
451  tags[i_card][3] = 2*numPtsPerFace; // total face dof
452 
453  tags[i_card_p1][0] = 2; // face dof
454  tags[i_card_p1][1] = i; // face id
455  tags[i_card_p1][2] = 2*j+1; // local face id
456  tags[i_card_p1][3] = 2*numPtsPerFace; // total face dof
457 
458  }
459  }
460  }
461 
462 
463  // internal dof, if needed
464  if (numPtsPerCell > 0) {
465  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
466  cellPoints( "Hcurl::Tet::In::cellPoints", numPtsPerCell , spaceDim );
467  PointTools::getLattice( cellPoints ,
468  this->basisCellTopology_ ,
469  order + 1 ,
470  1 ,
471  pointType );
472 
473  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
474  phisAtCellPoints("Hcurl::Tet::In::phisAtCellPoints", cardPn , numPtsPerCell );
475  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>( phisAtCellPoints , cellPoints , order, OPERATOR_VALUE );
476 
477  // copy values into right positions of V2
478  for (ordinal_type j=0;j<numPtsPerCell;j++) {
479 
480  const ordinal_type i_card = numEdges*numPtsPerEdge+2*numFaces*numPtsPerFace+spaceDim*j;
481 
482  for (ordinal_type k=0;k<cardPn;k++)
483  for (ordinal_type d=0;d<spaceDim;d++)
484  V2(i_card+d,d*cardPn+k) = phisAtCellPoints(k,j);
485 
486 
487  //save dof coordinates
488  for(ordinal_type d=0; d<spaceDim; ++d) {
489  for(ordinal_type dim=0; dim<spaceDim; ++dim) {
490  dofCoords(i_card+d,dim) = cellPoints(j,dim);
491  dofCoeffs(i_card+d,dim) = (d==dim);
492  }
493 
494  tags[i_card+d][0] = spaceDim; // elem dof
495  tags[i_card+d][1] = 0; // elem id
496  tags[i_card+d][2] = spaceDim*j+d; // local dof id
497  tags[i_card+d][3] = spaceDim*numPtsPerCell; // total vert dof
498  }
499  }
500  }
501 
502  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
503  // so we transpose on copy below.
504  const ordinal_type lwork = card*card;
505  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
506  vmat("Hcurl::Tet::In::vmat", card, card),
507  work1("Hcurl::Tet::In::work", lwork),
508  ipiv("Hcurl::Tet::In::ipiv", card);
509 
510  //vmat = V2*U;
511  for(ordinal_type i=0; i< card; ++i) {
512  for(ordinal_type j=0; j< card; ++j) {
513  scalarType s=0;
514  for(ordinal_type k=0; k< cardVecPn; ++k)
515  s += V2(i,k)*U(k,j);
516  vmat(i,j) = s;
517  }
518  }
519 
520  info = 0;
521 
522  lapack.GETRF(card, card,
523  vmat.data(), vmat.stride_1(),
524  (ordinal_type*)ipiv.data(),
525  &info);
526 
527  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
528  std::runtime_error ,
529  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRF returns nonzero info." );
530 
531  lapack.GETRI(card,
532  vmat.data(), vmat.stride_1(),
533  (ordinal_type*)ipiv.data(),
534  work1.data(), lwork,
535  &info);
536 
537  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
538  std::runtime_error ,
539  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRI returns nonzero info." );
540 
541  for (ordinal_type i=0;i<cardVecPn;++i) {
542  for (ordinal_type j=0;j<card;++j){
543  scalarType s=0;
544  for(ordinal_type k=0; k< card; ++k)
545  s += U(i,k)*vmat(k,j);
546  coeffs(i,j) = s;
547  }
548  }
549 
550  this->coeffs_ = Kokkos::create_mirror_view(typename SpT::memory_space(), coeffs);
551  Kokkos::deep_copy(this->coeffs_ , coeffs);
552 
553  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
554  Kokkos::deep_copy(this->dofCoords_, dofCoords);
555 
556  this->dofCoeffs_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoeffs);
557  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
558 
559 
560  // set tags
561  {
562  // Basis-dependent initializations
563  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
564  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
565  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
566 
567  OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
568 
569  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
570  // tags are constructed on host
571  this->setOrdinalTagData(this->tagToOrdinal_,
572  this->ordinalToTag_,
573  tagView,
574  this->basisCardinality_,
575  tagSize,
576  posScDim,
577  posScOrd,
578  posDfOrd);
579  }
580 }
581 } // namespace Intrepid2
582 #endif
static void mapToReferenceSubcell(Kokkos::DynRankView< refSubcellPointValueType, refSubcellPointProperties... > refSubcellPoints, const Kokkos::DynRankView< paramPointValueType, paramPointProperties... > paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
virtual ordinal_type getNumPoints() const
Returns the number of cubature points.
Kokkos::View< ordinal_type *, typename ExecSpaceType::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
static void getReferenceEdgeTangent(Kokkos::DynRankView< refEdgeTangentValueType, refEdgeTangentProperties... > refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
static void getReferenceFaceTangents(Kokkos::DynRankView< refFaceTanUValueType, refFaceTanUProperties... > refFaceTanU, Kokkos::DynRankView< refFaceTanVValueType, refFaceTanVProperties... > refFaceTanV, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
Basis_HCURL_TET_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Header file for the Intrepid2::CubatureDirectTetDefault class.
EPointType
Enumeration of types of point distributions in Intrepid.
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 (currently disabled for other ce...
Defines direct integration rules on a tetrahedron.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
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...