Skip to content

Commit

Permalink
Merge 'trilinos/Trilinos:develop' (c3b1a94) into 'tcad-charon/Trilino…
Browse files Browse the repository at this point in the history 
…s:develop' (8fc5840).

* trilinos-develop:
  Intrepid2: Added shards-compatible serendipity HGRAD elements and enable projections with wedges (trilinos#11928)
  Use space object in m_impl_handle.create_functor calls
  sacado, stokhos: add exec space instances to allocate_shared
  • Loading branch information
Charonops Jenkins Pipeline committed Jun 3, 2023
2 parents 8fc5840 + c3b1a94 commit ab899a0
Show file tree
Hide file tree
Showing 37 changed files with 3,417 additions and 684 deletions.
48 changes: 44 additions & 4 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HCURL_WEDGE.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,8 +229,8 @@ namespace Intrepid2
Note that getDofCoeffs() is supported only for Lagrangian bases.
*/
virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),std::make_pair(0,2));
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),2);
this->TensorBasis::getDofCoeffs(dofCoeffs1);
Kokkos::deep_copy(dofCoeffs2,0.0);
}
Expand Down Expand Up @@ -376,8 +376,8 @@ namespace Intrepid2
Note that getDofCoeffs() is supported only for Lagrangian bases.
*/
virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),std::make_pair(0,2));
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),2);
Kokkos::deep_copy(dofCoeffs1,0.0);
this->TensorBasis::getDofCoeffs(dofCoeffs2);
}
Expand Down Expand Up @@ -448,6 +448,46 @@ namespace Intrepid2
return name_.c_str();
}

/** \brief returns the basis associated to a subCell.
The bases of the subCell are the restriction to the subCell
of the bases of the parent cell.
TODO: test this method when different orders are used in different directions
\param [in] subCellDim - dimension of subCell
\param [in] subCellOrd - position of the subCell among of the subCells having the same dimension
\return pointer to the subCell basis of dimension subCellDim and position subCellOrd
*/
Teuchos::RCP<BasisBase>
getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
{
using LineBasis = HVOL_LINE;
using TriCurlBasis = HCURL_TRI;
using QuadCurlBasis = Basis_Derived_HCURL_QUAD<HGRAD_LINE,HVOL_LINE>;
if(subCellDim == 1) {
if(subCellOrd < 6) //edges associated to basal and upper triagular faces
return Teuchos::rcp( new LineBasis(order_xy_-1, pointType_) );
else
return Teuchos::rcp( new LineBasis(order_z_-1, pointType_) );
}
else if(subCellDim == 2) {
switch(subCellOrd) {
case 0:
return Teuchos::rcp( new QuadCurlBasis(order_xy_, order_z_, pointType_) );
case 1:
return Teuchos::rcp( new QuadCurlBasis(order_xy_, order_z_, pointType_) );
case 2:
return Teuchos::rcp( new QuadCurlBasis(order_z_, order_xy_, pointType_) );
case 3:
return Teuchos::rcp( new TriCurlBasis(order_xy_, pointType_) );
case 4:
return Teuchos::rcp( new TriCurlBasis(order_xy_, pointType_) );
default:
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellOrd is out of bounds");
}
}
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellDim is out of bounds");
}

/** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this.
\return Pointer to the new Basis object.
Expand Down
51 changes: 43 additions & 8 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HDIV_WEDGE.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,10 +197,10 @@ namespace Intrepid2
Note that getDofCoeffs() is supported only for Lagrangian bases.
*/
virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
Kokkos::deep_copy(dofCoeffs1,0.0);
this->TensorBasis::getDofCoeffs(dofCoeffs2);
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), std::make_pair(0,2));
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), 2);
this->TensorBasis::getDofCoeffs(dofCoeffs1);
Kokkos::deep_copy(dofCoeffs2,0.0);
}
};

