Skip to content

Commit

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

* trilinos-develop:
  Panzer MiniEM: Remove old code
  Panzer :: Evaluator for intrepid2's local L2 projection (trilinos#11811)
  Intrepid2: Symmetric quadrature rules and refactor of Projections (trilinos#12116)
  Zoltan2: Update GraphAdapters to use Kokkos data structures and improve API to return Host/Device Views
  • Loading branch information
@prwolfe
prwolfe committed Aug 14, 2023
2 parents c9cfda8 + 81cd5c6 commit 1dd3da3
Show file tree
Hide file tree
Showing 102 changed files with 32,370 additions and 4,997 deletions.
1 change: 0 additions & 1 deletion packages/intrepid2/assembly-examples/StructuredAssembly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,6 @@ Intrepid2::ScalarView<Scalar,DeviceType> performStructuredAssembly(Intrepid2::Ce
auto initialSetupTimer = Teuchos::TimeMonitor::getNewTimer("Initial Setup");
initialSetupTimer->start();
using namespace std;
using FunctionSpaceTools = FunctionSpaceTools<DeviceType>;
using IntegrationTools = IntegrationTools<DeviceType>;
// dimensions of the returned view are (C,F,F)

Expand Down
19 changes: 19 additions & 0 deletions packages/intrepid2/src/Cell/Intrepid2_CellTools.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1153,6 +1153,25 @@ namespace Intrepid2 {
const shards::CellTopology parentCell );


/** \brief Computes parameterization maps of 1- and 2-subcells of reference cells.
Overload of the previous function (see explanation above) where the subcell parametrization is used instead of
passing the parent cell topology.
\param refSubcellPoints [out] - rank-2 (P,D1) array with images of parameter space points
\param paramPoints [in] - rank-2 (P,D2) array with points in 1D or 2D parameter domain
\param subcellParametrization [in] - parametrization map of a subcell in the cell
\param subcellOrd [in] - subcell ordinal.
*/
template<typename refSubcellPointValueType, class ...refSubcellPointProperties,
typename paramPointValueType, class ...paramPointProperties>
static void
mapToReferenceSubcell( Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
const Kokkos::DynRankView<paramPointValueType,paramPointProperties...> paramPoints,
const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
const ordinal_type subcellOrd);


//============================================================================================//
// //
// Physical-to-reference frame mapping and its inverse //
Expand Down
61 changes: 37 additions & 24 deletions packages/intrepid2/src/Cell/Intrepid2_CellToolsDefRefToPhys.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,27 +112,25 @@ namespace Intrepid2 {
refSubcellViewType refSubcellPoints_;
const paramPointsViewType paramPoints_;
const subcellMapViewType subcellMap_;
ordinal_type subcellOrd_, dim_;
ordinal_type subcellOrd_;

KOKKOS_INLINE_FUNCTION
F_mapReferenceSubcell2( refSubcellViewType refSubcellPoints,
const paramPointsViewType paramPoints,
const subcellMapViewType subcellMap,
ordinal_type subcellOrd,
ordinal_type dim)
ordinal_type subcellOrd)
: refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
subcellOrd_(subcellOrd), dim_(dim){};
subcellOrd_(subcellOrd){};

KOKKOS_INLINE_FUNCTION
void operator()(const size_type pt) const {
void operator()(const size_type pt, const size_type d) const {

const auto u = paramPoints_(pt, 0);
const auto v = paramPoints_(pt, 1);

// map_dim(u,v) = c_0(dim) + c_1(dim)*u + c_2(dim)*v because both Quad and Tri ref faces are affine!
for (ordinal_type i=0;i<dim_;++i)
refSubcellPoints_(pt, i) = subcellMap_(subcellOrd_, i, 0) + ( subcellMap_(subcellOrd_, i, 1)*u +
subcellMap_(subcellOrd_, i, 2)*v );
refSubcellPoints_(pt, d) = subcellMap_(subcellOrd_, d, 0) + ( subcellMap_(subcellOrd_, d, 1)*u +
subcellMap_(subcellOrd_, d, 2)*v );
}
};

Expand All @@ -143,22 +141,20 @@ namespace Intrepid2 {
refSubcellViewType refSubcellPoints_;
const paramPointsViewType paramPoints_;
const subcellMapViewType subcellMap_;
ordinal_type subcellOrd_, dim_;
ordinal_type subcellOrd_;

KOKKOS_INLINE_FUNCTION
F_mapReferenceSubcell1( refSubcellViewType refSubcellPoints,
const paramPointsViewType paramPoints,
const subcellMapViewType subcellMap,
ordinal_type subcellOrd,
ordinal_type dim)
ordinal_type subcellOrd)
: refSubcellPoints_(refSubcellPoints), paramPoints_(paramPoints), subcellMap_(subcellMap),
subcellOrd_(subcellOrd), dim_(dim){};
subcellOrd_(subcellOrd){};

KOKKOS_INLINE_FUNCTION
void operator()(const size_type pt) const {
void operator()(const size_type pt, const size_type d) const {
const auto u = paramPoints_(pt, 0);
for (ordinal_type i=0;i<dim_;++i)
refSubcellPoints_(pt, i) = subcellMap_(subcellOrd_, i, 0) + ( subcellMap_(subcellOrd_, i, 1)*u );
refSubcellPoints_(pt, d) = subcellMap_(subcellOrd_, d, 0) + ( subcellMap_(subcellOrd_, d, 1)*u );
}
};
}
Expand Down Expand Up @@ -237,8 +233,8 @@ namespace Intrepid2 {
const ordinal_type subcellDim,
const ordinal_type subcellOrd,
const shards::CellTopology parentCell ) {
ordinal_type parentCellDim = parentCell.getDimension();
#ifdef HAVE_INTREPID2_DEBUG
ordinal_type parentCellDim = parentCell.getDimension();
INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");

Expand Down Expand Up @@ -267,6 +263,22 @@ namespace Intrepid2 {
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
#endif

// Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
const auto subcellMap = RefSubcellParametrization<DeviceType>::get(subcellDim, parentCell.getKey());

mapToReferenceSubcell(refSubcellPoints, paramPoints, subcellMap, subcellOrd);
}

template<typename DeviceType>
template<typename refSubcellPointValueType, class ...refSubcellPointProperties,
typename paramPointValueType, class ...paramPointProperties>
void
CellTools<DeviceType>::
mapToReferenceSubcell( Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
const Kokkos::DynRankView<paramPointValueType,paramPointProperties...> paramPoints,
const typename RefSubcellParametrization<DeviceType>::ConstViewType subcellParametrization,
const ordinal_type subcellOrd) {

constexpr bool are_accessible =
Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
typename decltype(refSubcellPoints)::memory_space>::accessible &&
Expand All @@ -275,29 +287,30 @@ namespace Intrepid2 {

static_assert(are_accessible, "CellTools<DeviceType>::mapToReferenceSubcell(..): input/output views' memory spaces are not compatible with DeviceType");

// Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
const auto subcellMap = RefSubcellParametrization<DeviceType>::get(subcellDim, parentCell.getKey());

const ordinal_type parentCellDim = subcellParametrization.extent(1);
const ordinal_type numPts = paramPoints.extent(0);
const ordinal_type subcellDim = paramPoints.extent(1);

//Note: this function has several template parameters and the compiler gets confused if using a lambda function
using FunctorType1 = FunctorCellTools::F_mapReferenceSubcell1<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellMap)>;
using FunctorType2 = FunctorCellTools::F_mapReferenceSubcell2<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellMap)>;
using FunctorType1 = FunctorCellTools::F_mapReferenceSubcell1<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization)>;
using FunctorType2 = FunctorCellTools::F_mapReferenceSubcell2<decltype(refSubcellPoints), decltype(paramPoints), decltype(subcellParametrization)>;

Kokkos::RangePolicy<ExecSpaceType> policy(0, numPts);

Kokkos::MDRangePolicy<ExecSpaceType,Kokkos::Rank<2>> rangePolicy({0,0},{numPts,parentCellDim});
// Apply the parametrization map to every point in parameter domain
switch (subcellDim) {
case 2: {
Kokkos::parallel_for( policy, FunctorType2(refSubcellPoints, paramPoints, subcellMap, subcellOrd, parentCellDim) );
Kokkos::parallel_for( rangePolicy, FunctorType2(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
break;
}
case 1: {
Kokkos::parallel_for( policy, FunctorType1(refSubcellPoints, paramPoints, subcellMap, subcellOrd, parentCellDim) );
Kokkos::parallel_for( rangePolicy, FunctorType1(refSubcellPoints, paramPoints, subcellParametrization, subcellOrd) );
break;
}
default: {}
}
}

}

#endif
20 changes: 8 additions & 12 deletions packages/intrepid2/src/Cell/Intrepid2_ProjectedGeometry.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,21 +93,19 @@ namespace Intrepid2
INTREPID2_TEST_FOR_EXCEPTION(projectedBasisNodes.extent_int(2) != spaceDim, std::invalid_argument, "projectedBasisNodes must have shape (C,F,D)");

using ExecutionSpace = typename DeviceType::execution_space;
using ProjectionTools = Experimental::ProjectionTools<ExecutionSpace>; // TODO: when ProjectionTools supports it, replace template argument with DeviceType
using ProjectionStruct = Experimental::ProjectionStruct<ExecutionSpace,PointScalar>; // TODO: when ProjectionTools supports it, replace template argument with DeviceType
using ProjectionTools = Experimental::ProjectionTools<DeviceType>;
using ProjectionStruct = Experimental::ProjectionStruct<DeviceType,PointScalar>;

ProjectionStruct projectionStruct;
ordinal_type targetQuadratureDegree(targetHGradBasis->getDegree()), targetDerivativeQuadratureDegree(targetHGradBasis->getDegree());
projectionStruct.createHGradProjectionStruct(targetHGradBasis, targetQuadratureDegree, targetDerivativeQuadratureDegree);

const ordinal_type numPoints = projectionStruct.getNumTargetEvalPoints();
const ordinal_type numGradPoints = projectionStruct.getNumTargetDerivEvalPoints();

ViewType evaluationPointsRefSpace ("ProjectedGeometry evaluation points ref space (value)", numCells, numPoints, spaceDim);
ViewType evaluationGradPointsRefSpace("ProjectedGeometry evaluation points ref space (gradient)", numCells, numGradPoints, spaceDim);


auto evaluationPointsRefSpace = projectionStruct.getAllEvalPoints();
auto evaluationGradPointsRefSpace = projectionStruct.getAllDerivEvalPoints();
const ordinal_type numPoints = evaluationPointsRefSpace.extent(0);
const ordinal_type numGradPoints = evaluationGradPointsRefSpace.extent(0);

auto elementOrientations = flatCellGeometry.getOrientations();
ProjectionTools::getHGradEvaluationPoints(evaluationPointsRefSpace, evaluationGradPointsRefSpace, elementOrientations, targetHGradBasis.get(), &projectionStruct);

// printFunctor1(elementOrientations, std::cout);

Expand Down Expand Up @@ -179,8 +177,6 @@ namespace Intrepid2
ProjectionTools::getHGradBasisCoeffs(projectedBasisNodesForComp,
evaluationValues,
transformedGradientData.getUnderlyingView(),
evaluationPointsRefSpace,
evaluationGradPointsRefSpace,
elementOrientations,
targetHGradBasis.get(),
&projectionStruct);
Expand Down
4 changes: 2 additions & 2 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HCURL_TRI_In_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@
#define __INTREPID2_HCURL_TRI_IN_FEM_DEF_HPP__

#include "Intrepid2_HGRAD_TRI_Cn_FEM_ORTH.hpp"
#include "Intrepid2_CubatureDirectTriDefault.hpp"
#include "Intrepid2_CubatureDirectTriSymmetric.hpp"

namespace Intrepid2 {

Expand Down Expand Up @@ -249,7 +249,7 @@ namespace Intrepid2 {

// now I need to integrate { (x,y) \times phi } against the big basis
// first, get a cubature rule.
CubatureDirectTriDefault<Kokkos::HostSpace::execution_space,scalarType,scalarType> myCub( 2 * order );
CubatureDirectTriSymmetric<Kokkos::HostSpace::execution_space,scalarType,scalarType> myCub( 2 * order );
Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubPoints("Hcurl::Tri::In::cubPoints", myCub.getNumPoints() , spaceDim );
Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubWeights("Hcurl::Tri::In::cubWeights", myCub.getNumPoints() );
myCub.getCubature( cubPoints , cubWeights );
Expand Down
4 changes: 2 additions & 2 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_TRI_In_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@
#define __INTREPID2_HDIV_TRI_IN_FEM_DEF_HPP__

#include "Intrepid2_HGRAD_TRI_Cn_FEM_ORTH.hpp"
#include "Intrepid2_CubatureDirectTriDefault.hpp"
#include "Intrepid2_CubatureDirectTriSymmetric.hpp"

namespace Intrepid2 {

Expand Down Expand Up @@ -247,7 +247,7 @@ Basis_HDIV_TRI_In_FEM( const ordinal_type order,

// now I need to integrate { (x,y) phi } against the big basis
// first, get a cubature rule.
CubatureDirectTriDefault<Kokkos::HostSpace::execution_space,scalarType,scalarType> myCub( 2 * order );
CubatureDirectTriSymmetric<Kokkos::HostSpace::execution_space,scalarType,scalarType> myCub( 2 * order );
Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubPoints("Hdiv::Tri::In::cubPoints", myCub.getNumPoints() , spaceDim );
Kokkos::DynRankView<scalarType,typename DT::execution_space::array_layout,Kokkos::HostSpace> cubWeights("Hdiv::Tri::In::cubWeights", myCub.getNumPoints() );
myCub.getCubature( cubPoints , cubWeights );
Expand Down
1 change: 1 addition & 0 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HDIV_WEDGE_I1_FEMDef.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -214,6 +214,7 @@ namespace Intrepid2 {
dofCoeffs(4,0) = 0.0; dofCoeffs(4,1) = 0.0; dofCoeffs(4,2) = 1.0;

this->dofCoeffs_ = Kokkos::create_mirror_view(typename DT::memory_space(), dofCoeffs);
Kokkos::deep_copy(this->dofCoeffs_, dofCoeffs);
}

}// namespace Intrepid2
Expand Down
2 changes: 2 additions & 0 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_HierarchicalBasis_HDIV_TRI.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,8 @@ namespace Intrepid2

using BasisBase = Basis<DeviceType,OutputScalar,PointScalar>;

using HostBasis = HierarchicalBasis_HDIV_TRI<typename Kokkos::HostSpace::device_type, OutputScalar, PointScalar, useCGBasis>;

using typename BasisBase::OrdinalTypeArray1DHost;
using typename BasisBase::OrdinalTypeArray2DHost;

Expand Down
12 changes: 6 additions & 6 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_TensorBasis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ namespace Intrepid2
{
template<ordinal_type spaceDim>
KOKKOS_INLINE_FUNCTION
ordinal_type getDkEnumeration(Kokkos::Array<int,spaceDim> &entries);
ordinal_type getDkEnumeration(const Kokkos::Array<int,spaceDim> &entries);

template<ordinal_type spaceDim>
KOKKOS_INLINE_FUNCTION
Expand Down Expand Up @@ -147,7 +147,7 @@ namespace Intrepid2

// ensure that we don't try to allocate an empty array…
constexpr ordinal_type sizeForSubArray = (spaceDim > 2) ? spaceDim - 1 : 1;
Kokkos::Array<int,sizeForSubArray> subEntries;
Kokkos::Array<int,sizeForSubArray> subEntries = {};

// the -1 in sub-entry enumeration value accounts for the fact that the entry is the one *after* (k0,0,…,0)
getDkEnumerationInverse<spaceDim-1>(subEntries, dkEnum - dkEnumFor_k0 - 1, operatorOrder - entries[0]);
Expand All @@ -164,14 +164,14 @@ namespace Intrepid2

template<>
KOKKOS_INLINE_FUNCTION
ordinal_type getDkEnumeration<1>(Kokkos::Array<int,1> &entries)
ordinal_type getDkEnumeration<1>(const Kokkos::Array<int,1> &entries)
{
return getDkEnumeration<1>(entries[0]);
}

template<ordinal_type spaceDim>
KOKKOS_INLINE_FUNCTION
ordinal_type getDkEnumeration(Kokkos::Array<int,spaceDim> &entries)
ordinal_type getDkEnumeration(const Kokkos::Array<int,spaceDim> &entries)
{
ordinal_type k_minus_k0 = 0; // sum of all the entries but the first

Expand Down Expand Up @@ -203,10 +203,10 @@ namespace Intrepid2
ordinal_type getDkTensorIndex(const ordinal_type dkEnum1, const ordinal_type operatorOrder1,
const ordinal_type dkEnum2, const ordinal_type operatorOrder2)
{
Kokkos::Array<int,spaceDim1> entries1;
Kokkos::Array<int,spaceDim1> entries1 = {};
getDkEnumerationInverse<spaceDim1>(entries1, dkEnum1, operatorOrder1);

Kokkos::Array<int,spaceDim2> entries2;
Kokkos::Array<int,spaceDim2> entries2 = {};
getDkEnumerationInverse<spaceDim2>(entries2, dkEnum2, operatorOrder2);

const int spaceDim = spaceDim1 + spaceDim2;
Expand Down
2 changes: 1 addition & 1 deletion packages/intrepid2/src/Discretization/FunctionSpaceTools/Intrepid2_FunctionSpaceTools.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -666,7 +666,7 @@ namespace Intrepid2 {
\f$DF\f$ is the Jacobian of the map \f$F\f$.
\param outputVals [out] - Output array with reference HCURL data, indexed by (C,P,D)
\param jacobianInverse [in] - Input array containing the Jacobian, indexed by (C,P,D,D)
\param jacobian [in] - Input array containing the Jacobian, indexed by (C,P,D,D)
\param inputVals [in] - Input array of HCURL data in physical space, indexed by (C,P,D).
*/
Expand Down
Loading

0 comments on commit 1dd3da3

Please sign in to comment.