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 
204  const ordinal_type card = this->basisCardinality_;
205 
206  const ordinal_type cardPn = Intrepid2::getPnCardinality<spaceDim>(order); // dim of (P_{n}) -- smaller space
207  const ordinal_type cardPnm1 = Intrepid2::getPnCardinality<spaceDim>(order-1); // dim of (P_{n-1}) -- smaller space
208  const ordinal_type cardPnm2 = Intrepid2::getPnCardinality<spaceDim>(order-2); // dim of (P_{n-2}) -- smaller space
209  const ordinal_type cardVecPn = spaceDim*cardPn; // dim of (P_{n})^2 -- larger space
210  const ordinal_type cardVecPnm1 = spaceDim*cardPnm1; // dim of (P_{n-1})^2 -- smaller space
211  const ordinal_type cardPnm1H = cardPnm1-cardPnm2; //Homogeneous polynomial of order (n-1)
212 
213 
214 
215  // Basis-dependent initializations
216  constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
217  constexpr ordinal_type maxCard = CardinalityHCurlTet(Parameters::MaxOrder);
218  ordinal_type tags[maxCard][tagSize];
219 
220  // points are computed in the host and will be copied
221  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
222  dofCoords("Hcurl::Tet::In::dofCoords", card, spaceDim);
223 
224  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
225  coeffs("Hcurl::Tet::In::coeffs", cardVecPn, card);
226 
227  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
228  dofCoeffs("Hcurl::Tet::In::dofCoeffs", card, spaceDim);
229 
230  // first, need to project the basis for RT space onto the
231  // orthogonal basis of degree n
232  // get coefficients of PkHx
233 
234  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace> //use LayoutLeft for Lapack
235  V1("Hcurl::Tet::In::V1", cardVecPn, cardVecPnm1 + spaceDim*cardPnm1H);
236 
237 
238  // these two loops get the first three sets of basis functions
239  for (ordinal_type i=0;i<cardPnm1;i++)
240  for (ordinal_type d=0;d<spaceDim;d++)
241  V1(i+d*cardPn,i+d*cardPnm1) = 1.0;
242 
243 
244  // now I need to integrate { (x,y) \times phi } against the big basis
245  // first, get a cubature rule.
247  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> cubPoints("Hcurl::Tet::In::cubPoints", myCub.getNumPoints() , spaceDim );
248  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> cubWeights("Hcurl::Tet::In::cubWeights", myCub.getNumPoints() );
249  myCub.getCubature( cubPoints , cubWeights );
250 
251  // tabulate the scalar orthonormal basis at cubature points
252  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> phisAtCubPoints("Hcurl::Tet::In::phisAtCubPoints", cardPn , myCub.getNumPoints() );
253  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtCubPoints, cubPoints, order, OPERATOR_VALUE);
254 
255  // 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)
256  for (ordinal_type i=0;i<cardPn;i++) {
257  for (ordinal_type j=0;j<cardPnm1H;j++) { // loop over homogeneous polynomials
258  for (ordinal_type d=0; d< spaceDim; ++d) {
259  scalarType integral(0);
260  for (ordinal_type k=0;k<myCub.getNumPoints();k++)
261  integral += cubWeights(k) * cubPoints(k,d)
262  * phisAtCubPoints(cardPnm2+j,k)
263  * phisAtCubPoints(i,k);
264  ordinal_type d1 = (d+1) % spaceDim, d2 = (d+2) % spaceDim;
265  V1(i+d2*cardPn,cardVecPnm1+d1*cardPnm1H + j) = -integral;
266  V1(i+d1*cardPn,cardVecPnm1+d2*cardPnm1H + j) = integral;
267  }
268  }
269  }
270 
271 
272 
273 
274 
275  // now I need to set up an SVD to get a basis for the space
276  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
277  S("Hcurl::Tet::In::S", cardVecPn,1),
278  U("Hcurl::Tet::In::U", cardVecPn, cardVecPn),
279  Vt("Hcurl::Tet::In::Vt", cardVecPn, cardVecPn),
280  work("Hcurl::Tet::In::work", 5*cardVecPn,1),
281  rWork("Hcurl::Tet::In::rW", 1,1);
282 
283 
284 
285  ordinal_type info = 0;
286  Teuchos::LAPACK<ordinal_type,scalarType> lapack;
287 
288 
289  lapack.GESVD( 'A',
290  'N',
291  V1.extent(0) ,
292  V1.extent(1) ,
293  V1.data() ,
294  V1.stride_1() ,
295  S.data() ,
296  U.data() ,
297  U.stride_1() ,
298  Vt.data() ,
299  Vt.stride_1() ,
300  work.data() ,
301  5*cardVecPn ,
302  rWork.data() ,
303  &info );
304 
305 
306 #ifdef HAVE_INTREPID2_DEBUG
307  ordinal_type num_nonzero_sv = 0;
308  for (int i=0;i<cardVecPn;i++)
309  num_nonzero_sv += (S(i,0) > tolerence());
310 
311  INTREPID2_TEST_FOR_EXCEPTION( num_nonzero_sv != card, std::invalid_argument,
312  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM( order, pointType), Matrix V1 should have rank equal to the cardinality of HCURL space");
313 #endif
314 
315  // next, apply the RT nodes (rows) to the basis for (P_n)^2 (columns)
316  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
317  V2("Hcurl::Tet::In::V2", card ,cardVecPn);
318 
319  const ordinal_type numEdges = this->basisCellTopology_.getEdgeCount();
320  const ordinal_type numFaces = this->basisCellTopology_.getFaceCount();
321 
322  // first numEdges * degree nodes are normals at each edge
323  // get the points on the line
324 
325  shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
326  shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
327 
328  const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
329  order+1 ,
330  1 );
331 
332  const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
333  order+1 ,
334  1 );
335 
336  const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
337  order+1 ,
338  1 );
339 
340  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> linePts("Hcurl::Tet::In::linePts", numPtsPerEdge , 1 );
341  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> triPts("Hcurl::Tet::In::triPts", numPtsPerFace , 2 );
342 
343  // construct lattice
344  const ordinal_type offset = 1;
345 
346 
347 
348  PointTools::getLattice( linePts,
349  edgeTop,
350  order+1, offset,
351  pointType );
352 
353  PointTools::getLattice( triPts,
354  faceTop,
355  order+1, offset,
356  pointType );
357 
358  // holds the image of the line points
359  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> edgePts("Hcurl::Tet::In::edgePts", numPtsPerEdge , spaceDim );
360  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> facePts("Hcurl::Tet::In::facePts", numPtsPerFace , spaceDim );
361  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> phisAtEdgePoints("Hcurl::Tet::In::phisAtEdgePoints", cardPn , numPtsPerEdge );
362  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> phisAtFacePoints("Hcurl::Tet::In::phisAtFacePoints", cardPn , numPtsPerFace);
363 
364  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> edgeTan("Hcurl::Tet::In::edgeTan", spaceDim );
365 
366  // these are tangents scaled by the appropriate edge lengths.
367  for (ordinal_type i=0;i<numEdges;i++) { // loop over edges
369  i ,
370  this->basisCellTopology_ );
371  /* multiply by measure of reference edge so that magnitude of the edgeTan is equal to the edge measure */
372  const scalarType refEdgeMeasure = 2.0;
373  for (ordinal_type j=0;j<spaceDim;j++)
374  edgeTan(j) *= refEdgeMeasure;
375 
377  linePts ,
378  1 ,
379  i ,
380  this->basisCellTopology_ );
381 
382  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtEdgePoints , edgePts, order, OPERATOR_VALUE);
383 
384  // loop over points (rows of V2)
385  for (ordinal_type j=0;j<numPtsPerEdge;j++) {
386 
387  const ordinal_type i_card = numPtsPerEdge*i+j;
388 
389  // loop over orthonormal basis functions (columns of V2)
390  for (ordinal_type k=0;k<cardPn;k++)
391  for (ordinal_type d=0;d<spaceDim;d++)
392  V2(i_card,k+d*cardPn) = edgeTan(d) * phisAtEdgePoints(k,j);
393 
394  //save dof coordinates and coefficients
395  for(ordinal_type k=0; k<spaceDim; ++k) {
396  dofCoords(i_card,k) = edgePts(j,k);
397  dofCoeffs(i_card,k) = edgeTan(k);
398  }
399 
400  tags[i_card][0] = 1; // edge dof
401  tags[i_card][1] = i; // edge id
402  tags[i_card][2] = j; // local dof id
403  tags[i_card][3] = numPtsPerEdge; // total vert dof
404 
405  }
406  }
407 
408  if(numPtsPerFace >0) {//handle faces if needed (order >1)
409  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> faceTan1("Hcurl::Tet::In::edgeTan", spaceDim );
410  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace> faceTan2("Hcurl::Tet::In::edgeTan", spaceDim );
411 
412  for (ordinal_type i=0;i<numFaces;i++) { // loop over faces
414  faceTan2,
415  i ,
416  this->basisCellTopology_ );
417 
419  triPts ,
420  2 ,
421  i ,
422  this->basisCellTopology_ );
423 
424  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtFacePoints , facePts, order, OPERATOR_VALUE);
425 
426  // loop over points (rows of V2)
427  for (ordinal_type j=0;j<numPtsPerFace;j++) {
428 
429  const ordinal_type i_card = numEdges*numPtsPerEdge+2*numPtsPerFace*i+2*j;
430  const ordinal_type i_card_p1 = i_card+1; // creating a temp otherwise nvcc gets confused
431 
432  // loop over orthonormal basis functions (columns of V2)
433  for (ordinal_type k=0;k<cardPn;k++)
434  for (ordinal_type d=0;d<spaceDim;d++) {
435  V2(i_card,k+d*cardPn) = faceTan1(d) * phisAtFacePoints(k,j);
436  V2(i_card_p1,k+d*cardPn) = faceTan2(d) * phisAtFacePoints(k,j);
437  }
438 
439  //save dof coordinates
440  for(ordinal_type k=0; k<spaceDim; ++k) {
441  dofCoords(i_card,k) = facePts(j,k);
442  dofCoords(i_card_p1,k) = facePts(j,k);
443  dofCoeffs(i_card,k) = faceTan1(k);
444  dofCoeffs(i_card_p1,k) = faceTan2(k);
445  }
446 
447  tags[i_card][0] = 2; // face dof
448  tags[i_card][1] = i; // face id
449  tags[i_card][2] = 2*j; // local face id
450  tags[i_card][3] = 2*numPtsPerFace; // total face dof
451 
452  tags[i_card_p1][0] = 2; // face dof
453  tags[i_card_p1][1] = i; // face id
454  tags[i_card_p1][2] = 2*j+1; // local face id
455  tags[i_card_p1][3] = 2*numPtsPerFace; // total face dof
456 
457  }
458  }
459  }
460 
461 
462  // internal dof, if needed
463  if (numPtsPerCell > 0) {
464  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
465  cellPoints( "Hcurl::Tet::In::cellPoints", numPtsPerCell , spaceDim );
466  PointTools::getLattice( cellPoints ,
467  this->basisCellTopology_ ,
468  order + 1 ,
469  1 ,
470  pointType );
471 
472  Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
473  phisAtCellPoints("Hcurl::Tet::In::phisAtCellPoints", cardPn , numPtsPerCell );
474  Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>( phisAtCellPoints , cellPoints , order, OPERATOR_VALUE );
475 
476  // copy values into right positions of V2
477  for (ordinal_type j=0;j<numPtsPerCell;j++) {
478 
479  const ordinal_type i_card = numEdges*numPtsPerEdge+2*numFaces*numPtsPerFace+spaceDim*j;
480 
481  for (ordinal_type k=0;k<cardPn;k++)
482  for (ordinal_type d=0;d<spaceDim;d++)
483  V2(i_card+d,d*cardPn+k) = phisAtCellPoints(k,j);
484 
485 
486  //save dof coordinates
487  for(ordinal_type d=0; d<spaceDim; ++d) {
488  for(ordinal_type dim=0; dim<spaceDim; ++dim) {
489  dofCoords(i_card+d,dim) = cellPoints(j,dim);
490  dofCoeffs(i_card+d,dim) = (d==dim);
491  }
492 
493  tags[i_card+d][0] = spaceDim; // elem dof
494  tags[i_card+d][1] = 0; // elem id
495  tags[i_card+d][2] = spaceDim*j+d; // local dof id
496  tags[i_card+d][3] = spaceDim*numPtsPerCell; // total vert dof
497  }
498  }
499  }
500 
501  // form Vandermonde matrix. Actually, this is the transpose of the VDM,
502  // so we transpose on copy below.
503  const ordinal_type lwork = card*card;
504  Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
505  vmat("Hcurl::Tet::In::vmat", card, card),
506  work1("Hcurl::Tet::In::work", lwork),
507  ipiv("Hcurl::Tet::In::ipiv", card);
508 
509  //vmat = V2*U;
510  for(ordinal_type i=0; i< card; ++i) {
511  for(ordinal_type j=0; j< card; ++j) {
512  scalarType s=0;
513  for(ordinal_type k=0; k< cardVecPn; ++k)
514  s += V2(i,k)*U(k,j);
515  vmat(i,j) = s;
516  }
517  }
518 
519  info = 0;
520 
521  lapack.GETRF(card, card,
522  vmat.data(), vmat.stride_1(),
523  (ordinal_type*)ipiv.data(),
524  &info);
525 
526  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
527  std::runtime_error ,
528  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRF returns nonzero info." );
529 
530  lapack.GETRI(card,
531  vmat.data(), vmat.stride_1(),
532  (ordinal_type*)ipiv.data(),
533  work1.data(), lwork,
534  &info);
535 
536  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
537  std::runtime_error ,
538  ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRI returns nonzero info." );
539 
540  for (ordinal_type i=0;i<cardVecPn;++i) {
541  for (ordinal_type j=0;j<card;++j){
542  scalarType s=0;
543  for(ordinal_type k=0; k< card; ++k)
544  s += U(i,k)*vmat(k,j);
545  coeffs(i,j) = s;
546  }
547  }
548 
549  this->coeffs_ = Kokkos::create_mirror_view(typename SpT::memory_space(), coeffs);
550  Kokkos::deep_copy(this->coeffs_ , coeffs);
551 
552  this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
553  Kokkos::deep_copy(this->dofCoords_, dofCoords);
554 
555  this->dofCoeffs_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoeffs);
556  Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
557 
558 
559  // set tags
560  {
561  // Basis-dependent initializations
562  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
563  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
564  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
565 
566  ordinal_type_array_1d_host tagView(&tags[0][0], card*tagSize);
567 
568  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
569  // tags are constructed on host
570  this->setOrdinalTagData(this->tagToOrdinal_,
571  this->ordinalToTag_,
572  tagView,
573  this->basisCardinality_,
574  tagSize,
575  posScDim,
576  posScOrd,
577  posDfOrd);
578  }
579 }
580 } // namespace Intrepid2
581 #endif
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.
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...
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.
Basis_HCURL_TET_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
Header file for the Intrepid2::CubatureDirectTetDefault class.
virtual void getCubature(pointViewType cubPoints, weightViewType cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Defines direct integration rules on a tetrahedron.
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.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
virtual ordinal_type getNumPoints() const
Returns the number of cubature points.
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...