Intrepid2
Intrepid2_HGRAD_QUAD_Cn_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_HGRAD_QUAD_CN_FEM_DEF_HPP__
50#define __INTREPID2_HGRAD_QUAD_CN_FEM_DEF_HPP__
51
52namespace Intrepid2 {
53
54 // -------------------------------------------------------------------------------------
55 namespace Impl {
56
57 template<EOperator opType>
58 template<typename OutputViewType,
59 typename inputViewType,
60 typename workViewType,
61 typename vinvViewType>
62 KOKKOS_INLINE_FUNCTION
63 void
64 Basis_HGRAD_QUAD_Cn_FEM::Serial<opType>::
65 getValues( OutputViewType output,
66 const inputViewType input,
67 workViewType work,
68 const vinvViewType vinv,
69 const ordinal_type operatorDn ) {
70 ordinal_type opDn = operatorDn;
71
72 const ordinal_type cardLine = vinv.extent(0);
73 const ordinal_type npts = input.extent(0);
74
75 typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
76 const auto input_x = Kokkos::subview(input, Kokkos::ALL(), range_type(0,1));
77 const auto input_y = Kokkos::subview(input, Kokkos::ALL(), range_type(1,2));
78
79 const int dim_s = get_dimension_scalar(work);
80 auto ptr0 = work.data();
81 auto ptr1 = work.data()+cardLine*npts*dim_s;
82 auto ptr2 = work.data()+2*cardLine*npts*dim_s;
83
84 typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
85 auto vcprop = Kokkos::common_view_alloc_prop(work);
86
87 switch (opType) {
88 case OPERATOR_VALUE: {
89 viewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
90 viewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
91 viewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
92
93 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
94 getValues(output_x, input_x, work_line, vinv);
95
96 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
97 getValues(output_y, input_y, work_line, vinv);
98
99 // tensor product
100 ordinal_type idx = 0;
101 for (ordinal_type j=0;j<cardLine;++j) // y
102 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
103 for (ordinal_type k=0;k<npts;++k)
104 output.access(idx,k) = output_x.access(i,k)*output_y.access(j,k);
105 break;
106 }
107 case OPERATOR_CURL: {
108 for (auto l=0;l<2;++l) {
109 viewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
110
111 viewType output_x, output_y;
112
113 typename workViewType::value_type s = 0.0;
114 if (l) {
115 // l = 1
116 output_x = viewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
117 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
118 getValues(output_x, input_x, work_line, vinv, 1);
119
120 output_y = viewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
121 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
122 getValues(output_y, input_y, work_line, vinv);
123
124 s = -1.0;
125 } else {
126 // l = 0
127 output_x = viewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
128 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
129 getValues(output_x, input_x, work_line, vinv);
130
131 output_y = viewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts, 1);
132 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
133 getValues(output_y, input_y, work_line, vinv, 1);
134
135 s = 1.0;
136 }
137
138 // tensor product (extra dimension of ouput x and y are ignored)
139 ordinal_type idx = 0;
140 for (ordinal_type j=0;j<cardLine;++j) // y
141 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
142 for (ordinal_type k=0;k<npts;++k)
143 output.access(idx,k,l) = s*output_x.access(i,k,0)*output_y.access(j,k,0);
144 }
145 break;
146 }
147 case OPERATOR_GRAD:
148 case OPERATOR_D1:
149 case OPERATOR_D2:
150 case OPERATOR_D3:
151 case OPERATOR_D4:
152 case OPERATOR_D5:
153 case OPERATOR_D6:
154 case OPERATOR_D7:
155 case OPERATOR_D8:
156 case OPERATOR_D9:
157 case OPERATOR_D10:
158 opDn = getOperatorOrder(opType);
159 case OPERATOR_Dn: {
160 const auto dkcard = opDn + 1;
161 for (auto l=0;l<dkcard;++l) {
162 viewType work_line(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
163
164 viewType output_x, output_y;
165
166 const auto mult_x = opDn - l;
167 const auto mult_y = l;
168
169 if (mult_x) {
170 output_x = viewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts, 1);
171 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
172 getValues(output_x, input_x, work_line, vinv, mult_x);
173 } else {
174 output_x = viewType(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
175 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
176 getValues(output_x, input_x, work_line, vinv);
177 }
178
179 if (mult_y) {
180 output_y = viewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts, 1);
181 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
182 getValues(output_y, input_y, work_line, vinv, mult_y);
183 } else {
184 output_y = viewType(Kokkos::view_wrap(ptr2, vcprop), cardLine, npts);
185 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
186 getValues(output_y, input_y, work_line, vinv);
187 }
188
189 // tensor product (extra dimension of ouput x and y are ignored)
190 ordinal_type idx = 0;
191 for (ordinal_type j=0;j<cardLine;++j) // y
192 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
193 for (ordinal_type k=0;k<npts;++k)
194 output.access(idx,k,l) = output_x.access(i,k,0)*output_y.access(j,k,0);
195 }
196 break;
197 }
198 default: {
199 INTREPID2_TEST_FOR_ABORT( true,
200 ">>> ERROR: (Intrepid2::Basis_HGRAD_QUAD_Cn_FEM::Serial::getValues) operator is not supported" );
201 }
202 }
203 }
204
205 template<typename DT, ordinal_type numPtsPerEval,
206 typename outputValueValueType, class ...outputValueProperties,
207 typename inputPointValueType, class ...inputPointProperties,
208 typename vinvValueType, class ...vinvProperties>
209 void
210 Basis_HGRAD_QUAD_Cn_FEM::
211 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
212 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
213 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
214 const EOperator operatorType ) {
215 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
216 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
217 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
218 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
219
220 // loopSize corresponds to cardinality
221 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
222 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
223 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
224 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
225
226 typedef typename inputPointViewType::value_type inputPointType;
227
228 const ordinal_type cardinality = outputValues.extent(0);
229 const ordinal_type cardLine = std::sqrt(cardinality);
230 const ordinal_type workSize = 3*cardLine;
231
232 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
233 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
234 workViewType work(Kokkos::view_alloc("Basis_HGRAD_QUAD_Cn_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
235
236 switch (operatorType) {
237 case OPERATOR_VALUE: {
238 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
239 OPERATOR_VALUE,numPtsPerEval> FunctorType;
240 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
241 break;
242 }
243 case OPERATOR_CURL: {
244 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
245 OPERATOR_CURL,numPtsPerEval> FunctorType;
246 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
247 break;
248 }
249 case OPERATOR_GRAD:
250 case OPERATOR_D1:
251 case OPERATOR_D2:
252 case OPERATOR_D3:
253 case OPERATOR_D4:
254 case OPERATOR_D5:
255 case OPERATOR_D6:
256 case OPERATOR_D7:
257 case OPERATOR_D8:
258 case OPERATOR_D9:
259 case OPERATOR_D10: {
260 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
261 OPERATOR_Dn,numPtsPerEval> FunctorType;
262 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work,
263 getOperatorOrder(operatorType)) );
264 break;
265 }
266 default: {
267 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
268 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented" );
269 // break;commented out because exception
270 }
271 }
272 }
273 }
274
275 // -------------------------------------------------------------------------------------
276 template<typename DT, typename OT, typename PT>
278 Basis_HGRAD_QUAD_Cn_FEM( const ordinal_type order,
279 const EPointType pointType ) {
280 // INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
281 // pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
282 // ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): pointType must be either equispaced or warpblend." );
283
284 // this should be in host
285 Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
286 const auto cardLine = lineBasis.getCardinality();
287
288 this->vinv_ = Kokkos::DynRankView<typename ScalarViewType::value_type,DT>("Hgrad::Quad::Cn::vinv", cardLine, cardLine);
289 lineBasis.getVandermondeInverse(this->vinv_);
290
291 this->basisCardinality_ = cardLine*cardLine;
292 this->basisDegree_ = order;
293 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
294 this->basisType_ = BASIS_FEM_LAGRANGIAN;
295 this->basisCoordinates_ = COORDINATES_CARTESIAN;
296 this->functionSpace_ = FUNCTION_SPACE_HGRAD;
297 pointType_ = pointType;
298
299 // initialize tags
300 {
301 // Basis-dependent initializations
302 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
303 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
304 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
305 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
306
307 // 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.
308 // 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.)
309 INTREPID2_TEST_FOR_EXCEPTION( order > Parameters::MaxOrder, std::invalid_argument, "polynomial order exceeds the max supported by this class");
310 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
311 constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
312 ordinal_type tags[maxCardLine*maxCardLine][4];
313
314 const ordinal_type vert[2][2] = { {0,1}, {3,2} }; //[y][x]
315
316 const ordinal_type edge_x[2] = {0,2};
317 const ordinal_type edge_y[2] = {3,1};
318 {
319 ordinal_type idx = 0;
320 for (ordinal_type j=0;j<cardLine;++j) { // y
321 const auto tag_y = lineBasis.getDofTag(j);
322 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
323 const auto tag_x = lineBasis.getDofTag(i);
324
325 if (tag_x(0) == 0 && tag_y(0) == 0) {
326 // vertices
327 tags[idx][0] = 0; // vertex dof
328 tags[idx][1] = vert[tag_y(1)][tag_x(1)]; // vertex id
329 tags[idx][2] = 0; // local dof id
330 tags[idx][3] = 1; // total number of dofs in this vertex
331 } else if (tag_x(0) == 1 && tag_y(0) == 0) {
332 // edge: x edge, y vert
333 tags[idx][0] = 1; // edge dof
334 tags[idx][1] = edge_x[tag_y(1)];
335 tags[idx][2] = tag_x(2); // local dof id
336 tags[idx][3] = tag_x(3); // total number of dofs in this vertex
337 } else if (tag_x(0) == 0 && tag_y(0) == 1) {
338 // edge: x vert, y edge
339 tags[idx][0] = 1; // edge dof
340 tags[idx][1] = edge_y[tag_x(1)];
341 tags[idx][2] = tag_y(2); // local dof id
342 tags[idx][3] = tag_y(3); // total number of dofs in this vertex
343 } else {
344 // interior
345 tags[idx][0] = 2; // interior dof
346 tags[idx][1] = 0;
347 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
348 tags[idx][3] = tag_x(3)*tag_y(3); // total number of dofs in this vertex
349 }
350 }
351 }
352 }
353
354 OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
355
356 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
357 // tags are constructed on host
358 this->setOrdinalTagData(this->tagToOrdinal_,
359 this->ordinalToTag_,
360 tagView,
361 this->basisCardinality_,
362 tagSize,
363 posScDim,
364 posScOrd,
365 posDfOrd);
366 }
367
368 // dofCoords on host and create its mirror view to device
369 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
370 dofCoordsHost("dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
371
372 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
373 dofCoordsLine("dofCoordsLine", cardLine, 1);
374
375 lineBasis.getDofCoords(dofCoordsLine);
376 auto dofCoordsLineHost = Kokkos::create_mirror_view(dofCoordsLine);
377 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
378 {
379 ordinal_type idx = 0;
380 for (ordinal_type j=0;j<cardLine;++j) { // y
381 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
382 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
383 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
384 }
385 }
386 }
387
388 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
389 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
390 }
391
392}// namespace Intrepid2
393
394#endif
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
Basis_HGRAD_QUAD_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
ordinal_type getCardinality() const
Returns cardinality of the basis.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.