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>
277 Basis_HGRAD_QUAD_Cn_FEM<DT,OT,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 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
308 constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
309 ordinal_type tags[maxCardLine*maxCardLine][4];
310
311 const ordinal_type vert[2][2] = { {0,1}, {3,2} }; //[y][x]
312
313 const ordinal_type edge_x[2] = {0,2};
314 const ordinal_type edge_y[2] = {3,1};
315 {
316 ordinal_type idx = 0;
317 for (ordinal_type j=0;j<cardLine;++j) { // y
318 const auto tag_y = lineBasis.getDofTag(j);
319 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
320 const auto tag_x = lineBasis.getDofTag(i);
321
322 if (tag_x(0) == 0 && tag_y(0) == 0) {
323 // vertices
324 tags[idx][0] = 0; // vertex dof
325 tags[idx][1] = vert[tag_y(1)][tag_x(1)]; // vertex id
326 tags[idx][2] = 0; // local dof id
327 tags[idx][3] = 1; // total number of dofs in this vertex
328 } else if (tag_x(0) == 1 && tag_y(0) == 0) {
329 // edge: x edge, y vert
330 tags[idx][0] = 1; // edge dof
331 tags[idx][1] = edge_x[tag_y(1)];
332 tags[idx][2] = tag_x(2); // local dof id
333 tags[idx][3] = tag_x(3); // total number of dofs in this vertex
334 } else if (tag_x(0) == 0 && tag_y(0) == 1) {
335 // edge: x vert, y edge
336 tags[idx][0] = 1; // edge dof
337 tags[idx][1] = edge_y[tag_x(1)];
338 tags[idx][2] = tag_y(2); // local dof id
339 tags[idx][3] = tag_y(3); // total number of dofs in this vertex
340 } else {
341 // interior
342 tags[idx][0] = 2; // interior dof
343 tags[idx][1] = 0;
344 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
345 tags[idx][3] = tag_x(3)*tag_y(3); // total number of dofs in this vertex
346 }
347 }
348 }
349 }
350
351 OrdinalTypeArray1DHost tagView(&tags[0][0], this->basisCardinality_*4);
352
353 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
354 // tags are constructed on host
355 this->setOrdinalTagData(this->tagToOrdinal_,
356 this->ordinalToTag_,
357 tagView,
358 this->basisCardinality_,
359 tagSize,
360 posScDim,
361 posScOrd,
362 posDfOrd);
363 }
364
365 // dofCoords on host and create its mirror view to device
366 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
367 dofCoordsHost("dofCoordsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
368
369 Kokkos::DynRankView<typename ScalarViewType::value_type,DT>
370 dofCoordsLine("dofCoordsLine", cardLine, 1);
371
372 lineBasis.getDofCoords(dofCoordsLine);
373 auto dofCoordsLineHost = Kokkos::create_mirror_view(dofCoordsLine);
374 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
375 {
376 ordinal_type idx = 0;
377 for (ordinal_type j=0;j<cardLine;++j) { // y
378 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
379 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
380 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
381 }
382 }
383 }
384
385 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
386 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
387 }
388
389}// namespace Intrepid2
390
391#endif
Intrepid2::getOperatorOrder
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Definition: Intrepid2_BasisDef.hpp:194
Intrepid2::Basis_HGRAD_LINE_Cn_FEM
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
Definition: Intrepid2_HGRAD_LINE_Cn_FEM.hpp:176
Intrepid2::Basis_HGRAD_LINE_Cn_FEM::getDofCoords
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
Definition: Intrepid2_HGRAD_LINE_Cn_FEM.hpp:226
Intrepid2::Basis_HGRAD_QUAD_Cn_FEM::Basis_HGRAD_QUAD_Cn_FEM
Basis_HGRAD_QUAD_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition: Intrepid2_HGRAD_QUAD_Cn_FEMDef.hpp:278
Intrepid2::Basis::getDofTag
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
Definition: Intrepid2_Basis.hpp:785
Intrepid2::Basis::getCardinality
ordinal_type getCardinality() const
Returns cardinality of the basis.
Definition: Intrepid2_Basis.hpp:656
Intrepid2::Parameters::MaxOrder
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.
Definition: Intrepid2_Types.hpp:123
Generated by doxygen 1.9.4