Intrepid2
Intrepid2_HCURL_QUAD_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_QUAD_IN_FEM_DEF_HPP__
50#define __INTREPID2_HCURL_QUAD_IN_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_HCURL_QUAD_In_FEM::Serial<opType>::
65 getValues( OutputViewType output,
66 const inputViewType input,
67 workViewType work,
68 const vinvViewType vinvLine,
69 const vinvViewType vinvBubble) {
70 const ordinal_type cardLine = vinvLine.extent(0);
71 const ordinal_type cardBubble = vinvBubble.extent(0);
72
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 workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
90 viewType outputLine(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts);
91 viewType outputBubble(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
92
93 // tensor product
94 ordinal_type idx = 0;
95 {
96 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
97 getValues(outputBubble, input_x, workLine, vinvBubble);
98
99 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
100 getValues(outputLine, input_y, workLine, vinvLine);
101
102 // x component (lineBasis(y) bubbleBasis(x))
103 const auto output_x = outputBubble;
104 const auto output_y = outputLine;
105
106 for (ordinal_type j=0;j<cardLine;++j) // y
107 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
108 for (ordinal_type k=0;k<npts;++k) {
109 output.access(idx,k,0) = output_x.access(i,k)*output_y.access(j,k);
110 output.access(idx,k,1) = 0.0;
111 }
112 }
113
114 {
115 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
116 getValues(outputBubble, input_y, workLine, vinvBubble);
117
118 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
119 getValues(outputLine, input_x, workLine, vinvLine);
120
121 // y component (bubbleBasis(y) lineBasis(x))
122 const auto output_x = outputLine;
123 const auto output_y = outputBubble;
124 for (ordinal_type j=0;j<cardBubble;++j) // y
125 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
126 for (ordinal_type k=0;k<npts;++k) {
127 output.access(idx,k,0) = 0.0;
128 output.access(idx,k,1) = output_x.access(i,k)*output_y.access(j,k);
129 }
130 }
131
132 break;
133 }
134 case OPERATOR_CURL: {
135 ordinal_type idx = 0;
136 { // x - component
137 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
138 // x bubble value
139 viewType output_x(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
140 // y line grad
141 viewType output_y(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
142
143 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
144 getValues(output_x, input_x, workLine, vinvBubble);
145
146 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
147 getValues(output_y, input_y, workLine, vinvLine, 1);
148
149 // tensor product (extra dimension of ouput x and y are ignored)
150 for (ordinal_type j=0;j<cardLine;++j) // y
151 for (ordinal_type i=0;i<cardBubble;++i,++idx) // x
152 for (ordinal_type k=0;k<npts;++k)
153 output.access(idx,k) = -output_x.access(i,k)*output_y.access(j,k,0);
154 }
155 { // y - component
156 viewType workLine(Kokkos::view_wrap(ptr0, vcprop), cardLine, npts);
157 // x line grad
158 viewType output_x(Kokkos::view_wrap(ptr1, vcprop), cardLine, npts,1);
159 // y bubble value
160 viewType output_y(Kokkos::view_wrap(ptr2, vcprop), cardBubble, npts);
161
162 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_VALUE>::
163 getValues(output_y, input_y, workLine, vinvBubble);
164
165 Impl::Basis_HGRAD_LINE_Cn_FEM::Serial<OPERATOR_Dn>::
166 getValues(output_x, input_x, workLine, vinvLine, 1);
167
168 // tensor product (extra dimension of ouput x and y are ignored)
169 for (ordinal_type j=0;j<cardBubble;++j) // y
170 for (ordinal_type i=0;i<cardLine;++i,++idx) // x
171 for (ordinal_type k=0;k<npts;++k)
172 output.access(idx,k) = output_x.access(i,k,0)*output_y.access(j,k);
173 }
174 break;
175 }
176 default: {
177 INTREPID2_TEST_FOR_ABORT( true,
178 ">>> ERROR: (Intrepid2::Basis_HCURL_QUAD_In_FEM::Serial::getValues) operator is not supported" );
179 }
180 }
181 }
182
183 template<typename DT, ordinal_type numPtsPerEval,
184 typename outputValueValueType, class ...outputValueProperties,
185 typename inputPointValueType, class ...inputPointProperties,
186 typename vinvValueType, class ...vinvProperties>
187 void
188 Basis_HCURL_QUAD_In_FEM::
189 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
190 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
191 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvLine,
192 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvBubble,
193 const EOperator operatorType ) {
194 typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
195 typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
196 typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
197 typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
198
199 // loopSize corresponds to cardinality
200 const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
201 const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
202 const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
203 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
204
205 typedef typename inputPointViewType::value_type inputPointType;
206
207 const ordinal_type cardinality = outputValues.extent(0);
208 //get basis order based on basis cardinality.
209 ordinal_type order = 0;
210 ordinal_type cardBubble; // = std::sqrt(cardinality/2);
211 ordinal_type cardLine; // = cardBubble+1;
212 do {
213 cardBubble = Intrepid2::getPnCardinality<1>(order);
214 cardLine = Intrepid2::getPnCardinality<1>(++order);
215 } while((2*cardBubble*cardLine != cardinality) && (order != Parameters::MaxOrder));
216
217 auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
218 typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
219
220 switch (operatorType) {
221 case OPERATOR_VALUE: {
222 auto workSize = Serial<OPERATOR_VALUE>::getWorkSizePerPoint(order);
223 workViewType work(Kokkos::view_alloc("Basis_HCURL_QUAD_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
224 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
225 OPERATOR_VALUE,numPtsPerEval> FunctorType;
226 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
227 break;
228 }
229 case OPERATOR_CURL: {
230 auto workSize = Serial<OPERATOR_CURL>::getWorkSizePerPoint(order);
231 workViewType work(Kokkos::view_alloc("Basis_HCURL_QUAD_In_FEM::getValues::work", vcprop), workSize, inputPoints.extent(0));
232 typedef Functor<outputValueViewType,inputPointViewType,vinvViewType, workViewType,
233 OPERATOR_CURL,numPtsPerEval> FunctorType;
234 Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinvLine, vinvBubble, work) );
235 break;
236 }
237 default: {
238 INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
239 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): Operator type not implemented" );
240 }
241 }
242 }
243 }
244
245 // -------------------------------------------------------------------------------------
246 template<typename DT, typename OT, typename PT>
247 Basis_HCURL_QUAD_In_FEM<DT,OT,PT>::
248 Basis_HCURL_QUAD_In_FEM( const ordinal_type order,
249 const EPointType pointType ) {
250
251 INTREPID2_TEST_FOR_EXCEPTION( !(pointType == POINTTYPE_EQUISPACED ||
252 pointType == POINTTYPE_WARPBLEND), std::invalid_argument,
253 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): pointType must be either equispaced or warpblend.");
254
255 // this should be in host
256 Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> lineBasis( order, pointType );
257 Basis_HGRAD_LINE_Cn_FEM<DT,OT,PT> bubbleBasis( order - 1, POINTTYPE_GAUSS );
258
259 const ordinal_type
260 cardLine = lineBasis.getCardinality(),
261 cardBubble = bubbleBasis.getCardinality();
262
263 this->vinvLine_ = Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space>("Hcurl::Quad::In::vinvLine", cardLine, cardLine);
264 this->vinvBubble_ = Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space>("Hcurl::Quad::In::vinvBubble", cardBubble, cardBubble);
265
266 lineBasis.getVandermondeInverse(this->vinvLine_);
267 bubbleBasis.getVandermondeInverse(this->vinvBubble_);
268
269 this->basisCardinality_ = 2*cardLine*cardBubble;
270 this->basisDegree_ = order;
271 this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
272 this->basisType_ = BASIS_FEM_LAGRANGIAN;
273 this->basisCoordinates_ = COORDINATES_CARTESIAN;
274 this->functionSpace_ = FUNCTION_SPACE_HCURL;
275 pointType_ = pointType;
276
277 // initialize tags
278 {
279 // Basis-dependent initializations
280 const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
281 const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
282 const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
283 const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
284
285 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
286 constexpr ordinal_type maxCardLine = Parameters::MaxOrder + 1;
287 constexpr ordinal_type maxCardBubble = Parameters::MaxOrder;
288 ordinal_type tags[2*maxCardLine*maxCardBubble][4];
289
290 const ordinal_type edge_x[2] = {0,2};
291 const ordinal_type edge_y[2] = {3,1};
292 {
293 ordinal_type idx = 0;
294
298
299 // since there are x/y components in the interior
300 // dof sum should be computed before the information
301 // is assigned to tags
302 const ordinal_type
303 intr_ndofs_per_direction = (cardLine-2)*cardBubble,
304 intr_ndofs = 2*intr_ndofs_per_direction;
305
306 // x component (lineBasis(y) bubbleBasis(x))
307 for (ordinal_type j=0;j<cardLine;++j) { // y
308 const auto tag_y = lineBasis.getDofTag(j);
309 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
310 const auto tag_x = bubbleBasis.getDofTag(i);
311
312 if (tag_x(0) == 1 && tag_y(0) == 0) {
313 // edge: x edge, y vert
314 tags[idx][0] = 1; // edge dof
315 tags[idx][1] = edge_x[tag_y(1)];
316 tags[idx][2] = tag_x(2); // local dof id
317 tags[idx][3] = tag_x(3); // total number of dofs in this vertex
318 } else {
319 // interior
320 tags[idx][0] = 2; // interior dof
321 tags[idx][1] = 0;
322 tags[idx][2] = tag_x(2) + tag_x(3)*tag_y(2); // local dof id
323 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
324 }
325 }
326 }
327
328 // y component (bubbleBasis(y) lineBasis(x))
329 for (ordinal_type j=0;j<cardBubble;++j) { // y
330 const auto tag_y = bubbleBasis.getDofTag(j);
331 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
332 const auto tag_x = lineBasis.getDofTag(i);
333
334 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] = intr_ndofs_per_direction + tag_x(2) + tag_x(3)*tag_y(2); // local dof id
345 tags[idx][3] = intr_ndofs; // total number of dofs in this vertex
346 }
347 }
348 }
349 INTREPID2_TEST_FOR_EXCEPTION( idx != this->basisCardinality_ , std::runtime_error,
350 ">>> ERROR (Basis_HCURL_QUAD_In_FEM): " \
351 "counted tag index is not same as cardinality." );
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 // dofCoeffs on host and create its mirror view to device
373 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
374 dofCoeffsHost("dofCoeffsHost", this->basisCardinality_, this->basisCellTopology_.getDimension());
375
376 Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space>
377 dofCoordsLine("dofCoordsLine", cardLine, 1),
378 dofCoordsBubble("dofCoordsBubble", cardBubble, 1);
379
380 lineBasis.getDofCoords(dofCoordsLine);
381 auto dofCoordsLineHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsLine);
382 Kokkos::deep_copy(dofCoordsLineHost, dofCoordsLine);
383
384 bubbleBasis.getDofCoords(dofCoordsBubble);
385 auto dofCoordsBubbleHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), dofCoordsBubble);
386 Kokkos::deep_copy(dofCoordsBubbleHost, dofCoordsBubble);
387
388 {
389 ordinal_type idx = 0;
390
391 // x component (lineBasis(y) bubbleBasis(x))
392 for (ordinal_type j=0;j<cardLine;++j) { // y
393 for (ordinal_type i=0;i<cardBubble;++i,++idx) { // x
394 dofCoordsHost(idx,0) = dofCoordsBubbleHost(i,0);
395 dofCoordsHost(idx,1) = dofCoordsLineHost(j,0);
396 dofCoeffsHost(idx,0) = 1.0;
397 }
398 }
399
400 // y component (bubbleBasis(y) lineBasis(x))
401 for (ordinal_type j=0;j<cardBubble;++j) { // y
402 for (ordinal_type i=0;i<cardLine;++i,++idx) { // x
403 dofCoordsHost(idx,0) = dofCoordsLineHost(i,0);
404 dofCoordsHost(idx,1) = dofCoordsBubbleHost(j,0);
405 dofCoeffsHost(idx,1) = 1.0;
406 }
407 }
408 }
409
410 this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoordsHost);
411 Kokkos::deep_copy(this->dofCoords_, dofCoordsHost);
412
413 this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffsHost);
414 Kokkos::deep_copy(this->dofCoeffs_, dofCoeffsHost);
415 }
416
417}
418
419#endif
Intrepid2::Basis_HCURL_QUAD_In_FEM::Basis_HCURL_QUAD_In_FEM
Basis_HCURL_QUAD_In_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Definition: Intrepid2_HCURL_QUAD_In_FEMDef.hpp:248
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::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