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
56namespace Intrepid2 {
57
58// -------------------------------------------------------------------------------------
59
60namespace Impl {
61
62template<EOperator opType>
63template<typename OutputViewType,
64typename inputViewType,
65typename workViewType,
66typename vinvViewType>
67KOKKOS_INLINE_FUNCTION
68void
69Basis_HCURL_TET_In_FEM::Serial<opType>::
70getValues( 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
139template<typename DT, ordinal_type numPtsPerEval,
140typename outputValueValueType, class ...outputValueProperties,
141typename inputPointValueType, class ...inputPointProperties,
142typename vinvValueType, class ...vinvProperties>
143void
144Basis_HCURL_TET_In_FEM::
145getValues( 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,typename DT::execution_space>::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// -------------------------------------------------------------------------------------
192template<typename DT, typename OT, typename PT>
194Basis_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_LAGRANGIAN;
202 this->basisCoordinates_ = COORDINATES_CARTESIAN;
203 this->functionSpace_ = FUNCTION_SPACE_HCURL;
204 pointType_ = pointType;
205 const ordinal_type card = this->basisCardinality_;
206
207 const ordinal_type cardPn = Intrepid2::getPnCardinality<spaceDim>(order); // dim of (P_{n}) -- smaller space
208 const ordinal_type cardPnm1 = Intrepid2::getPnCardinality<spaceDim>(order-1); // dim of (P_{n-1}) -- smaller space
209 const ordinal_type cardPnm2 = Intrepid2::getPnCardinality<spaceDim>(order-2); // dim of (P_{n-2}) -- smaller space
210 const ordinal_type cardVecPn = spaceDim*cardPn; // dim of (P_{n})^2 -- larger space
211 const ordinal_type cardVecPnm1 = spaceDim*cardPnm1; // dim of (P_{n-1})^2 -- smaller space
212 const ordinal_type cardPnm1H = cardPnm1-cardPnm2; //Homogeneous polynomial of order (n-1)
213
214 // 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.
215 // 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.)
216 INTREPID2_TEST_FOR_EXCEPTION( order > Parameters::MaxOrder, std::invalid_argument, "polynomial order exceeds the max supported by this class");
217
218 // Basis-dependent initializations
219 constexpr ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
220 constexpr ordinal_type maxCard = CardinalityHCurlTet(Parameters::MaxOrder);
221 ordinal_type tags[maxCard][tagSize];
222
223 // points are computed in the host and will be copied
224 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
225 dofCoords("Hcurl::Tet::In::dofCoords", card, spaceDim);
226
227 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
228 coeffs("Hcurl::Tet::In::coeffs", cardVecPn, card);
229
230 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
231 dofCoeffs("Hcurl::Tet::In::dofCoeffs", card, spaceDim);
232
233 // first, need to project the basis for RT space onto the
234 // orthogonal basis of degree n
235 // get coefficients of PkHx
236
237 Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace> //use LayoutLeft for Lapack
238 V1("Hcurl::Tet::In::V1", cardVecPn, cardVecPnm1 + spaceDim*cardPnm1H);
239
240
241 // these two loops get the first three sets of basis functions
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;
245
246
247 // now I need to integrate { (x,y) \times phi } against the big basis
248 // first, get a cubature rule.
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() );
252 myCub.getCubature( cubPoints , cubWeights );
253
254 // tabulate the scalar orthonormal basis at cubature points
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);
257
258 // 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)
259 for (ordinal_type i=0;i<cardPn;i++) {
260 for (ordinal_type j=0;j<cardPnm1H;j++) { // loop over homogeneous polynomials
261 for (ordinal_type d=0; d< spaceDim; ++d) {
262 scalarType integral(0);
263 for (ordinal_type k=0;k<myCub.getNumPoints();k++)
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;
270 }
271 }
272 }
273
274
275
276
277
278 // now I need to set up an SVD to get a basis for the space
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);
285
286
287
288 ordinal_type info = 0;
289 Teuchos::LAPACK<ordinal_type,scalarType> lapack;
290
291
292 lapack.GESVD( 'A',
293 'N',
294 V1.extent(0) ,
295 V1.extent(1) ,
296 V1.data() ,
297 V1.stride_1() ,
298 S.data() ,
299 U.data() ,
300 U.stride_1() ,
301 Vt.data() ,
302 Vt.stride_1() ,
303 work.data() ,
304 5*cardVecPn ,
305 rWork.data() ,
306 &info );
307
308
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());
313
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");
316#endif
317
318 // next, apply the RT nodes (rows) to the basis for (P_n)^2 (columns)
319 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
320 V2("Hcurl::Tet::In::V2", card ,cardVecPn);
321
322 const ordinal_type numEdges = this->basisCellTopology_.getEdgeCount();
323 const ordinal_type numFaces = this->basisCellTopology_.getFaceCount();
324
325 // first numEdges * degree nodes are normals at each edge
326 // get the points on the line
327
328 shards::CellTopology edgeTop(shards::getCellTopologyData<shards::Line<2> >() );
329 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
330
331 const int numPtsPerEdge = PointTools::getLatticeSize( edgeTop ,
332 order+1 ,
333 1 );
334
335 const int numPtsPerFace = PointTools::getLatticeSize( faceTop ,
336 order+1 ,
337 1 );
338
339 const int numPtsPerCell = PointTools::getLatticeSize( this->basisCellTopology_ ,
340 order+1 ,
341 1 );
342
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 );
345
346 // construct lattice
347 const ordinal_type offset = 1;
348
349
350
351 PointTools::getLattice( linePts,
352 edgeTop,
353 order+1, offset,
354 pointType );
355
357 faceTop,
358 order+1, offset,
359 pointType );
360
361 // holds the image of the line points
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);
366
367 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> edgeTan("Hcurl::Tet::In::edgeTan", spaceDim );
368
369 // these are tangents scaled by the appropriate edge lengths.
370 for (ordinal_type i=0;i<numEdges;i++) { // loop over edges
372 i ,
373 this->basisCellTopology_ );
374
376 linePts ,
377 1 ,
378 i ,
379 this->basisCellTopology_);
380
381 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace,Parameters::MaxNumPtsPerBasisEval>(phisAtEdgePoints , edgePts, order, OPERATOR_VALUE);
382
383 // loop over points (rows of V2)
384 for (ordinal_type j=0;j<numPtsPerEdge;j++) {
385
386 const ordinal_type i_card = numPtsPerEdge*i+j;
387
388 // loop over orthonormal basis functions (columns of V2)
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);
392
393 //save dof coordinates and coefficients
394 for(ordinal_type k=0; k<spaceDim; ++k) {
395 dofCoords(i_card,k) = edgePts(j,k);
396 dofCoeffs(i_card,k) = edgeTan(k);
397 }
398
399 tags[i_card][0] = 1; // edge dof
400 tags[i_card][1] = i; // edge id
401 tags[i_card][2] = j; // local dof id
402 tags[i_card][3] = numPtsPerEdge; // total vert dof
403
404 }
405 }
406
407 if(numPtsPerFace >0) {//handle faces if needed (order >1)
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 );
410
411 for (ordinal_type i=0;i<numFaces;i++) { // loop over faces
413 faceTan2,
414 i ,
415 this->basisCellTopology_ );
416
418 triPts ,
419 2 ,
420 i ,
421 this->basisCellTopology_ );
422
423 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace,Parameters::MaxNumPtsPerBasisEval>(phisAtFacePoints , facePts, order, OPERATOR_VALUE);
424
425 // loop over points (rows of V2)
426 for (ordinal_type j=0;j<numPtsPerFace;j++) {
427
428 const ordinal_type i_card = numEdges*numPtsPerEdge+2*numPtsPerFace*i+2*j;
429 const ordinal_type i_card_p1 = i_card+1; // creating a temp otherwise nvcc gets confused
430
431 // loop over orthonormal basis functions (columns of V2)
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);
436 }
437
438 //save dof coordinates
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);
444 }
445
446 tags[i_card][0] = 2; // face dof
447 tags[i_card][1] = i; // face id
448 tags[i_card][2] = 2*j; // local face id
449 tags[i_card][3] = 2*numPtsPerFace; // total face dof
450
451 tags[i_card_p1][0] = 2; // face dof
452 tags[i_card_p1][1] = i; // face id
453 tags[i_card_p1][2] = 2*j+1; // local face id
454 tags[i_card_p1][3] = 2*numPtsPerFace; // total face dof
455
456 }
457 }
458 }
459
460
461 // internal dof, if needed
462 if (numPtsPerCell > 0) {
463 Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace>
464 cellPoints( "Hcurl::Tet::In::cellPoints", numPtsPerCell , spaceDim );
465 PointTools::getLattice( cellPoints ,
466 this->basisCellTopology_ ,
467 order + 1 ,
468 1 ,
469 pointType );
470
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 );
474
475 // copy values into right positions of V2
476 for (ordinal_type j=0;j<numPtsPerCell;j++) {
477
478 const ordinal_type i_card = numEdges*numPtsPerEdge+2*numFaces*numPtsPerFace+spaceDim*j;
479
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);
483
484
485 //save dof coordinates
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);
490 }
491
492 tags[i_card+d][0] = spaceDim; // elem dof
493 tags[i_card+d][1] = 0; // elem id
494 tags[i_card+d][2] = spaceDim*j+d; // local dof id
495 tags[i_card+d][3] = spaceDim*numPtsPerCell; // total vert dof
496 }
497 }
498 }
499
500 // form Vandermonde matrix. Actually, this is the transpose of the VDM,
501 // so we transpose on copy below.
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);
507
508 //vmat = V2*U;
509 for(ordinal_type i=0; i< card; ++i) {
510 for(ordinal_type j=0; j< card; ++j) {
511 scalarType s=0;
512 for(ordinal_type k=0; k< cardVecPn; ++k)
513 s += V2(i,k)*U(k,j);
514 vmat(i,j) = s;
515 }
516 }
517
518 info = 0;
519
520 lapack.GETRF(card, card,
521 vmat.data(), vmat.stride_1(),
522 (ordinal_type*)ipiv.data(),
523 &info);
524
525 INTREPID2_TEST_FOR_EXCEPTION( info != 0,
526 std::runtime_error ,
527 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRF returns nonzero info." );
528
529 lapack.GETRI(card,
530 vmat.data(), vmat.stride_1(),
531 (ordinal_type*)ipiv.data(),
532 work1.data(), lwork,
533 &info);
534
535 INTREPID2_TEST_FOR_EXCEPTION( info != 0,
536 std::runtime_error ,
537 ">>> ERROR: (Intrepid2::Basis_HCURL_TET_In_FEM) lapack.GETRI returns nonzero info." );
538
539 for (ordinal_type i=0;i<cardVecPn;++i) {
540 for (ordinal_type j=0;j<card;++j){
541 scalarType s=0;
542 for(ordinal_type k=0; k< card; ++k)
543 s += U(i,k)*vmat(k,j);
544 coeffs(i,j) = s;
545 }
546 }
547
548 this->coeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), coeffs);
549 Kokkos::deep_copy(this->coeffs_ , coeffs);
550
551 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
552 Kokkos::deep_copy(this->dofCoords_, dofCoords);
553
554 this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffs);
555 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
556
557
558 // set tags
559 {
560 // Basis-dependent initializations
561 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
562 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
563 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
564
565 OrdinalTypeArray1DHost tagView(&tags[0][0], card*tagSize);
566
567 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
568 // tags are constructed on host
569 this->setOrdinalTagData(this->tagToOrdinal_,
570 this->ordinalToTag_,
571 tagView,
572 this->basisCardinality_,
573 tagSize,
574 posScDim,
575 posScOrd,
576 posDfOrd);
577 }
578}
579} // namespace Intrepid2
580#endif
Header file for the Intrepid2::CubatureDirectTetDefault class.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
Basis_HCURL_TET_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
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 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.
static void getReferenceFaceTangents(Kokkos::DynRankView< refFaceTanValueType, refFaceTanProperties... > refFaceTanU, Kokkos::DynRankView< refFaceTanValueType, refFaceTanProperties... > refFaceTanV, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
Defines direct integration rules on a tetrahedron.
virtual void getCubature(PointViewType cubPoints, weightViewType cubWeights) const override
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
virtual ordinal_type getNumPoints() const override
Returns the number of cubature points.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
static ordinal_type getLatticeSize(const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...
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...