50 #ifndef __INTREPID2_BASIS_HPP__ 51 #define __INTREPID2_BASIS_HPP__ 53 #include "Intrepid2_ConfigDefs.hpp" 60 #include "Kokkos_Vector.hpp" 61 #include "Shards_CellTopology.hpp" 62 #include <Teuchos_RCPDecl.hpp> 68 template<
typename DeviceType = void,
69 typename OutputType = double,
70 typename PointType =
double>
75 template <
typename DeviceType =
void,
typename OutputType =
double,
typename Po
intType =
double>
76 using BasisPtr = Teuchos::RCP<Basis<DeviceType,OutputType,PointType> >;
80 template <
typename OutputType =
double,
typename Po
intType =
double>
119 template<
typename Device,
120 typename outputValueType,
121 typename pointValueType>
190 typedef typename ScalarTraits<pointValueType>::scalar_type
scalarType;
260 template<
typename OrdinalTypeView3D,
261 typename OrdinalTypeView2D,
262 typename OrdinalTypeView1D>
264 OrdinalTypeView2D &ordinalToTag,
265 const OrdinalTypeView1D tags,
266 const ordinal_type basisCard,
267 const ordinal_type tagSize,
268 const ordinal_type posScDim,
269 const ordinal_type posScOrd,
270 const ordinal_type posDfOrd ) {
272 ordinalToTag = OrdinalTypeView2D(
"ordinalToTag", basisCard, tagSize);
275 Kokkos::deep_copy( ordinalToTag, -1 );
278 for (ordinal_type i=0;i<basisCard;++i)
279 for (ordinal_type j=0;j<tagSize;++j)
280 ordinalToTag(i, j) = tags(i*tagSize + j);
284 for (ordinal_type i=0;i<basisCard;++i)
285 if (maxScDim < tags(i*tagSize + posScDim))
286 maxScDim = tags(i*tagSize + posScDim);
290 for (ordinal_type i=0;i<basisCard;++i)
291 if (maxScOrd < tags(i*tagSize + posScOrd))
292 maxScOrd = tags(i*tagSize + posScOrd);
296 for (ordinal_type i=0;i<basisCard;++i)
297 if (maxDfOrd < tags(i*tagSize + posDfOrd))
298 maxDfOrd = tags(i*tagSize + posDfOrd);
302 tagToOrdinal = OrdinalTypeView3D(
"tagToOrdinal", maxScDim, maxScOrd, maxDfOrd);
305 Kokkos::deep_copy( tagToOrdinal, -1 );
308 for (ordinal_type i=0;i<basisCard;++i)
309 tagToOrdinal(tags(i*tagSize), tags(i*tagSize+1), tags(i*tagSize+2)) = i;
338 virtual~
Basis() =
default;
351 using OutputViewType = Kokkos::DynRankView<OutputValueType,Kokkos::LayoutStride,DeviceType>;
355 using PointViewType = Kokkos::DynRankView<PointValueType,Kokkos::LayoutStride,DeviceType>;
359 using ScalarViewType = Kokkos::DynRankView<scalarType,Kokkos::LayoutStride,DeviceType>;
365 Kokkos::DynRankView<OutputValueType,DeviceType>
allocateOutputView(
const int numPoints,
const EOperator operatorType = OPERATOR_VALUE)
const;
374 const bool operatorIsDk = (operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10);
375 const bool operatorSupported = (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV) || operatorIsDk;
376 INTREPID2_TEST_FOR_EXCEPTION(!operatorSupported, std::invalid_argument,
"operator is not supported by allocateBasisValues");
388 bool useVectorData = (dataView.rank() == 3);
424 const EOperator = OPERATOR_VALUE )
const {
425 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
426 ">>> ERROR (Basis::getValues): this method (FEM) is not supported or should be overridden accordingly by derived classes.");
444 const EOperator operatorType = OPERATOR_VALUE )
const {
450 if (outputValues.numTensorDataFamilies() > 0)
452 INTREPID2_TEST_FOR_EXCEPTION(outputValues.
tensorData(0).numTensorComponents() != 1, std::invalid_argument,
"default implementation of getValues() only supports outputValues with trivial tensor-product structure");
453 outputData = outputValues.
tensorData().getTensorComponent(0);
457 INTREPID2_TEST_FOR_EXCEPTION(outputValues.
vectorData().numComponents() != 1, std::invalid_argument,
"default implementation of getValues() only supports outputValues with trivial tensor-product structure");
458 INTREPID2_TEST_FOR_EXCEPTION(outputValues.
vectorData().getComponent(0).numTensorComponents() != 1, std::invalid_argument,
"default implementation of getValues() only supports outputValues with trivial tensor-product structure");
459 outputData = outputValues.
vectorData().getComponent(0).getTensorComponent(0);
489 const EOperator = OPERATOR_VALUE )
const {
490 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
491 ">>> ERROR (Basis::getValues): this method (FVM) is not supported or should be overridden accordingly by derived classes.");
501 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
502 ">>> ERROR (Basis::getDofCoords): this method is not supported or should be overridden accordingly by derived classes.");
516 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
517 ">>> ERROR (Basis::getDofCoeffs): this method is not supported or should be overridden accordingly by derived classes.");
529 INTREPID2_TEST_FOR_EXCEPTION(
basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
530 ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
532 int requestedDegreeLength = degrees.extent_int(0);
533 INTREPID2_TEST_FOR_EXCEPTION(degreeEntryLength != requestedDegreeLength, std::invalid_argument,
"length of degrees does not match the entries in fieldOrdinalPolynomialDegree_");
534 std::vector<int> fieldOrdinalsVector;
538 for (
int d=0; d<degreeEntryLength; d++)
542 if (matches) fieldOrdinalsVector.push_back(basisOrdinal);
545 for (
unsigned i=0; i<fieldOrdinalsVector.size(); i++)
547 fieldOrdinals(i) = fieldOrdinalsVector[i];
549 return fieldOrdinals;
565 INTREPID2_TEST_FOR_EXCEPTION(
basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
566 ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
568 for (
unsigned d=0; d<degrees.size(); d++)
570 degreesView(d) = degrees[d];
573 std::vector<int> fieldOrdinalsVector(fieldOrdinalsView.extent_int(0));
574 for (
int i=0; i<fieldOrdinalsView.extent_int(0); i++)
576 fieldOrdinalsVector[i] = fieldOrdinalsView(i);
578 return fieldOrdinalsVector;
589 INTREPID2_TEST_FOR_EXCEPTION(
basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
590 ">>> ERROR (Basis::getPolynomialDegreeOfField): this method is not supported for non-hierarchical bases.");
591 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal < 0, std::invalid_argument,
"field ordinal must be non-negative");
592 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal >=
fieldOrdinalPolynomialDegree_.extent_int(0), std::invalid_argument,
"field ordinal out of bounds");
596 for (
int d=0; d<polyDegreeLength; d++)
612 INTREPID2_TEST_FOR_EXCEPTION(
basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
613 ">>> ERROR (Basis::getPolynomialDegreeOfFieldAsVector): this method is not supported for non-hierarchical bases.");
615 std::vector<int> polynomialDegree(polynomialDegreeView.extent_int(0));
617 for (
unsigned d=0; d<polynomialDegree.size(); d++)
619 polynomialDegree[d] = polynomialDegreeView(d);
621 return polynomialDegree;
628 INTREPID2_TEST_FOR_EXCEPTION(
basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
629 ">>> ERROR (Basis::getPolynomialDegreeLength): this method is not supported for non-hierarchical bases.");
640 return "Intrepid2_Basis";
718 const ordinal_type subcOrd )
const {
719 if ( subcDim >= 0 && subcDim < static_cast<ordinal_type>(
tagToOrdinal_.extent(0)) &&
720 subcOrd >= 0 && subcOrd < static_cast<ordinal_type>(
tagToOrdinal_.extent(1)) )
723 if (firstDofOrdinal == -1)
return static_cast<ordinal_type
>(0);
725 return static_cast<ordinal_type
>(this->
getDofTag(firstDofOrdinal)[3]);
730 return static_cast<ordinal_type
>(0);
744 const ordinal_type subcOrd,
745 const ordinal_type subcDofOrd )
const {
747 #ifdef HAVE_INTREPID2_DEBUG 748 INTREPID2_TEST_FOR_EXCEPTION( subcDim < 0 || subcDim >= static_cast<ordinal_type>(
tagToOrdinal_.extent(0)), std::out_of_range,
749 ">>> ERROR (Basis::getDofOrdinal): subcDim is out of range");
750 INTREPID2_TEST_FOR_EXCEPTION( subcOrd < 0 || subcOrd >= static_cast<ordinal_type>(
tagToOrdinal_.extent(1)), std::out_of_range,
751 ">>> ERROR (Basis::getDofOrdinal): subcOrd is out of range");
752 INTREPID2_TEST_FOR_EXCEPTION( subcDofOrd < 0 || subcDofOrd >= static_cast<ordinal_type>(
tagToOrdinal_.extent(2)), std::out_of_range,
753 ">>> ERROR (Basis::getDofOrdinal): subcDofOrd is out of range");
755 ordinal_type r_val = -1;
756 if ( subcDim < static_cast<ordinal_type>(
tagToOrdinal_.extent(0)) &&
757 subcOrd < static_cast<ordinal_type>(
tagToOrdinal_.extent(1)) &&
758 subcDofOrd < static_cast<ordinal_type>(
tagToOrdinal_.extent(2)) )
760 #ifdef HAVE_INTREPID2_DEBUG 761 INTREPID2_TEST_FOR_EXCEPTION( r_val == -1, std::runtime_error,
762 ">>> ERROR (Basis::getDofOrdinal): Invalid DoF tag is found.");
786 #ifdef HAVE_INTREPID2_DEBUG 787 INTREPID2_TEST_FOR_EXCEPTION( dofOrd < 0 || dofOrd >= static_cast<ordinal_type>(
ordinalToTag_.extent(0)), std::out_of_range,
788 ">>> ERROR (Basis::getDofTag): dofOrd is out of range");
790 return Kokkos::subview(
ordinalToTag_, dofOrd, Kokkos::ALL());
823 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
824 ">>> ERROR (Basis::getSubCellRefBasis): this method is not supported or should be overridden accordingly by derived classes.");
833 INTREPID2_TEST_FOR_EXCEPTION(
true, std::logic_error,
834 ">>> ERROR (Basis::getHostBasis): this method is not supported or should be overridden accordingly by derived classes.");
861 KOKKOS_INLINE_FUNCTION
862 ordinal_type getFieldRank(
const EFunctionSpace spaceType);
899 KOKKOS_INLINE_FUNCTION
900 ordinal_type getOperatorRank(
const EFunctionSpace spaceType,
901 const EOperator operatorType,
902 const ordinal_type spaceDim);
909 KOKKOS_INLINE_FUNCTION
910 ordinal_type getOperatorOrder(
const EOperator operatorType);
912 template<EOperator operatorType>
913 KOKKOS_INLINE_FUNCTION
914 constexpr ordinal_type getOperatorOrder();
939 template<ordinal_type spaceDim>
940 KOKKOS_INLINE_FUNCTION
941 ordinal_type getDkEnumeration(
const ordinal_type xMult,
942 const ordinal_type yMult = -1,
943 const ordinal_type zMult = -1);
956 template<ordinal_type spaceDim>
957 KOKKOS_INLINE_FUNCTION
958 ordinal_type getPnEnumeration(
const ordinal_type p,
959 const ordinal_type q = 0,
960 const ordinal_type r = 0);
982 template<
typename value_type>
983 KOKKOS_INLINE_FUNCTION
984 void getJacobyRecurrenceCoeffs (
988 const ordinal_type alpha,
989 const ordinal_type beta ,
990 const ordinal_type n);
1029 KOKKOS_INLINE_FUNCTION
1030 ordinal_type getDkCardinality(
const EOperator operatorType,
1031 const ordinal_type spaceDim);
1033 template<EOperator operatorType, ordinal_type spaceDim>
1034 KOKKOS_INLINE_FUNCTION
1035 constexpr ordinal_type getDkCardinality();
1048 template<ordinal_type spaceDim>
1049 KOKKOS_INLINE_FUNCTION
1050 ordinal_type getPnCardinality (ordinal_type n);
1052 template<ordinal_type spaceDim, ordinal_type n>
1053 KOKKOS_INLINE_FUNCTION
1054 constexpr ordinal_type getPnCardinality ();
1073 template<
typename outputValueViewType,
1074 typename inputPointViewType>
1075 void getValues_HGRAD_Args(
const outputValueViewType outputValues,
1076 const inputPointViewType inputPoints,
1077 const EOperator operatorType,
1078 const shards::CellTopology cellTopo,
1079 const ordinal_type basisCard );
1091 template<
typename outputValueViewType,
1092 typename inputPointViewType>
1093 void getValues_HCURL_Args(
const outputValueViewType outputValues,
1094 const inputPointViewType inputPoints,
1095 const EOperator operatorType,
1096 const shards::CellTopology cellTopo,
1097 const ordinal_type basisCard );
1109 template<
typename outputValueViewType,
1110 typename inputPointViewType>
1111 void getValues_HDIV_Args(
const outputValueViewType outputValues,
1112 const inputPointViewType inputPoints,
1113 const EOperator operatorType,
1114 const shards::CellTopology cellTopo,
1115 const ordinal_type basisCard );
1127 template<
typename outputValueViewType,
1128 typename inputPointViewType>
1129 void getValues_HVOL_Args(
const outputValueViewType outputValues,
1130 const inputPointViewType inputPoints,
1131 const EOperator operatorType,
1132 const shards::CellTopology cellTopo,
1133 const ordinal_type basisCard );
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
virtual void getDofCoords(ScalarViewType) const
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
OutputScalar OutputValueType
Output value type for basis; default is double.
std::vector< int > getFieldOrdinalsForDegree(std::vector< int > °rees) const
For hierarchical bases, returns the field ordinals that have at most the specified degree in each dim...
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
OrdinalTypeArray1DHost getPolynomialDegreeOfField(int fieldOrdinal) const
For hierarchical bases, returns the polynomial degree (which may have multiple values in higher spati...
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
Kokkos::View< EBasis, DeviceType > EBasisViewType
View for basis type.
Kokkos::DynRankView< OutputValueType, DeviceType > allocateOutputView(const int numPoints, const EOperator operatorType=OPERATOR_VALUE) const
Allocate a View container suitable for passing to the getValues() variant that accepts Kokkos DynRank...
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
ECoordinates getCoordinateSystem() const
Returns the type of coordinate system for which the basis is defined.
EFunctionSpace functionSpace_
The function space in which the basis is defined.
KOKKOS_INLINE_FUNCTION enable_if_t< rank==1, const Kokkos::View< typename RankExpander< DataScalar, rank >::value_type, DeviceType > & > getUnderlyingView() const
Returns the underlying view. Throws an exception if the underlying view is not rank 1...
OutputValueType getDummyOutputValue()
Dummy array to receive input arguments.
KOKKOS_INLINE_FUNCTION std::enable_if< std::is_integral< iType >::value, int >::type extent_int(const iType &r) const
Returns the logical extent in the requested dimension.
ECoordinates
Enumeration of coordinate systems for geometrical entities (cells, points).
Kokkos::View< ordinal_type ***, DeviceType > OrdinalTypeArray3D
View type for 3d device array.
virtual void getValues(OutputViewType, const PointViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of an FVD basis evaluation on a physical cell.
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
virtual void getDofCoeffs(ScalarViewType) const
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...
Header function for Intrepid2::Util class and other utility functions.
ordinal_type getCardinality() const
Returns cardinality of the basis.
PointScalar PointValueType
Point value type for basis; default is double.
Kokkos::View< ordinal_type ***, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray3DHost
View type for 3d host array.
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
BasisPtr< typename Kokkos::HostSpace::device_type, OutputType, PointType > HostBasisPtr
Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type), allowing host access to input and output views, and ensuring host execution of basis evaluation.
OrdinalTypeArray1DHost getFieldOrdinalsForDegree(OrdinalTypeArray1DHost °rees) const
For hierarchical bases, returns the field ordinals that have at most the specified degree in each dim...
virtual void getValues(BasisValues< OutputValueType, DeviceType > outputValues, const TensorPoints< PointValueType, DeviceType > inputPoints, const EOperator operatorType=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell, using point and output value containers that allow pre...
ScalarView< PointScalar, DeviceType > allocateAndFillExpandedRawPointView() const
This method is for compatibility with existing methods that take raw point views. Note that in genera...
virtual BasisValues< OutputValueType, DeviceType > allocateBasisValues(TensorPoints< PointValueType, DeviceType > points, const EOperator operatorType=OPERATOR_VALUE) const
Allocate BasisValues container suitable for passing to the getValues() variant that takes a TensorPoi...
The data containers in Intrepid2 that support sum factorization and other reduced-data optimizations ...
EBasis
Enumeration of basis types for discrete spaces in Intrepid.
EBasis basisType_
Type of the basis.
const OrdinalTypeArray3DHost getAllDofOrdinal() const
DoF tag to ordinal data structure.
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
Kokkos::View< ordinal_type *, DeviceType > OrdinalTypeArray1D
View type for 1d device array.
PointValueType getDummyPointValue()
Dummy array to receive input arguments.
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
DeviceType DeviceType
(Kokkos) Device type on which Basis is templated. Does not necessarily return true for Kokkos::is_dev...
Contains definitions of custom data types in Intrepid2.
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
ordinal_type getDegree() const
Returns the degree of the basis.
TensorDataType & tensorData()
TensorData accessor for single-family scalar data.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
virtual bool requireOrientation() const
True if orientation is required.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...
ordinal_type getDofOrdinal(const ordinal_type subcDim, const ordinal_type subcOrd, const ordinal_type subcDofOrd) const
DoF tag to ordinal lookup.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
const OrdinalTypeArray2DHost getAllDofTags() const
Retrieves all DoF tags.
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.
shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation https://trili...
virtual const char * getName() const
Returns basis name.
std::vector< int > getPolynomialDegreeOfFieldAsVector(int fieldOrdinal) const
For hierarchical bases, returns the polynomial degree (which may have multiple values in higher spati...
Kokkos::View< ECoordinates, DeviceType > ECoordinatesViewType
View for coordinate system type.
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
Implementation file for the abstract base class Intrepid2::Basis.
const VectorDataType & vectorData() const
VectorData accessor.
EBasis getBasisType() const
Returns the basis type.
virtual BasisPtr< DeviceType, OutputValueType, PointValueType > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const
returns the basis associated to a subCell.
Kokkos::DynRankView< scalarType, DeviceType > dofCoords_
Coordinates of degrees-of-freedom for basis functions defined in physical space.
Kokkos::View< ordinal_type *, Kokkos::LayoutStride, Kokkos::HostSpace > OrdinalTypeArrayStride1DHost
View type for 1d host array.
Kokkos::View< ordinal_type, DeviceType > OrdinalViewType
View type for ordinal.
Header file for the data-wrapper class Intrepid2::BasisValues.
Kokkos::View< ordinal_type **, DeviceType > OrdinalTypeArray2D
View type for 2d device array.
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now...
Reference-space field values for a basis, designed to support typical vector-valued bases...
Kokkos::View< ordinal_type *, Kokkos::LayoutStride, DeviceType > OrdinalTypeArrayStride1D
View type for 1d device array.
int getPolynomialDegreeLength() const
For hierarchical bases, returns the number of entries required to specify the polynomial degree of a ...
EFunctionSpace getFunctionSpace() const
Returns the function space for the basis.
Kokkos::DynRankView< scalarType, DeviceType > dofCoeffs_
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...