void Intrepid2FieldPattern::getSubcellNodes(const shards::CellTopology & cellTopo,unsigned dim,unsigned subCell,
std::vector<unsigned> & nodes)
{
if(dim==0) {
nodes.push_back(subCell);
return;
}
// get all nodes on requested sub cell
unsigned subCellNodeCount = cellTopo.getNodeCount(dim,subCell);
for(unsigned node=0;node<subCellNodeCount;++node)
nodes.push_back(cellTopo.getNodeMap(dim,subCell,node));
// sort them so they are ordered correctly for "includes" call
std::sort(nodes.begin(),nodes.end());
}
Basis_Constant_FEM<SpT,OT,PT>::
Basis_Constant_FEM(const shards::CellTopology cellTopo) {
const ordinal_type spaceDim = cellTopo.getDimension();
this->basisCardinality_ = 1;
this->basisDegree_ = 0;
this->basisCellTopology_ = cellTopo;
this->basisType_ = Intrepid2::BASIS_FEM_DEFAULT;
this->basisCoordinates_ = Intrepid2::COORDINATES_CARTESIAN;
// initialize tags
{
// Basis-dependent intializations
const ordinal_type tagSize = 4; // size of DoF tag, i.e., number of fields in the tag
const ordinal_type posScDim = 0; // position in the tag, counting from 0, of the subcell dim
const ordinal_type posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal
const ordinal_type posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell
// An array with local DoF tags assigned to the basis functions, in the order of their local enumeration
ordinal_type tags[4] = { spaceDim, 0, 0, 1 };
ordinal_type_array_1d_host tagView(&tags[0], 4);
this->setOrdinalTagData(this->tagToOrdinal_,
this->ordinalToTag_,
tagView,
this->basisCardinality_,
tagSize,
posScDim,
posScOrd,
posDfOrd);
}
// dofCoords on host and create its mirror view to device
Kokkos::DynRankView<typename scalarViewType::value_type,typename SpT::array_layout,Kokkos::HostSpace>
dofCoords("dofCoordsHost", this->basisCardinality_, spaceDim), cellVerts("cellVerts", spaceDim);
CellTools<SpT>::getReferenceCellCenter(Kokkos::subview(dofCoords, 0, Kokkos::ALL()),
cellVerts,
cellTopo);
this->dofCoords_ = Kokkos::create_mirror_view(typename SpT::memory_space(), dofCoords);
Kokkos::deep_copy(this->dofCoords_, dofCoords);
}
void
CellTools<SpT>::
mapToReferenceSubcell( /**/ Kokkos::DynRankView<refSubcellPointValueType,refSubcellPointProperties...> refSubcellPoints,
const Kokkos::DynRankView<paramPointValueType,paramPointProperties...> paramPoints,
const ordinal_type subcellDim,
const ordinal_type subcellOrd,
const shards::CellTopology parentCell ) {
#ifdef HAVE_INTREPID2_DEBUG
INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): the specified cell topology does not have a reference cell.");
INTREPID2_TEST_FOR_EXCEPTION( subcellDim != 1 &&
subcellDim != 2, std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-dimensional subcells.");
INTREPID2_TEST_FOR_EXCEPTION( subcellOrd < 0 ||
subcellOrd >= parentCell.getSubcellCount(subcellDim), std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): subcell ordinal out of range.");
// refSubcellPoints is rank-2 (P,D1), D1 = cell dimension
INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.rank() != 2, std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints must have rank 2.");
INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.dimension(1) != parentCell.getDimension(), std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (1) does not match to parent cell dimension.");
// paramPoints is rank-2 (P,D2) with D2 = subcell dimension
INTREPID2_TEST_FOR_EXCEPTION( paramPoints.rank() != 2, std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints must have rank 2.");
INTREPID2_TEST_FOR_EXCEPTION( paramPoints.dimension(1) != subcellDim, std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): paramPoints dimension (1) does not match to subcell dimension.");
// cross check: refSubcellPoints and paramPoints: dimension 0 must match
INTREPID2_TEST_FOR_EXCEPTION( refSubcellPoints.dimension(0) < paramPoints.dimension(0), std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): refSubcellPoints dimension (0) does not match to paramPoints dimension(0).");
#endif
const auto cellDim = parentCell.getDimension();
const auto numPts = paramPoints.dimension(0);
// Get the subcell map, i.e., the coefficients of the parametrization function for the subcell
// can i get this map from devices ?
subcellParamViewType subcellMap;
getSubcellParametrization( subcellMap,
subcellDim,
parentCell );
// subcell parameterization should be small computation (numPts is small) and it should be decorated with
// kokkos inline... let's not do this yet
// Apply the parametrization map to every point in parameter domain
switch (subcellDim) {
case 2: {
for (size_type pt=0;pt<numPts;++pt) {
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 (size_type i=0;i<cellDim;++i)
refSubcellPoints(pt, i) = subcellMap(subcellOrd, i, 0) + ( subcellMap(subcellOrd, i, 1)*u +
subcellMap(subcellOrd, i, 2)*v );
}
break;
}
case 1: {
for (size_type pt=0;pt<numPts;++pt) {
const auto u = paramPoints(pt, 0);
for (size_type i=0;i<cellDim;++i)
refSubcellPoints(pt, i) = subcellMap(subcellOrd, i, 0) + ( subcellMap(subcellOrd, i, 1)*u );
}
break;
}
default: {
INTREPID2_TEST_FOR_EXCEPTION( subcellDim != 1 &&
subcellDim != 2, std::invalid_argument,
">>> ERROR (Intrepid2::CellTools::mapToReferenceSubcell): method defined only for 1 and 2-subcells");
}
}
}
void mapToPhysicalFrame(ArrayPhysPoint & physPoints,
const ArrayRefPoint & refPoints,
const ArrayCell & cellWorkset,
const shards::CellTopology & cellTopo,
const int & whichCell)
{
int spaceDim = (int)cellTopo.getDimension();
int numCells = cellWorkset.dimension(0);
//points can be rank-2 (P,D), or rank-3 (C,P,D)
int numPoints = (refPoints.rank() == 2) ? refPoints.dimension(0) : refPoints.dimension(1);
/* // Mapping is computed using an appropriate H(grad) basis function: define RCP to the base class
Teuchos::RCP<Basis<Scalar, FieldContainer<Scalar> > > HGRAD_Basis;
// Choose the H(grad) basis depending on the cell topology. \todo define maps for shells and beams
switch( cellTopo.getKey() ){
// Standard Base topologies (number of cellWorkset = number of vertices)
case shards::Line<2>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_LINE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Triangle<3>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Quadrilateral<4>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Tetrahedron<4>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Hexahedron<8>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Wedge<6>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Pyramid<5>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_PYR_C1_FEM<Scalar, FieldContainer<Scalar> >() );
break;
// Standard Extended topologies
case shards::Triangle<6>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TRI_C2_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Quadrilateral<9>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_QUAD_C2_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Tetrahedron<10>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_C2_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Tetrahedron<11>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_TET_COMP12_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Hexahedron<27>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_HEX_C2_FEM<Scalar, FieldContainer<Scalar> >() );
break;
case shards::Wedge<18>::key:
HGRAD_Basis = Teuchos::rcp( new Basis_HGRAD_WEDGE_C2_FEM<Scalar, FieldContainer<Scalar> >() );
break;
// These extended topologies are not used for mapping purposes
case shards::Quadrilateral<8>::key:
case shards::Hexahedron<20>::key:
case shards::Wedge<15>::key:
TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
break;
// Base and Extended Line, Beam and Shell topologies
case shards::Line<3>::key:
case shards::Beam<2>::key:
case shards::Beam<3>::key:
case shards::ShellLine<2>::key:
case shards::ShellLine<3>::key:
case shards::ShellTriangle<3>::key:
case shards::ShellTriangle<6>::key:
case shards::ShellQuadrilateral<4>::key:
case shards::ShellQuadrilateral<8>::key:
case shards::ShellQuadrilateral<9>::key:
TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported. ");
break;
default:
TEUCHOS_TEST_FOR_EXCEPTION( (true), std::invalid_argument,
">>> ERROR (Intrepid::CellTools::mapToPhysicalFrame): Cell topology not supported.");
}// switch
// Temp (F,P) array for the values of nodal basis functions at the reference points
int basisCardinality = HGRAD_Basis -> getCardinality();
//.........这里部分代码省略.........
请发表评论