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),
38// Mauro Perego (mperego@sandia.gov), or
39// Nate Roberts (nvrober@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50#ifndef __INTREPID2_BASIS_HPP__
51#define __INTREPID2_BASIS_HPP__
52
53#include "Intrepid2_ConfigDefs.hpp"
54#include "Intrepid2_Types.hpp"
55#include "Intrepid2_Utils.hpp"
56
60#include "Kokkos_Vector.hpp"
61#include "Shards_CellTopology.hpp"
62#include <Teuchos_RCPDecl.hpp>
63
64#include <vector>
65
66namespace Intrepid2 {
67
68template<typename DeviceType = void,
69 typename OutputType = double,
70 typename PointType = double>
71class Basis;
72
75template <typename DeviceType = void, typename OutputType = double, typename PointType = double>
76using BasisPtr = Teuchos::RCP<Basis<DeviceType,OutputType,PointType> >;
77
80template <typename OutputType = double, typename PointType = double>
82
119 template<typename Device,
120 typename outputValueType,
121 typename pointValueType>
122 class Basis {
123 public:
126 using DeviceType = Device;
127
130 using ExecutionSpace = typename DeviceType::execution_space;
131
132
135 using OutputValueType = outputValueType;
136
139 using PointValueType = pointValueType;
140
143 using OrdinalViewType = Kokkos::View<ordinal_type,DeviceType>;
144
147 using EBasisViewType = Kokkos::View<EBasis,DeviceType>;
148
151 using ECoordinatesViewType = Kokkos::View<ECoordinates,DeviceType>;
152
153 // ** tag interface
154 // - tag interface is not decorated with Kokkos inline so it should be allocated on hostspace
155
158 using OrdinalTypeArray1DHost = Kokkos::View<ordinal_type*,typename ExecutionSpace::array_layout,Kokkos::HostSpace>;
159
162 using OrdinalTypeArray2DHost = Kokkos::View<ordinal_type**,typename ExecutionSpace::array_layout,Kokkos::HostSpace>;
163
166 using OrdinalTypeArray3DHost = Kokkos::View<ordinal_type***,typename ExecutionSpace::array_layout,Kokkos::HostSpace>;
167
170 using OrdinalTypeArrayStride1DHost = Kokkos::View<ordinal_type*, Kokkos::LayoutStride, Kokkos::HostSpace>;
171
174 using OrdinalTypeArray1D = Kokkos::View<ordinal_type*,DeviceType>;
175
178 using OrdinalTypeArray2D = Kokkos::View<ordinal_type**,DeviceType>;
179
182 using OrdinalTypeArray3D = Kokkos::View<ordinal_type***,DeviceType>;
183
186 using OrdinalTypeArrayStride1D = Kokkos::View<ordinal_type*, Kokkos::LayoutStride, DeviceType>;
187
190 typedef typename ScalarTraits<pointValueType>::scalar_type scalarType;
191 protected:
192
195 ordinal_type basisCardinality_;
196
199 ordinal_type basisDegree_;
200
208 shards::CellTopology basisCellTopology_;
209
213
216 ECoordinates basisCoordinates_;
217
220 EFunctionSpace functionSpace_ = FUNCTION_SPACE_MAX;
221
224 //Kokkos::View<bool,DeviceType> basisTagsAreSet_;
225
238
251
263 template<typename OrdinalTypeView3D,
264 typename OrdinalTypeView2D,
265 typename OrdinalTypeView1D>
266 void setOrdinalTagData( OrdinalTypeView3D &tagToOrdinal,
267 OrdinalTypeView2D &ordinalToTag,
268 const OrdinalTypeView1D tags,
269 const ordinal_type basisCard,
270 const ordinal_type tagSize,
271 const ordinal_type posScDim,
272 const ordinal_type posScOrd,
273 const ordinal_type posDfOrd ) {
274 // Create ordinalToTag
275 ordinalToTag = OrdinalTypeView2D("ordinalToTag", basisCard, tagSize);
276
277 // Initialize with -1
278 Kokkos::deep_copy( ordinalToTag, -1 );
279
280 // Copy tags
281 for (ordinal_type i=0;i<basisCard;++i)
282 for (ordinal_type j=0;j<tagSize;++j)
283 ordinalToTag(i, j) = tags(i*tagSize + j);
284
285 // Find out dimension of tagToOrdinal
286 auto maxScDim = 0; // first dimension of tagToOrdinal
287 for (ordinal_type i=0;i<basisCard;++i)
288 if (maxScDim < tags(i*tagSize + posScDim))
289 maxScDim = tags(i*tagSize + posScDim);
290 ++maxScDim;
291
292 auto maxScOrd = 0; // second dimension of tagToOrdinal
293 for (ordinal_type i=0;i<basisCard;++i)
294 if (maxScOrd < tags(i*tagSize + posScOrd))
295 maxScOrd = tags(i*tagSize + posScOrd);
296 ++maxScOrd;
297
298 auto maxDfOrd = 0; // third dimension of tagToOrdinal
299 for (ordinal_type i=0;i<basisCard;++i)
300 if (maxDfOrd < tags(i*tagSize + posDfOrd))
301 maxDfOrd = tags(i*tagSize + posDfOrd);
302 ++maxDfOrd;
303
304 // Create tagToOrdinal
305 tagToOrdinal = OrdinalTypeView3D("tagToOrdinal", maxScDim, maxScOrd, maxDfOrd);
306
307 // Initialize with -1
308 Kokkos::deep_copy( tagToOrdinal, -1 );
309
310 // Overwrite elements of the array corresponding to tags with local DoF Id's, leave all other = -1
311 for (ordinal_type i=0;i<basisCard;++i)
312 tagToOrdinal(tags(i*tagSize), tags(i*tagSize+1), tags(i*tagSize+2)) = i;
313 }
314
315 // dof coords
318 Kokkos::DynRankView<scalarType,DeviceType> dofCoords_;
319
320 // dof coeffs
328 Kokkos::DynRankView<scalarType,DeviceType> dofCoeffs_;
329
338
351 public:
352
353 Basis() = default;
354 virtual~Basis() = default;
355
356 // receives input arguments
359 OutputValueType getDummyOutputValue() { return outputValueType(); }
360
363 PointValueType getDummyPointValue() { return pointValueType(); }
364
367 using OutputViewType = Kokkos::DynRankView<OutputValueType,Kokkos::LayoutStride,DeviceType>;
368
371 using PointViewType = Kokkos::DynRankView<PointValueType,Kokkos::LayoutStride,DeviceType>;
372
375 using ScalarViewType = Kokkos::DynRankView<scalarType,Kokkos::LayoutStride,DeviceType>;
376
381 Kokkos::DynRankView<OutputValueType,DeviceType> allocateOutputView( const int numPoints, const EOperator operatorType = OPERATOR_VALUE) const;
382
388 virtual BasisValues<OutputValueType,DeviceType> allocateBasisValues( TensorPoints<PointValueType,DeviceType> points, const EOperator operatorType = OPERATOR_VALUE) const
389 {
390 const bool operatorIsDk = (operatorType >= OPERATOR_D1) && (operatorType <= OPERATOR_D10);
391 const bool operatorSupported = (operatorType == OPERATOR_VALUE) || (operatorType == OPERATOR_GRAD) || (operatorType == OPERATOR_CURL) || (operatorType == OPERATOR_DIV) || operatorIsDk;
392 INTREPID2_TEST_FOR_EXCEPTION(!operatorSupported, std::invalid_argument, "operator is not supported by allocateBasisValues");
393
394 const int numPoints = points.extent_int(0);
395
396 using Scalar = OutputValueType;
397
398 auto dataView = allocateOutputView(numPoints, operatorType);
399 Data<Scalar,DeviceType> data(dataView);
400
401 bool useVectorData = (dataView.rank() == 3);
402
403 if (useVectorData)
404 {
405 VectorData<Scalar,DeviceType> vectorData(data);
406 return BasisValues<Scalar,DeviceType>(vectorData);
407 }
408 else
409 {
410 TensorData<Scalar,DeviceType> tensorData(data);
411 return BasisValues<Scalar,DeviceType>(tensorData);
412 }
413 }
414
433 virtual
434 void
435 getValues( OutputViewType /* outputValues */,
436 const PointViewType /* inputPoints */,
437 const EOperator /* operatorType */ = OPERATOR_VALUE ) const {
438 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
439 ">>> ERROR (Basis::getValues): this method (FEM) is not supported or should be overridden accordingly by derived classes.");
440 }
441
453 virtual
454 void
457 const EOperator operatorType = OPERATOR_VALUE ) const {
458 // note the extra allocation/copy here (this is one reason, among several, to override this method):
459 auto rawExpandedPoints = inputPoints.allocateAndFillExpandedRawPointView();
460
461 OutputViewType rawOutputView;
463 if (outputValues.numTensorDataFamilies() > 0)
464 {
465 INTREPID2_TEST_FOR_EXCEPTION(outputValues.tensorData(0).numTensorComponents() != 1, std::invalid_argument, "default implementation of getValues() only supports outputValues with trivial tensor-product structure");
466 outputData = outputValues.tensorData().getTensorComponent(0);
467 }
468 else if (outputValues.vectorData().isValid())
469 {
470 INTREPID2_TEST_FOR_EXCEPTION(outputValues.vectorData().numComponents() != 1, std::invalid_argument, "default implementation of getValues() only supports outputValues with trivial tensor-product structure");
471 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");
472 outputData = outputValues.vectorData().getComponent(0).getTensorComponent(0);
473 }
474
475 this->getValues(outputData.getUnderlyingView(), rawExpandedPoints, operatorType);
476 }
477
497 virtual
498 void
499 getValues( OutputViewType /* outputValues */,
500 const PointViewType /* inputPoints */,
501 const PointViewType /* cellVertices */,
502 const EOperator /* operatorType */ = OPERATOR_VALUE ) const {
503 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
504 ">>> ERROR (Basis::getValues): this method (FVM) is not supported or should be overridden accordingly by derived classes.");
505 }
506
507
511 virtual
512 void
513 getDofCoords( ScalarViewType /* dofCoords */ ) const {
514 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
515 ">>> ERROR (Basis::getDofCoords): this method is not supported or should be overridden accordingly by derived classes.");
516 }
517
526 virtual
527 void
528 getDofCoeffs( ScalarViewType /* dofCoeffs */ ) const {
529 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
530 ">>> ERROR (Basis::getDofCoeffs): this method is not supported or should be overridden accordingly by derived classes.");
531 }
532
541 {
542 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
543 ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
544 int degreeEntryLength = fieldOrdinalPolynomialDegree_.extent_int(1);
545 int requestedDegreeLength = degrees.extent_int(0);
546 INTREPID2_TEST_FOR_EXCEPTION(degreeEntryLength != requestedDegreeLength, std::invalid_argument, "length of degrees does not match the entries in fieldOrdinalPolynomialDegree_");
547 std::vector<int> fieldOrdinalsVector;
548 for (int basisOrdinal=0; basisOrdinal<fieldOrdinalPolynomialDegree_.extent_int(0); basisOrdinal++)
549 {
550 bool matches = true;
551 for (int d=0; d<degreeEntryLength; d++)
552 {
553 if (fieldOrdinalPolynomialDegree_(basisOrdinal,d) > degrees(d)) matches = false;
554 }
555 if (matches) fieldOrdinalsVector.push_back(basisOrdinal);
556 }
557 OrdinalTypeArray1DHost fieldOrdinals("fieldOrdinalsForDegree",fieldOrdinalsVector.size());
558 for (unsigned i=0; i<fieldOrdinalsVector.size(); i++)
559 {
560 fieldOrdinals(i) = fieldOrdinalsVector[i];
561 }
562 return fieldOrdinals;
563 }
564
573 {
574 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
575 ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
576 int degreeEntryLength = fieldOrdinalH1PolynomialDegree_.extent_int(1);
577 int requestedDegreeLength = degrees.extent_int(0);
578 INTREPID2_TEST_FOR_EXCEPTION(degreeEntryLength != requestedDegreeLength, std::invalid_argument, "length of degrees does not match the entries in fieldOrdinalPolynomialDegree_");
579 std::vector<int> fieldOrdinalsVector;
580 for (int basisOrdinal=0; basisOrdinal<fieldOrdinalH1PolynomialDegree_.extent_int(0); basisOrdinal++)
581 {
582 bool matches = true;
583 for (int d=0; d<degreeEntryLength; d++)
584 {
585 if (fieldOrdinalH1PolynomialDegree_(basisOrdinal,d) > degrees(d)) matches = false;
586 }
587 if (matches) fieldOrdinalsVector.push_back(basisOrdinal);
588 }
589 OrdinalTypeArray1DHost fieldOrdinals("fieldOrdinalsForH1Degree",fieldOrdinalsVector.size());
590 for (unsigned i=0; i<fieldOrdinalsVector.size(); i++)
591 {
592 fieldOrdinals(i) = fieldOrdinalsVector[i];
593 }
594 return fieldOrdinals;
595 }
596
608 std::vector<int> getFieldOrdinalsForDegree(std::vector<int> &degrees) const
609 {
610 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
611 ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
612 OrdinalTypeArray1DHost degreesView("degrees",degrees.size());
613 for (unsigned d=0; d<degrees.size(); d++)
614 {
615 degreesView(d) = degrees[d];
616 }
617 auto fieldOrdinalsView = getFieldOrdinalsForDegree(degreesView);
618 std::vector<int> fieldOrdinalsVector(fieldOrdinalsView.extent_int(0));
619 for (int i=0; i<fieldOrdinalsView.extent_int(0); i++)
620 {
621 fieldOrdinalsVector[i] = fieldOrdinalsView(i);
622 }
623 return fieldOrdinalsVector;
624 }
625
636 std::vector<int> getFieldOrdinalsForH1Degree(std::vector<int> &degrees) const
637 {
638 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
639 ">>> ERROR (Basis::getFieldOrdinalsForDegree): this method is not supported for non-hierarchical bases.");
640 OrdinalTypeArray1DHost degreesView("degrees",degrees.size());
641 for (unsigned d=0; d<degrees.size(); d++)
642 {
643 degreesView(d) = degrees[d];
644 }
645 auto fieldOrdinalsView = getFieldOrdinalsForH1Degree(degreesView);
646 std::vector<int> fieldOrdinalsVector(fieldOrdinalsView.extent_int(0));
647 for (int i=0; i<fieldOrdinalsView.extent_int(0); i++)
648 {
649 fieldOrdinalsVector[i] = fieldOrdinalsView(i);
650 }
651 return fieldOrdinalsVector;
652 }
653
661 {
662 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
663 ">>> ERROR (Basis::getPolynomialDegreeOfField): this method is not supported for non-hierarchical bases.");
664 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal < 0, std::invalid_argument, "field ordinal must be non-negative");
665 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal >= fieldOrdinalPolynomialDegree_.extent_int(0), std::invalid_argument, "field ordinal out of bounds");
666
667 int polyDegreeLength = getPolynomialDegreeLength();
668 OrdinalTypeArray1DHost polyDegree("polynomial degree", polyDegreeLength);
669 for (int d=0; d<polyDegreeLength; d++)
670 {
671 polyDegree(d) = fieldOrdinalPolynomialDegree_(fieldOrdinal,d);
672 }
673 return polyDegree;
674 }
675
683 {
684 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
685 ">>> ERROR (Basis::getPolynomialDegreeOfField): this method is not supported for non-hierarchical bases.");
686 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal < 0, std::invalid_argument, "field ordinal must be non-negative");
687 INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinal >= fieldOrdinalH1PolynomialDegree_.extent_int(0), std::invalid_argument, "field ordinal out of bounds");
688
689 int polyDegreeLength = getPolynomialDegreeLength();
690 OrdinalTypeArray1DHost polyDegree("polynomial degree", polyDegreeLength);
691 for (int d=0; d<polyDegreeLength; d++)
692 {
693 polyDegree(d) = fieldOrdinalH1PolynomialDegree_(fieldOrdinal,d);
694 }
695 return polyDegree;
696 }
697
705 std::vector<int> getPolynomialDegreeOfFieldAsVector(int fieldOrdinal) const
706 {
707 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
708 ">>> ERROR (Basis::getPolynomialDegreeOfFieldAsVector): this method is not supported for non-hierarchical bases.");
709 auto polynomialDegreeView = getPolynomialDegreeOfField(fieldOrdinal);
710 std::vector<int> polynomialDegree(polynomialDegreeView.extent_int(0));
711
712 for (unsigned d=0; d<polynomialDegree.size(); d++)
713 {
714 polynomialDegree[d] = polynomialDegreeView(d);
715 }
716 return polynomialDegree;
717 }
718
726 std::vector<int> getH1PolynomialDegreeOfFieldAsVector(int fieldOrdinal) const
727 {
728 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
729 ">>> ERROR (Basis::getPolynomialDegreeOfFieldAsVector): this method is not supported for non-hierarchical bases.");
730 auto polynomialDegreeView = getH1PolynomialDegreeOfField(fieldOrdinal);
731 std::vector<int> polynomialDegree(polynomialDegreeView.extent_int(0));
732
733 for (unsigned d=0; d<polynomialDegree.size(); d++)
734 {
735 polynomialDegree[d] = polynomialDegreeView(d);
736 }
737 return polynomialDegree;
738 }
739
743 {
744 INTREPID2_TEST_FOR_EXCEPTION( basisType_ != BASIS_FEM_HIERARCHICAL, std::logic_error,
745 ">>> ERROR (Basis::getPolynomialDegreeLength): this method is not supported for non-hierarchical bases.");
746 return fieldOrdinalPolynomialDegree_.extent_int(1);
747 }
748
753 virtual
754 const char*
755 getName() const {
756 return "Intrepid2_Basis";
757 }
758
761 virtual
762 bool
764 return false;
765 }
766
771 ordinal_type
773 return basisCardinality_;
774 }
775
776
781 ordinal_type
782 getDegree() const {
783 return basisDegree_;
784 }
785
790 EFunctionSpace
792 return functionSpace_;
793 }
794
800 shards::CellTopology
802 return basisCellTopology_;
803 }
804
805
810 EBasis
811 getBasisType() const {
812 return basisType_;
813 }
814
815
820 ECoordinates
822 return basisCoordinates_;
823 }
824
832 ordinal_type
833 getDofCount( const ordinal_type subcDim,
834 const ordinal_type subcOrd ) const {
835 if ( subcDim >= 0 && subcDim < static_cast<ordinal_type>(tagToOrdinal_.extent(0)) &&
836 subcOrd >= 0 && subcOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(1)) )
837 {
838 int firstDofOrdinal = tagToOrdinal_(subcDim, subcOrd, 0); // will be -1 if no dofs for subcell
839 if (firstDofOrdinal == -1) return static_cast<ordinal_type>(0);
840 // otherwise, lookup its tag and return the dof count stored there
841 return static_cast<ordinal_type>(this->getDofTag(firstDofOrdinal)[3]);
842 }
843 else
844 {
845 // subcDim and/or subcOrd out of bounds -> no dofs associated with subcell
846 return static_cast<ordinal_type>(0);
847 }
848 }
849
858 ordinal_type
859 getDofOrdinal( const ordinal_type subcDim,
860 const ordinal_type subcOrd,
861 const ordinal_type subcDofOrd ) const {
862 // this should be abort and able to be called as a device function
863#ifdef HAVE_INTREPID2_DEBUG
864 INTREPID2_TEST_FOR_EXCEPTION( subcDim < 0 || subcDim >= static_cast<ordinal_type>(tagToOrdinal_.extent(0)), std::out_of_range,
865 ">>> ERROR (Basis::getDofOrdinal): subcDim is out of range");
866 INTREPID2_TEST_FOR_EXCEPTION( subcOrd < 0 || subcOrd >= static_cast<ordinal_type>(tagToOrdinal_.extent(1)), std::out_of_range,
867 ">>> ERROR (Basis::getDofOrdinal): subcOrd is out of range");
868 INTREPID2_TEST_FOR_EXCEPTION( subcDofOrd < 0 || subcDofOrd >= static_cast<ordinal_type>(tagToOrdinal_.extent(2)), std::out_of_range,
869 ">>> ERROR (Basis::getDofOrdinal): subcDofOrd is out of range");
870#endif
871 ordinal_type r_val = -1;
872 if ( subcDim < static_cast<ordinal_type>(tagToOrdinal_.extent(0)) &&
873 subcOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(1)) &&
874 subcDofOrd < static_cast<ordinal_type>(tagToOrdinal_.extent(2)) )
875 r_val = tagToOrdinal_(subcDim, subcOrd, subcDofOrd);
876#ifdef HAVE_INTREPID2_DEBUG
877 INTREPID2_TEST_FOR_EXCEPTION( r_val == -1, std::runtime_error,
878 ">>> ERROR (Basis::getDofOrdinal): Invalid DoF tag is found.");
879#endif
880 return r_val;
881 }
882
885 virtual int getNumTensorialExtrusions() const
886 {
887 return 0;
888 }
889
893 return tagToOrdinal_;
894 }
895
896
908 getDofTag( const ordinal_type dofOrd ) const {
909#ifdef HAVE_INTREPID2_DEBUG
910 INTREPID2_TEST_FOR_EXCEPTION( dofOrd < 0 || dofOrd >= static_cast<ordinal_type>(ordinalToTag_.extent(0)), std::out_of_range,
911 ">>> ERROR (Basis::getDofTag): dofOrd is out of range");
912#endif
913 return Kokkos::subview(ordinalToTag_, dofOrd, Kokkos::ALL());
914 }
915
916
927 return ordinalToTag_;
928 }
929
945 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const {
946 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
947 ">>> ERROR (Basis::getSubCellRefBasis): this method is not supported or should be overridden accordingly by derived classes.");
948 }
949
954 ordinal_type getDomainDimension() const
955 {
956 return this->getBaseCellTopology().getDimension() + this->getNumTensorialExtrusions();
957 }
958
964 getHostBasis() const {
965 INTREPID2_TEST_FOR_EXCEPTION( true, std::logic_error,
966 ">>> ERROR (Basis::getHostBasis): this method is not supported or should be overridden accordingly by derived classes.");
967 }
968
969 }; // class Basis
970
971 //--------------------------------------------------------------------------------------------//
972 // //
973 // Helper functions of the Basis class //
974 // //
975 //--------------------------------------------------------------------------------------------//
976
977 //
978 // functions for orders, cardinality, etc.
979 //
980
981
993 KOKKOS_INLINE_FUNCTION
994 ordinal_type getFieldRank(const EFunctionSpace spaceType);
995
1031 KOKKOS_INLINE_FUNCTION
1032 ordinal_type getOperatorRank(const EFunctionSpace spaceType,
1033 const EOperator operatorType,
1034 const ordinal_type spaceDim);
1035
1041 KOKKOS_INLINE_FUNCTION
1042 ordinal_type getOperatorOrder(const EOperator operatorType);
1043
1044 template<EOperator operatorType>
1045 KOKKOS_INLINE_FUNCTION
1046 constexpr ordinal_type getOperatorOrder();
1047
1071 template<ordinal_type spaceDim>
1072 KOKKOS_INLINE_FUNCTION
1073 ordinal_type getDkEnumeration(const ordinal_type xMult,
1074 const ordinal_type yMult = -1,
1075 const ordinal_type zMult = -1);
1076
1077
1088 template<ordinal_type spaceDim>
1089 KOKKOS_INLINE_FUNCTION
1090 ordinal_type getPnEnumeration(const ordinal_type p,
1091 const ordinal_type q = 0,
1092 const ordinal_type r = 0);
1093
1094
1095
1114template<typename value_type>
1115KOKKOS_INLINE_FUNCTION
1116void getJacobyRecurrenceCoeffs (
1117 value_type &an,
1118 value_type &bn,
1119 value_type &cn,
1120 const ordinal_type alpha,
1121 const ordinal_type beta ,
1122 const ordinal_type n);
1123
1124
1125
1126
1127
1128 // /** \brief Returns multiplicities of dx, dy, and dz based on the enumeration of the partial
1129 // derivative, its order and the space dimension. Inverse of the getDkEnumeration() method.
1130
1131 // \param partialMult [out] - array with the multiplicities f dx, dy and dz
1132 // \param derivativeEnum [in] - enumeration of the partial derivative
1133 // \param operatorType [in] - k-th partial derivative Dk
1134 // \param spaceDim [in] - space dimension
1135 // */
1136 // template<typename OrdinalArrayType>
1137 // KOKKOS_INLINE_FUNCTION
1138 // void getDkMultiplicities(OrdinalArrayType partialMult,
1139 // const ordinal_type derivativeEnum,
1140 // const EOperator operatorType,
1141 // const ordinal_type spaceDim);
1142
1161 KOKKOS_INLINE_FUNCTION
1162 ordinal_type getDkCardinality(const EOperator operatorType,
1163 const ordinal_type spaceDim);
1164
1165 template<EOperator operatorType, ordinal_type spaceDim>
1166 KOKKOS_INLINE_FUNCTION
1167 constexpr ordinal_type getDkCardinality();
1168
1169
1170
1180 template<ordinal_type spaceDim>
1181 KOKKOS_INLINE_FUNCTION
1182 ordinal_type getPnCardinality (ordinal_type n);
1183
1184 template<ordinal_type spaceDim, ordinal_type n>
1185 KOKKOS_INLINE_FUNCTION
1186 constexpr ordinal_type getPnCardinality ();
1187
1188
1189
1190 //
1191 // Argument check
1192 //
1193
1194
1205 template<typename outputValueViewType,
1206 typename inputPointViewType>
1207 void getValues_HGRAD_Args( const outputValueViewType outputValues,
1208 const inputPointViewType inputPoints,
1209 const EOperator operatorType,
1210 const shards::CellTopology cellTopo,
1211 const ordinal_type basisCard );
1212
1223 template<typename outputValueViewType,
1224 typename inputPointViewType>
1225 void getValues_HCURL_Args( const outputValueViewType outputValues,
1226 const inputPointViewType inputPoints,
1227 const EOperator operatorType,
1228 const shards::CellTopology cellTopo,
1229 const ordinal_type basisCard );
1230
1241 template<typename outputValueViewType,
1242 typename inputPointViewType>
1243 void getValues_HDIV_Args( const outputValueViewType outputValues,
1244 const inputPointViewType inputPoints,
1245 const EOperator operatorType,
1246 const shards::CellTopology cellTopo,
1247 const ordinal_type basisCard );
1248
1259 template<typename outputValueViewType,
1260 typename inputPointViewType>
1261 void getValues_HVOL_Args( const outputValueViewType outputValues,
1262 const inputPointViewType inputPoints,
1263 const EOperator operatorType,
1264 const shards::CellTopology cellTopo,
1265 const ordinal_type basisCard );
1266
1267}// namespace Intrepid2
1268
1269#include <Intrepid2_BasisDef.hpp>
1270
1271//--------------------------------------------------------------------------------------------//
1272// //
1273// D O C U M E N T A T I O N P A G E S //
1274// //
1275//--------------------------------------------------------------------------------------------//
1377#endif
Implementation file for the abstract base class Intrepid2::Basis.
Header file for the data-wrapper class Intrepid2::BasisValues.
BasisPtr< typename Kokkos::HostSpace::device_type, OutputType, PointType > HostBasisPtr
Pointer to a Basis whose device type is on the host (Kokkos::HostSpace::device_type),...
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
Definition of cell topology information.
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
Contains definitions of custom data types in Intrepid2.
Header function for Intrepid2::Util class and other utility functions.
The data containers in Intrepid2 that support sum factorization and other reduced-data optimizations ...
const VectorDataType & vectorData() const
VectorData accessor.
TensorDataType & tensorData()
TensorData accessor for single-family scalar data.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
OrdinalTypeArray1DHost getFieldOrdinalsForDegree(OrdinalTypeArray1DHost &degrees) const
For hierarchical bases, returns the field ordinals that have at most the specified degree in each dim...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
OrdinalTypeArray1DHost getFieldOrdinalsForH1Degree(OrdinalTypeArray1DHost &degrees) const
For hierarchical bases, returns the field ordinals that have at most the specified H^1 degree in each...
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
const OrdinalTypeArrayStride1DHost getDofTag(const ordinal_type dofOrd) const
DoF ordinal to DoF tag lookup.
outputValueType OutputValueType
Output value type for basis; default is double.
ordinal_type getDofOrdinal(const ordinal_type subcDim, const ordinal_type subcOrd, const ordinal_type subcDofOrd) const
DoF tag to ordinal lookup.
ordinal_type getDofCount(const ordinal_type subcDim, const ordinal_type subcOrd) const
DoF count for specified subcell.
Kokkos::View< ordinal_type ***, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray3DHost
View type for 3d host array.
const OrdinalTypeArray2DHost getAllDofTags() const
Retrieves all DoF tags.
std::vector< int > getFieldOrdinalsForH1Degree(std::vector< int > &degrees) const
For hierarchical bases, returns the field ordinals that have at most the specified H^1 degree in each...
std::vector< int > getH1PolynomialDegreeOfFieldAsVector(int fieldOrdinal) const
For hierarchical bases, returns the polynomial degree (which may have multiple values in higher spati...
EBasis basisType_
Type of the basis.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
Kokkos::View< ordinal_type *, Kokkos::LayoutStride, DeviceType > OrdinalTypeArrayStride1D
View type for 1d device array.
Kokkos::View< ordinal_type *, Kokkos::LayoutStride, Kokkos::HostSpace > OrdinalTypeArrayStride1DHost
View type for 1d host array.
ordinal_type getDegree() const
Returns the degree of the basis.
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.
virtual HostBasisPtr< OutputValueType, PointValueType > getHostBasis() const
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
ordinal_type getCardinality() const
Returns cardinality of the basis.
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...
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
OutputValueType getDummyOutputValue()
Dummy array to receive input arguments.
ECoordinates getCoordinateSystem() const
Returns the type of coordinate system for which the basis is defined.
OrdinalTypeArray2DHost ordinalToTag_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
PointValueType getDummyPointValue()
Dummy array to receive input arguments.
int getPolynomialDegreeLength() const
For hierarchical bases, returns the number of entries required to specify the polynomial degree of a ...
virtual bool requireOrientation() const
True if orientation is required.
Kokkos::View< ordinal_type, DeviceType > OrdinalViewType
View type for ordinal.
const OrdinalTypeArray3DHost getAllDofOrdinal() const
DoF tag to ordinal data structure.
Device DeviceType
(Kokkos) Device type on which Basis is templated. Does not necessarily return true for Kokkos::is_dev...
Kokkos::DynRankView< scalarType, DeviceType > dofCoords_
Coordinates of degrees-of-freedom for basis functions defined in physical space.
Kokkos::View< ECoordinates, DeviceType > ECoordinatesViewType
View for coordinate system type.
ordinal_type getDomainDimension() const
Returns the spatial dimension of the domain of the basis; this is equal to getBaseCellTopology()....
OrdinalTypeArray2DHost fieldOrdinalH1PolynomialDegree_
H^1 polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
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...
pointValueType PointValueType
Point value type for basis; default is double.
Kokkos::View< ordinal_type **, DeviceType > OrdinalTypeArray2D
View type for 2d device array.
Kokkos::DynRankView< scalarType, DeviceType > dofCoeffs_
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...
OrdinalTypeArray1DHost getPolynomialDegreeOfField(int fieldOrdinal) const
For hierarchical bases, returns the polynomial degree (which may have multiple values in higher spati...
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
virtual void getDofCoords(ScalarViewType) const
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
ordinal_type basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
OrdinalTypeArray3DHost tagToOrdinal_
DoF tag to ordinal lookup table.
EFunctionSpace getFunctionSpace() const
Returns the function space for the basis.
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...
virtual void getValues(OutputViewType, const PointViewType, const EOperator=OPERATOR_VALUE) const
Evaluation of a FEM basis on a reference cell.
virtual void getDofCoeffs(ScalarViewType) const
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...
EBasis getBasisType() const
Returns the basis type.
OrdinalTypeArray1DHost getH1PolynomialDegreeOfField(int fieldOrdinal) const
For hierarchical bases, returns the polynomial degree (which may have multiple values in higher spati...
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
Kokkos::View< ordinal_type ***, DeviceType > OrdinalTypeArray3D
View type for 3d device array.
shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation https://trili...
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package for definition of b...
Kokkos::View< ordinal_type *, DeviceType > OrdinalTypeArray1D
View type for 1d 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.
typename DeviceType::execution_space ExecutionSpace
(Kokkos) Execution space for basis.
virtual const char * getName() const
Returns basis name.
std::vector< int > getFieldOrdinalsForDegree(std::vector< int > &degrees) const
For hierarchical bases, returns the field ordinals that have at most the specified degree in each dim...
OrdinalTypeArray2DHost fieldOrdinalPolynomialDegree_
Polynomial degree for each degree of freedom. Only defined for hierarchical bases right now....
virtual BasisPtr< DeviceType, OutputValueType, PointValueType > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const
returns the basis associated to a subCell.
virtual int getNumTensorialExtrusions() const
returns the number of tensorial extrusions relative to the cell topology returned by getBaseCellTopol...
EFunctionSpace functionSpace_
The function space in which the basis is defined.
std::vector< int > getPolynomialDegreeOfFieldAsVector(int fieldOrdinal) const
For hierarchical bases, returns the polynomial degree (which may have multiple values in higher spati...
Kokkos::View< EBasis, DeviceType > EBasisViewType
View for basis type.
Wrapper around a Kokkos::View that allows data that is constant or repeating in various logical dimen...
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.
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
View-like interface to tensor data; tensor components are stored separately and multiplied together a...
KOKKOS_INLINE_FUNCTION const Data< Scalar, DeviceType > & getTensorComponent(const ordinal_type &r) const
Returns the requested tensor component.
KOKKOS_INLINE_FUNCTION ordinal_type numTensorComponents() const
Return the number of tensorial components.
View-like interface to tensor points; point components are stored separately; the appropriate coordin...
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.
ScalarView< PointScalar, DeviceType > allocateAndFillExpandedRawPointView() const
This method is for compatibility with existing methods that take raw point views. Note that in genera...
Reference-space field values for a basis, designed to support typical vector-valued bases.
KOKKOS_INLINE_FUNCTION int numComponents() const
returns the number of components
KOKKOS_INLINE_FUNCTION const TensorData< Scalar, DeviceType > & getComponent(const int &componentOrdinal) const
Single-argument component accessor for the axial-component or the single-family case; in this case,...
KOKKOS_INLINE_FUNCTION constexpr bool isValid() const
returns true for containers that have data; false for those that don't (e.g., those that have been co...