Expand Down Expand Up @@ -329,10 +329,10 @@ namespace Intrepid2
Note that getDofCoeffs() is supported only for Lagrangian bases.
*/
virtual void getDofCoeffs( ScalarViewType dofCoeffs ) const override {
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),0);
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(),1);
this->TensorBasis::getDofCoeffs(dofCoeffs1);
Kokkos::deep_copy(dofCoeffs2,0.0);
auto dofCoeffs1 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), std::make_pair(0,2));
auto dofCoeffs2 = Kokkos::subview(dofCoeffs,Kokkos::ALL(), 2);
Kokkos::deep_copy(dofCoeffs1,0.0);
this->TensorBasis::getDofCoeffs(dofCoeffs2);
}
};

Expand Down Expand Up @@ -402,6 +402,41 @@ namespace Intrepid2
return name_.c_str();
}


/** \brief returns the basis associated to a subCell.
The bases of the subCell are the restriction to the subCell
of the bases of the parent cell.
TODO: test this method when different orders are used in different directions
\param [in] subCellDim - dimension of subCell
\param [in] subCellOrd - position of the subCell among of the subCells having the same dimension
\return pointer to the subCell basis of dimension subCellDim and position subCellOrd
*/
Teuchos::RCP<BasisBase>
getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
{
using QuadBasis = Basis_Derived_HVOL_QUAD<HVOL_LINE>;
using TriBasis = HVOL_TRI;

if(subCellDim == 2) {
switch(subCellOrd) {
case 0:
return Teuchos::rcp( new QuadBasis(order_xy_-1, order_z_-1, pointType_) );
case 1:
return Teuchos::rcp( new QuadBasis(order_xy_-1, order_z_-1, pointType_) );
case 2:
return Teuchos::rcp( new QuadBasis(order_z_-1, order_xy_-1, pointType_) );
case 3:
return Teuchos::rcp( new TriBasis(order_xy_-1, pointType_) );
case 4:
return Teuchos::rcp( new TriBasis(order_xy_-1, pointType_) );
default:
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellOrd is out of bounds");
}
}
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellDim is out of bounds");
}

/** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this.
\return Pointer to the new Basis object.
Expand Down
41 changes: 22 additions & 19 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_DerivedBasis_HGRAD_WEDGE.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@
#define Intrepid2_DerivedBasis_HGRAD_WEDGE_h

#include "Intrepid2_TensorBasis.hpp"
#include "Intrepid2_DerivedBasis_HGRAD_QUAD.hpp"

namespace Intrepid2
{
Expand Down Expand Up @@ -137,27 +138,29 @@ namespace Intrepid2
Teuchos::RCP<BasisBase>
getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
{
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Method not yet implemented");
if(subCellDim == 2) {
if(subCellDim == 1) {
if(subCellOrd < 6) //edges associated to basal and upper triagular faces
return Teuchos::rcp( new LineBasis(order_xy_, pointType_) );
else
return Teuchos::rcp( new LineBasis(order_z_, pointType_) );
}
else if(subCellDim == 2) {
switch(subCellOrd) {
// TODO: check this subcell numbering -- which are the tris, and which are the quads?
// (one way: use the shards CellTopology to look this up…)
case 0:
case 1:
return Teuchos::rcp( new TriBasis(order_xy_, pointType_) );
case 2:
case 3:
case 4:
// TODO: check what order the poly order arguments should go in...
// (one way: use the shards CellTopology to look this up…)
// TODO: define QuadBasis somehow, somewhere. (Probably use DerivedBasis_HGRAD_QUAD.)
// return Teuchos::rcp( new QuadBasis(order_xy_, order_z_, pointType_) );
return Teuchos::null;
case 0:
return Teuchos::rcp( new Intrepid2::Basis_Derived_HGRAD_QUAD<LineBasis>(order_xy_, order_z_, pointType_) );
case 1:
return Teuchos::rcp( new Intrepid2::Basis_Derived_HGRAD_QUAD<LineBasis>(order_xy_, order_z_, pointType_) );
case 2:
return Teuchos::rcp( new Intrepid2::Basis_Derived_HGRAD_QUAD<LineBasis>(order_z_, order_xy_, pointType_) );
case 3:
return Teuchos::rcp( new TriBasis(order_xy_, pointType_) );
case 4:
return Teuchos::rcp( new TriBasis(order_xy_, pointType_) );
default:
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellOrd is out of bounds");
}
return Teuchos::null;
}

INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
}
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"subCellDim is out of bounds");
}

/** \brief Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_type, but is otherwise identical to this.
Expand Down
80 changes: 53 additions & 27 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HGRAD_HEX_C2_FEM.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,12 +54,20 @@

namespace Intrepid2 {

/** \class Intrepid2::Basis_HGRAD_HEX_C2_FEM
/** \class Intrepid2::Basis_HGRAD_HEX_DEG2_FEM
\brief Implementation of the default H(grad)-compatible FEM basis of degree 2 on Hexahedron cell
Implements Lagrangian basis of degree 2 on the reference Hexahedron cell. The basis has
cardinality 27 and spans a COMPLETE tri-quadratic polynomial space. Basis functions are dual
to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows:
Implements Lagrangian basis of degree 2 on the reference Hexahedron cell.
When the serendipity template argument is false, the basis has
cardinality 27 and spans a COMPLETE tri-quadratic polynomial space.
Note, Basis_HGRAD_HEX_C2_FEM = Basis_HGRAD_HEX_DEG2_FEM<false>
When the serendipity template argument is true, the basis has
cardinality 20 and spans an INCOMPLETE tri-quadratic polynomial space (Dofs are associated only to vertices and edges).
Note, Basis_HGRAD_HEX_I2_Serendipity_FEM = Basis_HGRAD_HEX_DEG2_FEM<true>
Basis functions are dual to a unisolvent set of degrees-of-freedom (DoF) defined and enumerated as follows:
\verbatim
=================================================================================================
Expand Down Expand Up @@ -139,13 +147,14 @@ namespace Intrepid2 {
namespace Impl {

/**
\brief See Intrepid2::Basis_HGRAD_HEX_C2_FEM
\brief See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM
*/
class Basis_HGRAD_HEX_C2_FEM {
template<bool serendipity>
class Basis_HGRAD_HEX_DEG2_FEM {
public:
typedef struct Hexahedron<27> cell_topology_type;
typedef struct Hexahedron<serendipity ? 20 : 27> cell_topology_type;
/**
\brief See Intrepid2::Basis_HGRAD_HEX_C2_FEM
\brief See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM
*/
template<EOperator opType>
struct Serial {
Expand All @@ -167,7 +176,7 @@ namespace Intrepid2 {
const EOperator operatorType);

/**
\brief See Intrepid2::Basis_HGRAD_HEX_C2_FEM
\brief See Intrepid2::Basis_HGRAD_HEX_DEG2_FEM
*/
template<typename outputValueViewType,
typename inputPointViewType,
Expand Down Expand Up @@ -209,26 +218,27 @@ namespace Intrepid2 {
opType != OPERATOR_D3 &&
opType != OPERATOR_D4 &&
opType != OPERATOR_MAX,
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::Serial::getValues) operator is not supported");
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::Serial::getValues) operator is not supported");
}
}
}
};
};
}

