48 #ifndef __INTREPID2_HVOL_TET_CN_FEM_DEF_HPP__
49 #define __INTREPID2_HVOL_TET_CN_FEM_DEF_HPP__
59 template<EOperator opType>
60 template<
typename outputViewType,
61 typename inputViewType,
62 typename workViewType,
63 typename vinvViewType>
64 KOKKOS_INLINE_FUNCTION
66 Basis_HVOL_TET_Cn_FEM::Serial<opType>::
67 getValues( outputViewType output,
68 const inputViewType input,
70 const vinvViewType vinv ) {
72 constexpr ordinal_type spaceDim = 3;
74 card = vinv.extent(0),
75 npts = input.extent(0);
78 ordinal_type order = 0;
80 if (card == Intrepid2::getPnCardinality<spaceDim>(p)) {
86 typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
87 auto vcprop = Kokkos::common_view_alloc_prop(work);
88 auto ptr = work.data();
91 case OPERATOR_VALUE: {
92 const viewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts);
93 workViewType dummyView;
95 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
96 Serial<opType>::getValues(phis, input, dummyView, order);
98 for (ordinal_type i=0;i<card;++i)
99 for (ordinal_type j=0;j<npts;++j) {
100 output.access(i,j) = 0.0;
101 for (ordinal_type k=0;k<card;++k)
102 output.access(i,j) += vinv(k,i)*phis.access(k,j);
108 const viewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts, spaceDim);
109 ptr += card*npts*spaceDim*get_dimension_scalar(work);
110 const viewType workView(Kokkos::view_wrap(ptr, vcprop), card, npts, spaceDim+1);
111 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
112 Serial<opType>::getValues(phis, input, workView, order);
114 for (ordinal_type i=0;i<card;++i)
115 for (ordinal_type j=0;j<npts;++j)
116 for (ordinal_type k=0;k<spaceDim;++k) {
117 output.access(i,j,k) = 0.0;
118 for (ordinal_type l=0;l<card;++l)
119 output.access(i,j,k) += vinv(l,i)*phis.access(l,j,k);
132 const ordinal_type dkcard = getDkCardinality<opType,spaceDim>();
133 const viewType phis(Kokkos::view_wrap(ptr, vcprop), card, npts, dkcard);
134 workViewType dummyView;
136 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::
137 Serial<opType>::getValues(phis, input, dummyView, order);
139 for (ordinal_type i=0;i<card;++i)
140 for (ordinal_type j=0;j<npts;++j)
141 for (ordinal_type k=0;k<dkcard;++k) {
142 output.access(i,j,k) = 0.0;
143 for (ordinal_type l=0;l<card;++l)
144 output.access(i,j,k) += vinv(l,i)*phis.access(l,j,k);
149 INTREPID2_TEST_FOR_ABORT(
true,
150 ">>> ERROR (Basis_HVOL_TET_Cn_FEM): Operator type not implemented");
155 template<
typename SpT, ordinal_type numPtsPerEval,
156 typename outputValueValueType,
class ...outputValueProperties,
157 typename inputPointValueType,
class ...inputPointProperties,
158 typename vinvValueType,
class ...vinvProperties>
160 Basis_HVOL_TET_Cn_FEM::
161 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
162 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
163 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
164 const EOperator operatorType) {
165 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
166 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
167 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
168 typedef typename ExecSpace<typename inputPointViewType::execution_space,SpT>::ExecSpaceType ExecSpaceType;
171 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
172 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
173 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
174 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
176 typedef typename inputPointViewType::value_type inputPointType;
178 const ordinal_type cardinality = outputValues.extent(0);
179 const ordinal_type spaceDim = 3;
181 ordinal_type order = 0;
182 while((Intrepid2::getPnCardinality<spaceDim>(++order) != cardinality) && (order !=
Parameters::MaxOrder));
184 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
185 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
187 switch (operatorType) {
188 case OPERATOR_VALUE: {
189 auto bufferSize = Basis_HVOL_TET_Cn_FEM::Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
190 workViewType work(Kokkos::view_alloc(
"Basis_HVOL_TET_Cn_FEM::getValues::work", vcprop), bufferSize, inputPoints.extent(0));
191 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
192 OPERATOR_VALUE,numPtsPerEval> FunctorType;
193 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
198 auto bufferSize = Basis_HVOL_TET_Cn_FEM::Serial<OPERATOR_D1>::getWorkSizePerPoint(order);
199 workViewType work(Kokkos::view_alloc(
"Basis_HVOL_TET_Cn_FEM::getValues::work", vcprop), bufferSize, inputPoints.extent(0));
200 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
201 OPERATOR_D1,numPtsPerEval> FunctorType;
202 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
206 auto bufferSize = Basis_HVOL_TET_Cn_FEM::Serial<OPERATOR_D2>::getWorkSizePerPoint(order);
207 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
208 OPERATOR_D2,numPtsPerEval> FunctorType;
209 workViewType work(Kokkos::view_alloc(
"Basis_HVOL_TET_Cn_FEM::getValues::work", vcprop), bufferSize, inputPoints.extent(0));
210 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
221 INTREPID2_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
222 ">>> ERROR (Basis_HVOL_TET_Cn_FEM): Operator type not implemented" );
229 template<
typename SpT,
typename OT,
typename PT>
232 const EPointType pointType ) {
233 constexpr ordinal_type spaceDim = 3;
235 this->basisCardinality_ = Intrepid2::getPnCardinality<spaceDim>(order);
236 this->basisDegree_ = order;
237 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
238 this->basisType_ = BASIS_FEM_FIAT;
239 this->basisCoordinates_ = COORDINATES_CARTESIAN;
241 const ordinal_type card = this->basisCardinality_;
244 Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
245 dofCoords(
"HVOL::Tet::Cn::dofCoords", card, spaceDim);
248 const ordinal_type offset = 1;
250 this->basisCellTopology_,
251 order+spaceDim+offset, offset,
254 this->dofCoords_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), dofCoords);
255 Kokkos::deep_copy(this->dofCoords_, dofCoords);
259 const ordinal_type lwork = card*card;
260 Kokkos::DynRankView<scalarType,Kokkos::LayoutLeft,Kokkos::HostSpace>
261 vmat(
"HVOL::Tet::Cn::vmat", card, card),
262 work(
"HVOL::Tet::Cn::work", lwork),
263 ipiv(
"HVOL::Tet::Cn::ipiv", card);
265 Impl::Basis_HGRAD_TET_Cn_FEM_ORTH::getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>(vmat, dofCoords, order, OPERATOR_VALUE);
267 ordinal_type info = 0;
268 Teuchos::LAPACK<ordinal_type,scalarType> lapack;
270 lapack.GETRF(card, card,
271 vmat.data(), vmat.stride_1(),
272 (ordinal_type*)ipiv.data(),
275 INTREPID2_TEST_FOR_EXCEPTION( info != 0,
277 ">>> ERROR: (Intrepid2::Basis_HVOL_TET_Cn_FEM) lapack.GETRF returns nonzero info." );
280 vmat.data(), vmat.stride_1(),
281 (ordinal_type*)ipiv.data(),
285 INTREPID2_TEST_FOR_EXCEPTION( info != 0,
287 ">>> ERROR: (Intrepid2::Basis_HVOL_TET_Cn_FEM) lapack.GETRI returns nonzero info." );
290 Kokkos::DynRankView<scalarType,typename SpT::array_layout,Kokkos::HostSpace>
291 vinv(
"HVOL::Line::Cn::vinv", card, card);
293 for (ordinal_type i=0;i<card;++i)
294 for (ordinal_type j=0;j<card;++j)
295 vinv(i,j) = vmat(j,i);
297 this->vinv_ = Kokkos::create_mirror_view(
typename SpT::memory_space(), vinv);
298 Kokkos::deep_copy(this->vinv_ , vinv);
303 constexpr ordinal_type tagSize = 4;
304 const ordinal_type posScDim = 0;
305 const ordinal_type posScOrd = 1;
306 const ordinal_type posDfOrd = 2;
308 constexpr ordinal_type maxCard = Intrepid2::getPnCardinality<spaceDim, Parameters::MaxOrder>();
309 ordinal_type tags[maxCard][tagSize];
312 numElemDof = this->basisCardinality_;
315 ordinal_type elemId = 0;
316 for (ordinal_type i=0;i<this->basisCardinality_;++i) {
318 tags[i][0] = spaceDim;
320 tags[i][2] = elemId++;
321 tags[i][3] = numElemDof;
328 this->setOrdinalTagData(this->tagToOrdinal_,
331 this->basisCardinality_,
Basis_HVOL_TET_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.