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
23namespace Intrepid2 {
24
25// -------------------------------------------------------------------------------------
26
27namespace Impl {
28
29template<EOperator OpType>
30template<typename OutputViewType,
31typename InputViewType,
32typename WorkViewType,
33typename VinvViewType>
34KOKKOS_INLINE_FUNCTION
35void
36Basis_HCURL_TET_In_FEM::Serial<OpType>::
37getValues( 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
106template<typename DT, ordinal_type numPtsPerEval,
107typename outputValueValueType, class ...outputValueProperties,
108typename inputPointValueType, class ...inputPointProperties,
109typename vinvValueType, class ...vinvProperties>
110void
111Basis_HCURL_TET_In_FEM::
112getValues(
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// -------------------------------------------------------------------------------------
161template<typename DT, typename OT, typename PT>
162Basis_HCURL_TET_In_FEM<DT,OT,PT>::
163Basis_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
566template<typename DT, typename OT, typename PT>
567void
568Basis_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
578template<typename DT, typename OT, typename PT>
579KOKKOS_INLINE_FUNCTION
580void
581Basis_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), [=, &coeffs_ = this->coeffs_] (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, coeffs_ );
608 });
609 break;
610 case OPERATOR_CURL:
611 Kokkos::parallel_for (Kokkos::TeamThreadRange (team_member, numPoints), [=, &coeffs_ = this->coeffs_] (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, 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_CubatureDirectTetDefault.hpp
Header file for the Intrepid2::CubatureDirectTetDefault class.
Intrepid2_HGRAD_TET_Cn_FEM_ORTH.hpp
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
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::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::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::Basis::PointViewType
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Definition Intrepid2_Basis.hpp:338
Intrepid2::Basis< DeviceType, outputValueType, pointValueType >::scalarType
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
Definition Intrepid2_Basis.hpp:157
Intrepid2::Basis::OrdinalTypeArray1DHost
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
Definition Intrepid2_Basis.hpp:125
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::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::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
Intrepid2::CubatureDirectTetDefault
Defines direct integration rules on a tetrahedron.
Definition Intrepid2_CubatureDirectTetDefault.hpp:31
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::CubatureDirect::getNumPoints
virtual ordinal_type getNumPoints() const override
Returns the number of cubature points.
Definition Intrepid2_CubatureDirect.hpp:164
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::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::Impl::Basis_HCURL_TET_In_FEM::Serial
See Intrepid2::Basis_HCURL_TET_In_FEM.
Definition Intrepid2_HCURL_TET_In_FEM.hpp:74
Generated by doxygen 1.9.8