Intrepid2
Intrepid2_Basis.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_BASIS_HPP__
50 #define __INTREPID2_BASIS_HPP__
51 
52 #include "Intrepid2_ConfigDefs.hpp"
53 #include "Intrepid2_Types.hpp"
54 #include "Intrepid2_Utils.hpp"
55 
57 #include "Shards_CellTopology.hpp"
58 
59 namespace Intrepid2 {
60 
90  template<typename ExecSpaceType = void,
91  typename outputValueType = double,
92  typename pointValueType = double>
93  class Basis {
94  public:
95 
98  typedef Kokkos::View<ordinal_type,ExecSpaceType> ordinal_view_type;
99 
102  typedef Kokkos::View<EBasis,ExecSpaceType> ebasis_view_type;
103 
106  typedef Kokkos::View<ECoordinates,ExecSpaceType> ecoordinates_view_type;
107 
108  // ** tag interface
109  // - tag interface is not decorated with Kokkos inline so it should be allocated on hostspace
110 
113  typedef Kokkos::View<ordinal_type* ,typename ExecSpaceType::array_layout,Kokkos::HostSpace> ordinal_type_array_1d_host;
114 
117  typedef Kokkos::View<ordinal_type** ,typename ExecSpaceType::array_layout,Kokkos::HostSpace> ordinal_type_array_2d_host;
118 
121  typedef Kokkos::View<ordinal_type***,typename ExecSpaceType::array_layout,Kokkos::HostSpace> ordinal_type_array_3d_host;
122 
125  typedef Kokkos::View<ordinal_type* , Kokkos::LayoutStride, Kokkos::HostSpace> ordinal_type_array_stride_1d_host;
126 
129  typedef Kokkos::View<ordinal_type* ,ExecSpaceType> ordinal_type_array_1d;
130 
133  typedef Kokkos::View<ordinal_type** ,ExecSpaceType> ordinal_type_array_2d;
134 
137  typedef Kokkos::View<ordinal_type***,ExecSpaceType> ordinal_type_array_3d;
138 
141  typedef Kokkos::View<ordinal_type* , Kokkos::LayoutStride, ExecSpaceType> ordinal_type_array_stride_1d;
142 
145  typedef typename ScalarTraits<pointValueType>::scalar_type scalarType;
146 
147  protected:
148 
151  ordinal_type basisCardinality_;
152 
155  ordinal_type basisDegree_;
156 
161  shards::CellTopology basisCellTopology_;
162 
166 
170 
173  //Kokkos::View<bool,ExecSpaceType> basisTagsAreSet_;
174 
187 
200 
212  template<typename OrdinalTypeView3D,
213  typename OrdinalTypeView2D,
214  typename OrdinalTypeView1D>
215  void setOrdinalTagData( OrdinalTypeView3D &tagToOrdinal,
216  OrdinalTypeView2D &ordinalToTag,
217  const OrdinalTypeView1D tags,
218  const ordinal_type basisCard,
219  const ordinal_type tagSize,
220  const ordinal_type posScDim,
221  const ordinal_type posScOrd,
222  const ordinal_type posDfOrd ) {
223  // Create ordinalToTag
224  ordinalToTag = OrdinalTypeView2D("ordinalToTag", basisCard, tagSize);
225 
226  // Initialize with -1
227  Kokkos::deep_copy( ordinalToTag, -1 );
228 
229  // Copy tags
230  for (ordinal_type i=0;i<basisCard;++i)
231  for (ordinal_type j=0;j<tagSize;++j)
232  ordinalToTag(i, j) = tags(i*tagSize + j);
233 
234  // Find out dimension of tagToOrdinal
235  auto maxScDim = 0; // first dimension of tagToOrdinal
236  for (ordinal_type i=0;i<basisCard;++i)
237  if (maxScDim < tags(i*tagSize + posScDim))
238  maxScDim = tags(i*tagSize + posScDim);
239  ++maxScDim;
240 
241  auto maxScOrd = 0; // second dimension of tagToOrdinal
242  for (ordinal_type i=0;i<basisCard;++i)
243  if (maxScOrd < tags(i*tagSize + posScOrd))
244  maxScOrd = tags(i*tagSize + posScOrd);
245  ++maxScOrd;
246 
247  auto maxDfOrd = 0; // third dimension of tagToOrdinal
248  for (ordinal_type i=0;i<basisCard;++i)
249  if (maxDfOrd < tags(i*tagSize + posDfOrd))
250  maxDfOrd = tags(i*tagSize + posDfOrd);
251  ++maxDfOrd;
252 
253  // Create tagToOrdinal
254  tagToOrdinal = OrdinalTypeView3D("tagToOrdinal", maxScDim, maxScOrd, maxDfOrd);
255 
256  // Initialize with -1
257  Kokkos::deep_copy( tagToOrdinal, -1 );
258 
259  // Overwrite elements of the array corresponding to tags with local DoF Id's, leave all other = -1
260  for (ordinal_type i=0;i<basisCard;++i)
261  tagToOrdinal(tags(i*tagSize), tags(i*tagSize+1), tags(i*tagSize+2)) = i;
262  }
263 
264  // dof coords
267  Kokkos::DynRankView<scalarType,ExecSpaceType> dofCoords_;
268 
269  // dof coeffs
277  Kokkos::DynRankView<scalarType,ExecSpaceType> dofCoeffs_;
278 
279  public:
280 
281  Basis() = default;
282  virtual~Basis() = default;
283 
284  // receives input arguments
287  outputValueType getDummyOutputValue() { return outputValueType(); }
288 
291  pointValueType getDummyPointValue() { return pointValueType(); }
292 
295  typedef Kokkos::DynRankView<outputValueType,Kokkos::LayoutStride,ExecSpaceType> outputViewType;
296 
299  typedef Kokkos::DynRankView<pointValueType,Kokkos::LayoutStride,ExecSpaceType> pointViewType;
300 
303  typedef Kokkos::DynRankView<scalarType,Kokkos::LayoutStride,ExecSpaceType> scalarViewType;
304 
323  virtual
324  void
325  getValues( outputViewType outputValues,
326  const pointViewType inputPoints,
327  const EOperator operatorType = OPERATOR_VALUE ) const {
328  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
329  ">>> ERROR (Basis::getValues): this method (FEM) is not supported or should be over-riden accordingly by derived classes.");
330  }
331 
351  virtual
352  void
353  getValues( outputViewType outputValues,
354  const pointViewType inputPoints,
355  const pointViewType cellVertices,
356  const EOperator operatorType = OPERATOR_VALUE ) const {
357  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
358  ">>> ERROR (Basis::getValues): this method (FVM) is not supported or should be over-riden accordingly by derived classes.");
359  }
360 
361 
365  virtual
366  void
367  getDofCoords( scalarViewType dofCoords ) const {
368  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
369  ">>> ERROR (Basis::getDofCoords): this method is not supported or should be over-riden accordingly by derived classes.");
370  }
371 
380  virtual
381  void
382  getDofCoeffs( scalarViewType dofCoeffs ) const {
383  INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
384  ">>> ERROR (Basis::getDofCoeffs): this method is not supported or should be over-riden accordingly by derived classes.");
385  }
386 
387 
392  virtual
393  const char*
394  getName() const {
395  return "Intrepid2_Basis";
396  }
397 
400  virtual
401  bool
403  return false;
404  }
405 
410  ordinal_type
411  getCardinality() const {
412  return basisCardinality_;
413  }
414 
415 
420  ordinal_type
421  getDegree() const {
422  return basisDegree_;
423  }
424 
425 
431  shards::CellTopology
433  return basisCellTopology_;
434  }
435 
436 
441  EBasis
442  getBasisType() const {
443  return basisType_;
444  }
445 
446 
453  return basisCoordinates_;
454  }
455 
463  ordinal_type
464  getDofCount( const ordinal_type subcDim,
465  const ordinal_type subcOrd ) const {
466  if ( subcDim >= 0 && subcDim < static_cast<ordinal_type>(tagToOrdinal_.extent(0)) &&
467  subcOrd >= 0 && subcOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(1)) )
468  {
469  int firstDofOrdinal = tagToOrdinal_(subcDim, subcOrd, 0); // will be -1 if no dofs for subcell
470  if (firstDofOrdinal == -1) return static_cast<ordinal_type>(0);
471  // otherwise, lookup its tag and return the dof count stored there
472  return static_cast<ordinal_type>(this->getDofTag(firstDofOrdinal)[3]);
473  }
474  else
475  {
476  // subcDim and/or subcOrd out of bounds -> no dofs associated with subcell
477  return static_cast<ordinal_type>(0);
478  }
479  }
480 
489  ordinal_type
490  getDofOrdinal( const ordinal_type subcDim,
491  const ordinal_type subcOrd,
492  const ordinal_type subcDofOrd ) const {
493  // this should be abort and able to be called as a device function
494 #ifdef HAVE_INTREPID2_DEBUG
495  INTREPID2_TEST_FOR_EXCEPTION( subcDim < 0 || subcDim >= static_cast<ordinal_type>(tagToOrdinal_.extent(0)), std::out_of_range,
496  ">>> ERROR (Basis::getDofOrdinal): subcDim is out of range");
497  INTREPID2_TEST_FOR_EXCEPTION( subcOrd < 0 || subcOrd >= static_cast<ordinal_type>(tagToOrdinal_.extent(1)), std::out_of_range,
498  ">>> ERROR (Basis::getDofOrdinal): subcOrd is out of range");
499  INTREPID2_TEST_FOR_EXCEPTION( subcDofOrd < 0 || subcDofOrd >= static_cast<ordinal_type>(tagToOrdinal_.extent(2)), std::out_of_range,
500  ">>> ERROR (Basis::getDofOrdinal): subcDofOrd is out of range");
501 #endif
502  ordinal_type r_val = -1;
503  if ( subcDim < static_cast<ordinal_type>(tagToOrdinal_.extent(0)) &&
504  subcOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(1)) &&
505  subcDofOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(2)) )
506  r_val = tagToOrdinal_(subcDim, subcOrd, subcDofOrd);
507 #ifdef HAVE_INTREPID2_DEBUG
508  INTREPID2_TEST_FOR_EXCEPTION( r_val == -1, std::runtime_error,
509  ">>> ERROR (Basis::getDofOrdinal): Invalid DoF tag is found.");
510 #endif
511  return r_val;
512  }
513 
517  return tagToOrdinal_;
518  }
519 
520 
532  getDofTag( const ordinal_type dofOrd ) const {
533 #ifdef HAVE_INTREPID2_DEBUG
534  INTREPID2_TEST_FOR_EXCEPTION( dofOrd < 0 || dofOrd >= static_cast<ordinal_type>(ordinalToTag_.extent(0)), std::out_of_range,
535  ">>> ERROR (Basis::getDofTag): dofOrd is out of range");
536 #endif
537  return Kokkos::subview(ordinalToTag_, dofOrd, Kokkos::ALL());
538  }
539 
540 
550  getAllDofTags() const {
551  return ordinalToTag_;
552  }
553 
554  }; // class Basis
555 
556 
557  //--------------------------------------------------------------------------------------------//
558  // //
559  // Helper functions of the Basis class //
560  // //
561  //--------------------------------------------------------------------------------------------//
562 
563  //
564  // functions for orders, cardinality, etc.
565  //
566 
567 
579  KOKKOS_INLINE_FUNCTION
580  ordinal_type getFieldRank(const EFunctionSpace spaceType);
581 
617  KOKKOS_INLINE_FUNCTION
618  ordinal_type getOperatorRank(const EFunctionSpace spaceType,
619  const EOperator operatorType,
620  const ordinal_type spaceDim);
621 
627  KOKKOS_INLINE_FUNCTION
628  ordinal_type getOperatorOrder(const EOperator operatorType);
629 
630  template<EOperator operatorType>
631  KOKKOS_INLINE_FUNCTION
632  constexpr ordinal_type getOperatorOrder();
633 
657  template<ordinal_type spaceDim>
658  KOKKOS_INLINE_FUNCTION
659  ordinal_type getDkEnumeration(const ordinal_type xMult,
660  const ordinal_type yMult = -1,
661  const ordinal_type zMult = -1);
662 
663 
674  template<ordinal_type spaceDim>
675  KOKKOS_INLINE_FUNCTION
676  ordinal_type getPnEnumeration(const ordinal_type p,
677  const ordinal_type q = 0,
678  const ordinal_type r = 0);
679 
680 
681 
700 template<typename value_type>
701 KOKKOS_INLINE_FUNCTION
702 void getJacobyRecurrenceCoeffs (
703  value_type &an,
704  value_type &bn,
705  value_type &cn,
706  const ordinal_type alpha,
707  const ordinal_type beta ,
708  const ordinal_type n);
709 
710 
711 
712 
713 
714  // /** \brief Returns multiplicities of dx, dy, and dz based on the enumeration of the partial
715  // derivative, its order and the space dimension. Inverse of the getDkEnumeration() method.
716 
717  // \param partialMult [out] - array with the multiplicities f dx, dy and dz
718  // \param derivativeEnum [in] - enumeration of the partial derivative
719  // \param operatorType [in] - k-th partial derivative Dk
720  // \param spaceDim [in] - space dimension
721  // */
722  // template<typename OrdinalArrayType>
723  // KOKKOS_INLINE_FUNCTION
724  // void getDkMultiplicities(OrdinalArrayType partialMult,
725  // const ordinal_type derivativeEnum,
726  // const EOperator operatorType,
727  // const ordinal_type spaceDim);
728 
747  KOKKOS_INLINE_FUNCTION
748  ordinal_type getDkCardinality(const EOperator operatorType,
749  const ordinal_type spaceDim);
750 
751  template<EOperator operatorType, ordinal_type spaceDim>
752  KOKKOS_INLINE_FUNCTION
753  constexpr ordinal_type getDkCardinality();
754 
755 
756 
766  template<ordinal_type spaceDim>
767  KOKKOS_INLINE_FUNCTION
768  ordinal_type getPnCardinality (ordinal_type n);
769 
770  template<ordinal_type spaceDim, ordinal_type n>
771  KOKKOS_INLINE_FUNCTION
772  constexpr ordinal_type getPnCardinality ();
773 
774 
775 
776  //
777  // Argument check
778  //
779 
780 
791  template<typename outputValueViewType,
792  typename inputPointViewType>
793  void getValues_HGRAD_Args( const outputValueViewType outputValues,
794  const inputPointViewType inputPoints,
795  const EOperator operatorType,
796  const shards::CellTopology cellTopo,
797  const ordinal_type basisCard );
798 
809  template<typename outputValueViewType,
810  typename inputPointViewType>
811  void getValues_HCURL_Args( const outputValueViewType outputValues,
812  const inputPointViewType inputPoints,
813  const EOperator operatorType,
814  const shards::CellTopology cellTopo,
815  const ordinal_type basisCard );
816 
827  template<typename outputValueViewType,
828  typename inputPointViewType>
829  void getValues_HDIV_Args( const outputValueViewType outputValues,
830  const inputPointViewType inputPoints,
831  const EOperator operatorType,
832  const shards::CellTopology cellTopo,
833  const ordinal_type basisCard );
834 
845  template<typename outputValueViewType,
846  typename inputPointViewType>
847  void getValues_HVOL_Args( const outputValueViewType outputValues,
848  const inputPointViewType inputPoints,
849  const EOperator operatorType,
850  const shards::CellTopology cellTopo,
851  const ordinal_type basisCard );
852 
853 }// namespace Intrepid2
854 
855 #include <Intrepid2_BasisDef.hpp>
856 
857 //--------------------------------------------------------------------------------------------//
858 // //
859 // D O C U M E N T A T I O N P A G E S //
860 // //
861 //--------------------------------------------------------------------------------------------//
963 #endif
Kokkos::View< ordinal_type *, Kokkos::LayoutStride, ExecSpaceType > ordinal_type_array_stride_1d
View type for 1d device array.
Kokkos::View< ordinal_type ***, ExecSpaceType > ordinal_type_array_3d
View type for 3d device array.
const ordinal_type_array_2d_host getAllDofTags() const
Retrieves all DoF tags.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, ExecSpaceType > scalarViewType
View type for scalars.
Kokkos::View< ordinal_type **,ExecSpaceType > ordinal_type_array_2d
View type for 2d device array.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
ordinal_type getCardinality() const
Returns cardinality of the basis.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Kokkos::DynRankView< outputValueType, Kokkos::LayoutStride, ExecSpaceType > outputViewType
View type for basis value output.
Kokkos::DynRankView< scalarType, ExecSpaceType > dofCoords_
Coordinates of degrees-of-freedom for basis functions defined in physical space.
const ordinal_type_array_stride_1d_host getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
virtual void getDofCoords(scalarViewType dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation https://trili...
ECoordinates
Enumeration of coordinate systems for geometrical entities (cells, points).
EBasis getBasisType() const
Returns the basis type.
Kokkos::View< ordinal_type ***, typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_3d_host
View type for 3d host array.
outputValueType getDummyOutputValue()
Dummy array to receive input arguments.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
Kokkos::View< ECoordinates, ExecSpaceType > ecoordinates_view_type
View for coordinate system type.
Header function for Intrepid2::Util class and other utility functions.
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
Kokkos::View< ordinal_type *, Kokkos::LayoutStride, Kokkos::HostSpace > ordinal_type_array_stride_1d_host
View type for 1d host array.
ordinal_type_array_2d_host ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual bool requireOrientation() const
True if orientation is required.
const ordinal_type_array_3d_host getAllDofOrdinal() const
DoF tag to ordinal data structure.
EBasis basisType_
Type of the basis.
EBasis
Enumeration of basis types for discrete spaces in Intrepid.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
Kokkos::View< ordinal_type *,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_1d_host
View type for 1d host array.
Definition of cell topology information.
Kokkos::View< ordinal_type *,ExecSpaceType > ordinal_type_array_1d
View type for 1d device array.
Kokkos::DynRankView< pointValueType, Kokkos::LayoutStride, ExecSpaceType > pointViewType
View type for input points.
Contains definitions of custom data types in Intrepid2.
void setOrdinalTagData(OrdinalTypeView3D &tagToOrdinal, OrdinalTypeView2D &ordinalToTag, const OrdinalTypeView1D tags, const ordinal_type basisCard, const ordinal_type tagSize, const ordinal_type posScDim, const ordinal_type posScOrd, const ordinal_type posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
ordinal_type_array_3d_host tagToOrdinal_
DoF tag to ordinal lookup table.
ordinal_type getDofOrdinal(const ordinal_type subcDim, const ordinal_type subcOrd, const ordinal_type subcDofOrd) const
DoF tag to ordinal lookup.
ordinal_type getDegree() const
Returns the degree of the basis.
virtual void getValues(outputViewType outputValues, const pointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
ECoordinates getCoordinateSystem() const
Returns the type of coordinate system for which the basis is defined.
Kokkos::View< EBasis, ExecSpaceType > ebasis_view_type
View for basis type.
Kokkos::DynRankView< scalarType, ExecSpaceType > dofCoeffs_
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...
Kokkos::View< ordinal_type, ExecSpaceType > ordinal_view_type
View type for ordinal.
virtual void getValues(outputViewType outputValues, const pointViewType inputPoints, const pointViewType cellVertices, const EOperator operatorType=OPERATOR_VALUE) const
Evaluation of an FVD basis evaluation on a physical cell.
Implementation file for the abstract base class Intrepid2::Basis.
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
virtual void getDofCoeffs(scalarViewType dofCoeffs) const
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...
pointValueType getDummyPointValue()
Dummy array to receive input arguments.
virtual const char * getName() const
Returns basis name.
Kokkos::View< ordinal_type **,typename ExecSpaceType::array_layout, Kokkos::HostSpace > ordinal_type_array_2d_host
View type for 2d host array.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...