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