195 const EPointType pointType ) {
197 constexpr ordinal_type spaceDim = 3;
198 this->basisCardinality_ = CardinalityHCurlTet(order);
199 this->basisDegree_ = order;
200 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
201 this->basisType_ = BASIS_FEM_LAGRANGIAN;
202 this->basisCoordinates_ = COORDINATES_CARTESIAN;
203 this->functionSpace_ = FUNCTION_SPACE_HCURL;
204 pointType_ = pointType;
205 const ordinal_type card = this->basisCardinality_;
207 const ordinal_type cardPn = Intrepid2::getPnCardinality<spaceDim>(order);
208 const ordinal_type cardPnm1 = Intrepid2::getPnCardinality<spaceDim>(order-1);
209 const ordinal_type cardPnm2 = Intrepid2::getPnCardinality<spaceDim>(order-2);
210 const ordinal_type cardVecPn = spaceDim*cardPn;
211 const ordinal_type cardVecPnm1 = spaceDim*cardPnm1;
212 const ordinal_type cardPnm1H = cardPnm1-cardPnm2;
216 INTREPID2_TEST_FOR_EXCEPTION( order >
Parameters::MaxOrder, std::invalid_argument,
"polynomial order exceeds the max supported by this class");
219 constexpr ordinal_type tagSize = 4;
221 ordinal_type tags[maxCard][tagSize];
224 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
225 dofCoords(
"Hcurl::Tet::In::dofCoords", card, spaceDim);
227 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
228 coeffs(
"Hcurl::Tet::In::coeffs", cardVecPn, card);
230 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
231 dofCoeffs(
"Hcurl::Tet::In::dofCoeffs", card, spaceDim);
237 Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
238 V1(
"Hcurl::Tet::In::V1", cardVecPn, cardVecPnm1 + spaceDim*cardPnm1H);
242 for (ordinal_type i=0;i<cardPnm1;i++)
243 for (ordinal_type d=0;d<spaceDim;d++)
244 V1(i+d*cardPn,i+d*cardPnm1) = 1.0;
250 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubPoints(
"Hcurl::Tet::In::cubPoints", myCub.
getNumPoints() , spaceDim );
251 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubWeights(
"Hcurl::Tet::In::cubWeights", myCub.
getNumPoints() );
255 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> phisAtCubPoints(
"Hcurl::Tet::In::phisAtCubPoints", cardPn , myCub.
getNumPoints() );
256 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(phisAtCubPoints, cubPoints, order, OPERATOR_VALUE);
259 for (ordinal_type i=0;i<cardPn;i++) {
260 for (ordinal_type j=0;j<cardPnm1H;j++) {
261 for (ordinal_type d=0; d< spaceDim; ++d) {
262 scalarType integral(0);
264 integral += cubWeights(k) * cubPoints(k,d)
265 * phisAtCubPoints(cardPnm2+j,k)
266 * phisAtCubPoints(i,k);
267 ordinal_type d1 = (d+1) % spaceDim, d2 = (d+2) % spaceDim;
268 V1(i+d2*cardPn,cardVecPnm1+d1*cardPnm1H + j) = -integral;
269 V1(i+d1*cardPn,cardVecPnm1+d2*cardPnm1H + j) = integral;
279 Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
280 S(
"Hcurl::Tet::In::S", cardVecPn,1),
281 U(
"Hcurl::Tet::In::U", cardVecPn, cardVecPn),
282 Vt(
"Hcurl::Tet::In::Vt", cardVecPn, cardVecPn),
283 work(
"Hcurl::Tet::In::work", 5*cardVecPn,1),
284 rWork(
"Hcurl::Tet::In::rW", 1,1);
288 ordinal_type info = 0;
289 Teuchos::LAPACK<ordinal_type,scalarType> lapack;
309#ifdef HAVE_INTREPID2_DEBUG
310 ordinal_type num_nonzero_sv = 0;
311 for (
int i=0;i<cardVecPn;i++)
312 num_nonzero_sv += (S(i,0) > tolerence());
314 INTREPID2_TEST_FOR_EXCEPTION( num_nonzero_sv != card, std::invalid_argument,
315 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM( order, pointType), Matrix V1 should have rank equal to the cardinality of HCURL space");
319 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
320 V2(
"Hcurl::Tet::In::V2", card ,cardVecPn);
322 const ordinal_type numEdges = this->basisCellTopology_.getEdgeCount();
323 const ordinal_type numFaces = this->basisCellTopology_.getFaceCount();
328 shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
329 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
343 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> linePts(
"Hcurl::Tet::In::linePts", numPtsPerEdge , 1 );
344 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> triPts(
"Hcurl::Tet::In::triPts", numPtsPerFace , 2 );
347 const ordinal_type offset = 1;
362 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> edgePts(
"Hcurl::Tet::In::edgePts", numPtsPerEdge , spaceDim );
363 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> facePts(
"Hcurl::Tet::In::facePts", numPtsPerFace , spaceDim );
364 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> phisAtEdgePoints(
"Hcurl::Tet::In::phisAtEdgePoints", cardPn , numPtsPerEdge );
365 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> phisAtFacePoints(
"Hcurl::Tet::In::phisAtFacePoints", cardPn , numPtsPerFace);
367 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> edgeTan(
"Hcurl::Tet::In::edgeTan", spaceDim );
370 for (ordinal_type i=0;i<numEdges;i++) {
373 this->basisCellTopology_ );
379 this->basisCellTopology_);
381 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace,Parameters::MaxNumPtsPerBasisEval>(phisAtEdgePoints , edgePts, order, OPERATOR_VALUE);
384 for (ordinal_type j=0;j<numPtsPerEdge;j++) {
386 const ordinal_type i_card = numPtsPerEdge*i+j;
389 for (ordinal_type k=0;k<cardPn;k++)
390 for (ordinal_type d=0;d<spaceDim;d++)
391 V2(i_card,k+d*cardPn) = edgeTan(d) * phisAtEdgePoints(k,j);
394 for(ordinal_type k=0; k<spaceDim; ++k) {
395 dofCoords(i_card,k) = edgePts(j,k);
396 dofCoeffs(i_card,k) = edgeTan(k);
402 tags[i_card][3] = numPtsPerEdge;
407 if(numPtsPerFace >0) {
408 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> faceTan1(
"Hcurl::Tet::In::edgeTan", spaceDim );
409 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> faceTan2(
"Hcurl::Tet::In::edgeTan", spaceDim );
411 for (ordinal_type i=0;i<numFaces;i++) {
415 this->basisCellTopology_ );
421 this->basisCellTopology_ );
423 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace,Parameters::MaxNumPtsPerBasisEval>(phisAtFacePoints , facePts, order, OPERATOR_VALUE);
426 for (ordinal_type j=0;j<numPtsPerFace;j++) {
428 const ordinal_type i_card = numEdges*numPtsPerEdge+2*numPtsPerFace*i+2*j;
429 const ordinal_type i_card_p1 = i_card+1;
432 for (ordinal_type k=0;k<cardPn;k++)
433 for (ordinal_type d=0;d<spaceDim;d++) {
434 V2(i_card,k+d*cardPn) = faceTan1(d) * phisAtFacePoints(k,j);
435 V2(i_card_p1,k+d*cardPn) = faceTan2(d) * phisAtFacePoints(k,j);
439 for(ordinal_type k=0; k<spaceDim; ++k) {
440 dofCoords(i_card,k) = facePts(j,k);
441 dofCoords(i_card_p1,k) = facePts(j,k);
442 dofCoeffs(i_card,k) = faceTan1(k);
443 dofCoeffs(i_card_p1,k) = faceTan2(k);
448 tags[i_card][2] = 2*j;
449 tags[i_card][3] = 2*numPtsPerFace;
451 tags[i_card_p1][0] = 2;
452 tags[i_card_p1][1] = i;
453 tags[i_card_p1][2] = 2*j+1;
454 tags[i_card_p1][3] = 2*numPtsPerFace;
462 if (numPtsPerCell > 0) {
463 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
464 cellPoints(
"Hcurl::Tet::In::cellPoints", numPtsPerCell , spaceDim );
466 this->basisCellTopology_ ,
471 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
472 phisAtCellPoints(
"Hcurl::Tet::In::phisAtCellPoints", cardPn , numPtsPerCell );
473 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>( phisAtCellPoints , cellPoints , order, OPERATOR_VALUE );
476 for (ordinal_type j=0;j<numPtsPerCell;j++) {
478 const ordinal_type i_card = numEdges*numPtsPerEdge+2*numFaces*numPtsPerFace+spaceDim*j;
480 for (ordinal_type k=0;k<cardPn;k++)
481 for (ordinal_type d=0;d<spaceDim;d++)
482 V2(i_card+d,d*cardPn+k) = phisAtCellPoints(k,j);
486 for(ordinal_type d=0; d<spaceDim; ++d) {
487 for(ordinal_type dim=0; dim<spaceDim; ++dim) {
488 dofCoords(i_card+d,dim) = cellPoints(j,dim);
489 dofCoeffs(i_card+d,dim) = (d==dim);
492 tags[i_card+d][0] = spaceDim;
493 tags[i_card+d][1] = 0;
494 tags[i_card+d][2] = spaceDim*j+d;
495 tags[i_card+d][3] = spaceDim*numPtsPerCell;
502 const ordinal_type lwork = card*card;
503 Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
504 vmat(
"Hcurl::Tet::In::vmat", card, card),
505 work1(
"Hcurl::Tet::In::work", lwork),
506 ipiv(
"Hcurl::Tet::In::ipiv", card);
509 for(ordinal_type i=0; i< card; ++i) {
510 for(ordinal_type j=0; j< card; ++j) {
512 for(ordinal_type k=0; k< cardVecPn; ++k)
520 lapack.GETRF(card, card,
521 vmat.data(), vmat.stride_1(),
522 (ordinal_type*)ipiv.data(),
525 INTREPID2_TEST_FOR_EXCEPTION( info != 0,
527 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRF returns nonzero info." );
530 vmat.data(), vmat.stride_1(),
531 (ordinal_type*)ipiv.data(),
535 INTREPID2_TEST_FOR_EXCEPTION( info != 0,
537 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRI returns nonzero info." );
539 for (ordinal_type i=0;i<cardVecPn;++i) {
540 for (ordinal_type j=0;j<card;++j){
542 for(ordinal_type k=0; k< card; ++k)
543 s += U(i,k)*vmat(k,j);
548 this->coeffs_ = Kokkos::create_mirror_view(
typename DT::memory_space(), coeffs);
549 Kokkos::deep_copy(this->coeffs_ , coeffs);
551 this->dofCoords_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoords);
552 Kokkos::deep_copy(this->dofCoords_, dofCoords);
554 this->dofCoeffs_ = Kokkos::create_mirror_view(
typename DT::memory_space(), dofCoeffs);
555 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
561 const ordinal_type posScDim = 0;
562 const ordinal_type posScOrd = 1;
563 const ordinal_type posDfOrd = 2;
565 OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
569 this->setOrdinalTagData(this->tagToOrdinal_,
572 this->basisCardinality_,