template<typename DeviceType = void,
typename outputValueType = double,
typename pointValueType = double>
class Basis_HGRAD_HEX_C2_FEM : public Basis<DeviceType,outputValueType,pointValueType> {
template<bool serendipity,
typename DeviceType,
typename outputValueType,
typename pointValueType>
class Basis_HGRAD_HEX_DEG2_FEM : public Basis<DeviceType,outputValueType,pointValueType> {
public:
using OrdinalTypeArray1DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray1DHost;
using OrdinalTypeArray2DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray2DHost;
using OrdinalTypeArray3DHost = typename Basis<DeviceType,outputValueType,pointValueType>::OrdinalTypeArray3DHost;

/** \brief Constructor.
*/
Basis_HGRAD_HEX_C2_FEM();
Basis_HGRAD_HEX_DEG2_FEM();

using OutputViewType = typename Basis<DeviceType,outputValueType,pointValueType>::OutputViewType;
using PointViewType = typename Basis<DeviceType,outputValueType,pointValueType>::PointViewType;
Expand All @@ -249,10 +259,16 @@ namespace Intrepid2 {
this->getBaseCellTopology(),
this->getCardinality() );
#endif
Impl::Basis_HGRAD_HEX_C2_FEM::
getValues<DeviceType>( outputValues,
inputPoints,
operatorType );
if constexpr (serendipity)
Impl::Basis_HGRAD_HEX_DEG2_FEM<true>::
getValues<DeviceType>( outputValues,
inputPoints,
operatorType );
else
Impl::Basis_HGRAD_HEX_DEG2_FEM<false>::
getValues<DeviceType>( outputValues,
inputPoints,
operatorType );
}

virtual
Expand All @@ -261,13 +277,13 @@ namespace Intrepid2 {
#ifdef HAVE_INTREPID2_DEBUG
// Verify rank of output array.
INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument,
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getDofCoords) rank = 2 required for dofCoords array");
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getDofCoords) rank = 2 required for dofCoords array");
// Verify 0th dimension of output array.
INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(dofCoords.extent(0)) != this->getCardinality(), std::invalid_argument,
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array");
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array");
// Verify 1st dimension of output array.
INTREPID2_TEST_FOR_EXCEPTION( dofCoords.extent(1) != this->getBaseCellTopology().getDimension(), std::invalid_argument,
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array");
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array");
#endif
Kokkos::deep_copy(dofCoords, this->dofCoords_);
}
Expand All @@ -278,18 +294,21 @@ namespace Intrepid2 {
#ifdef HAVE_INTREPID2_DEBUG
// Verify rank of output array.
INTREPID2_TEST_FOR_EXCEPTION( dofCoeffs.rank() != 1, std::invalid_argument,
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array");
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array");
// Verify 0th dimension of output array.
INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(dofCoeffs.extent(0)) != this->getCardinality(), std::invalid_argument,
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_C2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array");
">>> ERROR: (Intrepid2::Basis_HGRAD_HEX_DEG2_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array");
#endif
Kokkos::deep_copy(dofCoeffs, 1.0);
}

virtual
const char*
getName() const override {
return "Intrepid2_HGRAD_HEX_C2_FEM";
if constexpr (serendipity)
return "Intrepid2_HGRAD_HEX_I2_Serendipity_FEM";
else
return "Intrepid2_HGRAD_HEX_C2_FEM";
}

/** \brief returns the basis associated to a subCell.
Expand All @@ -307,16 +326,23 @@ namespace Intrepid2 {
Basis_HGRAD_LINE_C2_FEM<DeviceType,outputValueType,pointValueType>());
} else if(subCellDim == 2) {
return Teuchos::rcp(new
Basis_HGRAD_QUAD_C2_FEM<DeviceType,outputValueType,pointValueType>());
Basis_HGRAD_QUAD_DEG2_FEM<serendipity, DeviceType, outputValueType, pointValueType>());
}
INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
}

BasisPtr<typename Kokkos::HostSpace::device_type,outputValueType,pointValueType>
getHostBasis() const override{
return Teuchos::rcp(new Basis_HGRAD_HEX_C2_FEM<typename Kokkos::HostSpace::device_type,outputValueType,pointValueType>());
return Teuchos::rcp(new Basis_HGRAD_HEX_DEG2_FEM<serendipity, typename Kokkos::HostSpace::device_type,outputValueType,pointValueType>());
}
};

template<typename DeviceType = void, typename outputValueType = double, typename pointValueType = double>
using Basis_HGRAD_HEX_C2_FEM = Basis_HGRAD_HEX_DEG2_FEM<false, DeviceType, outputValueType, pointValueType>;

template<typename DeviceType = void, typename outputValueType = double, typename pointValueType = double>
using Basis_HGRAD_HEX_I2_Serendipity_FEM = Basis_HGRAD_HEX_DEG2_FEM<true, DeviceType, outputValueType, pointValueType>;

}// namespace Intrepid2

#include "Intrepid2_HGRAD_HEX_C2_FEMDef.hpp"
Expand Down
Loading

0 comments on commit ab899a0

Please sign in to comment.