Intrepid2
Intrepid2_HGRAD_LINE_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_LINE_CN_FEM_DEF_HPP__
50 #define __INTREPID2_HGRAD_LINE_CN_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_HGRAD_LINE_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 card = vinv.extent(0);
73  const ordinal_type npts = input.extent(0);
74 
75  const ordinal_type order = card - 1;
76  const double alpha = 0.0, beta = 0.0;
77 
78  typedef typename Kokkos::DynRankView<typename workViewType::value_type, typename workViewType::memory_space> viewType;
79  auto vcprop = Kokkos::common_view_alloc_prop(work);
80 
81  switch (opType) {
82  case OPERATOR_VALUE: {
83  viewType phis(Kokkos::view_wrap(work.data(), vcprop), card, npts);
84 
85  Impl::Basis_HGRAD_LINE_Cn_FEM_JACOBI::
86  Serial<opType>::getValues(phis, input, order, alpha, beta);
87 
88  for (ordinal_type i=0;i<card;++i)
89  for (ordinal_type j=0;j<npts;++j) {
90  output.access(i,j) = 0.0;
91  for (ordinal_type k=0;k<card;++k)
92  output.access(i,j) += vinv(k,i)*phis.access(k,j);
93  }
94  break;
95  }
96  case OPERATOR_GRAD:
97  case OPERATOR_D1:
98  case OPERATOR_D2:
99  case OPERATOR_D3:
100  case OPERATOR_D4:
101  case OPERATOR_D5:
102  case OPERATOR_D6:
103  case OPERATOR_D7:
104  case OPERATOR_D8:
105  case OPERATOR_D9:
106  case OPERATOR_D10:
107  opDn = getOperatorOrder(opType);
108  case OPERATOR_Dn: {
109  // dkcard is always 1 for 1D element
110  const ordinal_type dkcard = 1;
111  viewType phis(Kokkos::view_wrap(work.data(), vcprop), card, npts, dkcard);
112  Impl::Basis_HGRAD_LINE_Cn_FEM_JACOBI::
113  Serial<opType>::getValues(phis, input, order, alpha, beta, opDn);
114 
115  for (ordinal_type i=0;i<card;++i)
116  for (ordinal_type j=0;j<npts;++j)
117  for (ordinal_type k=0;k<dkcard;++k) {
118  output.access(i,j,k) = 0.0;
119  for (ordinal_type l=0;l<card;++l)
120  output.access(i,j,k) += vinv(l,i)*phis.access(l,j,k);
121  }
122  break;
123  }
124  default: {
125  INTREPID2_TEST_FOR_ABORT( true,
126  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_Cn_FEM::Serial::getValues) operator is not supported." );
127  }
128  }
129  }
130 
131 
132  template<typename DT, ordinal_type numPtsPerEval,
133  typename outputValueValueType, class ...outputValueProperties,
134  typename inputPointValueType, class ...inputPointProperties,
135  typename vinvValueType, class ...vinvProperties>
136  void
137  Basis_HGRAD_LINE_Cn_FEM::
138  getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
139  const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
140  const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
141  const EOperator operatorType ) {
142  typedef Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValueViewType;
143  typedef Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPointViewType;
144  typedef Kokkos::DynRankView<vinvValueType, vinvProperties...> vinvViewType;
145  typedef typename ExecSpace<typename inputPointViewType::execution_space,typename DT::execution_space>::ExecSpaceType ExecSpaceType;
146 
147  // loopSize corresponds to cardinality
148  const auto loopSizeTmp1 = (inputPoints.extent(0)/numPtsPerEval);
149  const auto loopSizeTmp2 = (inputPoints.extent(0)%numPtsPerEval != 0);
150  const auto loopSize = loopSizeTmp1 + loopSizeTmp2;
151  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
152 
153  typedef typename inputPointViewType::value_type inputPointType;
154 
155  const ordinal_type cardinality = outputValues.extent(0);
156 
157  auto vcprop = Kokkos::common_view_alloc_prop(inputPoints);
158  typedef typename Kokkos::DynRankView< inputPointType, typename inputPointViewType::memory_space> workViewType;
159  workViewType work(Kokkos::view_alloc("Basis_HGRAD_LINE_Cn_FEM::getValues::work", vcprop), cardinality, inputPoints.extent(0));
160 
161  switch (operatorType) {
162  case OPERATOR_VALUE: {
163  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
164  OPERATOR_VALUE,numPtsPerEval> FunctorType;
165  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work) );
166  break;
167  }
168  case OPERATOR_GRAD:
169  case OPERATOR_D1:
170  case OPERATOR_D2:
171  case OPERATOR_D3:
172  case OPERATOR_D4:
173  case OPERATOR_D5:
174  case OPERATOR_D6:
175  case OPERATOR_D7:
176  case OPERATOR_D8:
177  case OPERATOR_D9:
178  case OPERATOR_D10: {
179  typedef Functor<outputValueViewType,inputPointViewType,vinvViewType,workViewType,
180  OPERATOR_Dn,numPtsPerEval> FunctorType;
181  Kokkos::parallel_for( policy, FunctorType(outputValues, inputPoints, vinv, work,
182  getOperatorOrder(operatorType)) );
183  break;
184  }
185  default: {
186  INTREPID2_TEST_FOR_EXCEPTION( true , std::invalid_argument,
187  ">>> ERROR (Basis_HGRAD_LINE_Cn_FEM): Operator type not implemented" );
188  //break; commented out because this always throws
189  }
190  }
191  }
192  }
193 
194  // -------------------------------------------------------------------------------------
195  template<typename DT, typename OT, typename PT>
197  Basis_HGRAD_LINE_Cn_FEM( const ordinal_type order,
198  const EPointType pointType ) {
199  this->pointType_ = pointType;
200  this->basisCardinality_ = order+1;
201  this->basisDegree_ = order;
202  this->basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >() );
203  this->basisType_ = BASIS_FEM_LAGRANGIAN;
204  this->basisCoordinates_ = COORDINATES_CARTESIAN;
205  this->functionSpace_ = FUNCTION_SPACE_HGRAD;
206 
207  const ordinal_type card = this->basisCardinality_;
208 
209  // points are computed in the host and will be copied
210  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
211  dofCoords("Hgrad::Line::Cn::dofCoords", card, 1);
212 
213  //Default is Equispaced
214  auto pointT = (pointType == POINTTYPE_DEFAULT) ? POINTTYPE_EQUISPACED : pointType;
215 
216  switch (pointT) {
217  case POINTTYPE_EQUISPACED:
218  case POINTTYPE_WARPBLEND: {
219  // lattice ordering
220  {
221  const ordinal_type offset = 0;
222  PointTools::getLattice( dofCoords,
223  this->basisCellTopology_,
224  order, offset,
225  pointT );
226 
227  }
228  // topological order
229  // {
230  // // two vertices
231  // dofCoords(0,0) = -1.0;
232  // dofCoords(1,0) = 1.0;
233 
234  // // internal points
235  // typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
236  // auto pts = Kokkos::subview(dofCoords, range_type(2, card), Kokkos::ALL());
237 
238  // const auto offset = 1;
239  // PointTools::getLattice( pts,
240  // this->basisCellTopology_,
241  // order, offset,
242  // pointType );
243  // }
244  break;
245  }
246  case POINTTYPE_GAUSS: {
247  // internal points only
248  PointTools::getGaussPoints( dofCoords,
249  order );
250  break;
251  }
252  default: {
253  INTREPID2_TEST_FOR_EXCEPTION( !isValidPointType(pointT),
254  std::invalid_argument ,
255  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_Cn_FEM) invalid pointType." );
256  }
257  }
258 
259  this->dofCoords_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoords);
260  Kokkos::deep_copy(this->dofCoords_, dofCoords);
261 
262  // form Vandermonde matrix; actually, this is the transpose of the VDM,
263  // this matrix is used in LAPACK so it should be column major and left layout
264  const ordinal_type lwork = card*card;
265  Kokkos::DynRankView<typename ScalarViewType::value_type,Kokkos::LayoutLeft,Kokkos::HostSpace>
266  vmat("Hgrad::Line::Cn::vmat", card, card),
267  work("Hgrad::Line::Cn::work", lwork),
268  ipiv("Hgrad::Line::Cn::ipiv", card);
269 
270  const double alpha = 0.0, beta = 0.0;
271  Impl::Basis_HGRAD_LINE_Cn_FEM_JACOBI::
272  getValues<Kokkos::HostSpace::execution_space,Parameters::MaxNumPtsPerBasisEval>
273  (vmat, dofCoords, order, alpha, beta, OPERATOR_VALUE);
274 
275  ordinal_type info = 0;
276  Teuchos::LAPACK<ordinal_type,typename ScalarViewType::value_type> lapack;
277 
278  lapack.GETRF(card, card,
279  vmat.data(), vmat.stride_1(),
280  (ordinal_type*)ipiv.data(),
281  &info);
282 
283  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
284  std::runtime_error ,
285  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_Cn_FEM) lapack.GETRF returns nonzero info." );
286 
287  lapack.GETRI(card,
288  vmat.data(), vmat.stride_1(),
289  (ordinal_type*)ipiv.data(),
290  work.data(), lwork,
291  &info);
292 
293  INTREPID2_TEST_FOR_EXCEPTION( info != 0,
294  std::runtime_error ,
295  ">>> ERROR: (Intrepid2::Basis_HGRAD_LINE_Cn_FEM) lapack.GETRI returns nonzero info." );
296 
297  // create host mirror
298  Kokkos::DynRankView<typename ScalarViewType::value_type,typename DT::execution_space::array_layout,Kokkos::HostSpace>
299  vinv("Hgrad::Line::Cn::vinv", card, card);
300 
301  for (ordinal_type i=0;i<card;++i)
302  for (ordinal_type j=0;j<card;++j)
303  vinv(i,j) = vmat(j,i);
304 
305  this->vinv_ = Kokkos::create_mirror_view(typename DT::memory_space(), vinv);
306  Kokkos::deep_copy(this->vinv_ , vinv);
307 
308  // initialize tags
309  {
310  const bool is_vertex_included = (pointType != POINTTYPE_GAUSS);
311 
312  // Basis-dependent initializations
313  const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
314  const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
315  const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
316  const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
317 
318 
319  ordinal_type tags[Parameters::MaxOrder+1][4];
320 
321  // now we check the points for association
322  if (is_vertex_included) {
323  // lattice order
324  {
325  const auto v0 = 0;
326  tags[v0][0] = 0; // vertex dof
327  tags[v0][1] = 0; // vertex id
328  tags[v0][2] = 0; // local dof id
329  tags[v0][3] = 1; // total number of dofs in this vertex
330 
331  const ordinal_type iend = card - 2;
332  for (ordinal_type i=0;i<iend;++i) {
333  const auto e = i + 1;
334  tags[e][0] = 1; // edge dof
335  tags[e][1] = 0; // edge id
336  tags[e][2] = i; // local dof id
337  tags[e][3] = iend; // total number of dofs in this edge
338  }
339 
340  const auto v1 = card -1;
341  tags[v1][0] = 0; // vertex dof
342  tags[v1][1] = 1; // vertex id
343  tags[v1][2] = 0; // local dof id
344  tags[v1][3] = 1; // total number of dofs in this vertex
345  }
346 
347  // topological order
348  // {
349  // tags[0][0] = 0; // vertex dof
350  // tags[0][1] = 0; // vertex id
351  // tags[0][2] = 0; // local dof id
352  // tags[0][3] = 1; // total number of dofs in this vertex
353 
354  // tags[1][0] = 0; // vertex dof
355  // tags[1][1] = 1; // vertex id
356  // tags[1][2] = 0; // local dof id
357  // tags[1][3] = 1; // total number of dofs in this vertex
358 
359  // const ordinal_type iend = card - 2;
360  // for (ordinal_type i=0;i<iend;++i) {
361  // const auto ii = i + 2;
362  // tags[ii][0] = 1; // edge dof
363  // tags[ii][1] = 0; // edge id
364  // tags[ii][2] = i; // local dof id
365  // tags[ii][3] = iend; // total number of dofs in this edge
366  // }
367  // }
368  } else {
369  for (ordinal_type i=0;i<card;++i) {
370  tags[i][0] = 1; // edge dof
371  tags[i][1] = 0; // edge id
372  tags[i][2] = i; // local dof id
373  tags[i][3] = card; // total number of dofs in this edge
374  }
375  }
376 
377  OrdinalTypeArray1DHost tagView(&tags[0][0], card*4);
378 
379  // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays:
380  // tags are constructed on host
381  this->setOrdinalTagData(this->tagToOrdinal_,
382  this->ordinalToTag_,
383  tagView,
384  this->basisCardinality_,
385  tagSize,
386  posScDim,
387  posScOrd,
388  posDfOrd);
389  }
390  }
391 
392 }// namespace Intrepid2
393 
394 #endif
395 
396 
397 
398 
399 
400 
401 
402 
403 
404 
405 
406 
407 
408 
409 
KOKKOS_FORCEINLINE_FUNCTION bool isValidPointType(const EPointType pointType)
Verifies validity of a point type enum.
static void getGaussPoints(Kokkos::DynRankView< pointValueType, pointProperties... > points, const ordinal_type order)
EPointType
Enumeration of types of point distributions in Intrepid.
static void getLattice(Kokkos::DynRankView< pointValueType, pointProperties... > points, const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference simplex, quadrilateral or hexahedron (cu...
KOKKOS_INLINE_FUNCTION ordinal_type getOperatorOrder(const EOperator operatorType)
Returns order of an operator.
Basis_HGRAD_LINE_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
static constexpr ordinal_type MaxOrder
The maximum reconstruction order.