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