1 #ifndef INTREPID_HGRAD_HEX_CN_FEMDEF_HPP 
    2 #define INTREPID_HGRAD_HEX_CN_FEMDEF_HPP 
   52   template<
class Scalar, 
class ArrayScalar>
 
   56                                                                       const ArrayScalar &pts_x ,
 
   57                                                                       const ArrayScalar &pts_y ,
 
   58                                                                       const ArrayScalar &pts_z ):
 
   59     ptsx_( pts_x.dimension(0) , 1 ),
 
   60     ptsy_( pts_y.dimension(0) , 1 ),
 
   61     ptsz_( pts_z.dimension(0) , 1 )
 
   63     for (
int i=0;i<pts_x.dimension(0);i++)
 
   65         ptsx_(i,0) = pts_x(i,0);
 
   67     for (
int i=0;i<pts_y.dimension(0);i++)
 
   69         ptsy_(i,0) = pts_y(i,0);
 
   71     for (
int i=0;i<pts_z.dimension(0);i++)
 
   73         ptsz_(i,0) = pts_z(i,0);
 
   76     Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
 
   83     this->setBases( bases );
 
   86     if (orderx >= ordery && orderx >= orderz ) {
 
   89     else if (ordery >= orderx && ordery >= orderz) {
 
   95     this -> 
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
 
  101   template<
class Scalar, 
class ArrayScalar>
 
  103                                                                       const EPointType & pointType ):
 
  104     ptsx_( order+1 , 1 ),
 
  105     ptsy_( order+1 , 1 ),
 
  108     Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
 
  113     bases[0][1] = bases[0][0];
 
  114     bases[0][2] = bases[0][0];
 
  116     this->setBases( bases );
 
  120     this -> 
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
 
  126     EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND;
 
  127     PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ ,
 
  128                                                             shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
 
  132     for (
int i=0;i<order+1;i++)
 
  134         ptsy_(i,0) = ptsx_(i,0);
 
  135         ptsz_(i,0) = ptsx_(i,0);
 
  140   template<
class Scalar, 
class ArrayScalar>
 
  151     std::vector<int> tags( tagSize * this->getCardinality() );
 
  154     for (
int i=0;i<this->getCardinality();i++) {
 
  156        tags[tagSize*i+1] = 0;
 
  157        tags[tagSize*i+2] = i;
 
  158        tags[tagSize*i+3] = this->getCardinality();
 
  168     const std::vector<std::vector<int> >& xdoftags = xBasis_.
getAllDofTags();
 
  169     const std::vector<std::vector<int> >& ydoftags = yBasis_.
getAllDofTags();
 
  170     const std::vector<std::vector<int> >& zdoftags = zBasis_.
getAllDofTags();
 
  172     std::map<int,std::map<int,int> > total_dof_per_entity;
 
  173     std::map<int,std::map<int,int> > current_dof_per_entity;
 
  176     for (
int i=0;i<8;i++) {
 
  177       total_dof_per_entity[0][i] = 0;
 
  178       current_dof_per_entity[0][i] = 0;
 
  181     for (
int i=0;i<12;i++) {
 
  182       total_dof_per_entity[1][i] = 0;
 
  183       current_dof_per_entity[1][i] = 0;
 
  186     for (
int i=0;i<6;i++) {
 
  187       total_dof_per_entity[2][i] = 0;
 
  188       current_dof_per_entity[2][i] = 0;
 
  191     total_dof_per_entity[3][0] = 0;
 
  196       const int zdim = zdoftags[k][0];
 
  197       const int zent = zdoftags[k][1];
 
  199         const int ydim = ydoftags[j][0];
 
  200         const int yent = ydoftags[j][1];
 
  202           const int xdim = xdoftags[i][0];
 
  203           const int xent = xdoftags[i][1];
 
  207           total_dof_per_entity[dofdim][dofent] += 1;
 
  215       const int zdim = zdoftags[k][0];
 
  216       const int zent = zdoftags[k][1];
 
  218         const int ydim = ydoftags[j][0];
 
  219         const int yent = ydoftags[j][1];
 
  221           const int xdim = xdoftags[i][0];
 
  222           const int xent = xdoftags[i][1];
 
  226           tags[4*tagcur] = dofdim;
 
  227           tags[4*tagcur+1] = dofent;
 
  228           tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
 
  229           current_dof_per_entity[dofdim][dofent]++;
 
  230           tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
 
  236     Intrepid::setOrdinalTagData(
this -> tagToOrdinal_,
 
  237                                 this -> ordinalToTag_,
 
  239                                 this -> basisCardinality_,
 
  247   template<
class Scalar, 
class ArrayScalar>
 
  249                                                                const ArrayScalar &inputPoints ,
 
  250                                                                const EOperator operatorType )
 const  
  252 #ifdef HAVE_INTREPID_DEBUG 
  253     Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
 
  256                                                         this -> getBaseCellTopology(),
 
  257                                                         this -> getCardinality() );
 
  269     for (
int i=0;i<inputPoints.dimension(0);i++) {
 
  270       xInputPoints(i,0) = inputPoints(i,0);
 
  271       yInputPoints(i,0) = inputPoints(i,1);
 
  272       zInputPoints(i,0) = inputPoints(i,2);
 
  275     switch (operatorType) {
 
  282         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
 
  283         yBasis_.
getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
 
  284         zBasis_.
getValues(zBasisValues,zInputPoints,OPERATOR_VALUE);
 
  289             for (
int i=0;i<xBasis_.getCardinality();i++) {
 
  290               for (
int l=0;l<inputPoints.dimension(0);l++) {
 
  291                 outputValues(bfcur,l) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l);
 
  309         xBasis_.getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
 
  310         yBasis_.
getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
 
  311         zBasis_.
getValues(zBasisValues,zInputPoints,OPERATOR_VALUE);
 
  312         xBasis_.getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
 
  313         yBasis_.
getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);       
 
  314         zBasis_.
getValues(zBasisDerivs,zInputPoints,OPERATOR_D1);       
 
  319             for (
int i=0;i<xBasis_.getCardinality();i++) {
 
  320               for (
int l=0;l<inputPoints.dimension(0);l++) {
 
  321                 outputValues(bfcur,l,0) = xBasisDerivs(i,l,0) * yBasisValues(j,l) * zBasisValues(k,l);
 
  322                 outputValues(bfcur,l,1) = xBasisValues(i,l) * yBasisDerivs(j,l,0) * zBasisValues(k,l);
 
  323                 outputValues(bfcur,l,2) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisDerivs(k,l,0);
 
  345         Teuchos::Array<int> partialMult;
 
  347         for (
int d=0;d<getDkCardinality(operatorType,3);d++) {
 
  348           getDkMultiplicities( partialMult , d , operatorType , 3 );
 
  349           if (partialMult[0] == 0) {
 
  350             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
 
  351             xBasis_.getValues( xBasisValues , xInputPoints, OPERATOR_VALUE );
 
  354             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0),1);
 
  355             EOperator xop = (EOperator) ( (
int) OPERATOR_D1 + partialMult[0] - 1 );
 
  356             xBasis_.getValues( xBasisValues , xInputPoints, xop );
 
  357             xBasisValues.resize(xBasis_.getCardinality(),xInputPoints.dimension(0));
 
  359           if (partialMult[1] == 0) {
 
  360             yBasisValues.resize(yBasis_.
getCardinality(),yInputPoints.dimension(0));
 
  361             yBasis_.
getValues( yBasisValues , yInputPoints, OPERATOR_VALUE );
 
  364             yBasisValues.resize(yBasis_.
getCardinality(),yInputPoints.dimension(0),1);
 
  365             EOperator yop = (EOperator) ( (
int) OPERATOR_D1 + partialMult[1] - 1 );
 
  366             yBasis_.
getValues( yBasisValues , yInputPoints, yop );
 
  367             yBasisValues.resize(yBasis_.
getCardinality(),yInputPoints.dimension(0));
 
  369           if (partialMult[2] == 0) {
 
  370             zBasisValues.resize(zBasis_.
getCardinality(),zInputPoints.dimension(0));
 
  371             zBasis_.
getValues( zBasisValues , zInputPoints, OPERATOR_VALUE );
 
  374             zBasisValues.resize(zBasis_.
getCardinality(),zInputPoints.dimension(0),1);
 
  375             EOperator zop = (EOperator) ( (
int) OPERATOR_D1 + partialMult[2] - 1 );
 
  376             zBasis_.
getValues( zBasisValues , zInputPoints, zop );
 
  377             zBasisValues.resize(zBasis_.
getCardinality(),zInputPoints.dimension(0));
 
  384               for (
int i=0;i<xBasis_.getCardinality();i++) {
 
  385                 for (
int l=0;l<inputPoints.dimension(0);l++) {
 
  386                   outputValues(bfcur,l,d) = xBasisValues(i,l) * yBasisValues(j,l) * zBasisValues(k,l);
 
  396         TEUCHOS_TEST_FOR_EXCEPTION( 
true , std::invalid_argument,
 
  397                             ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): Operator type not implemented");
 
  402   template<
class Scalar,
class ArrayScalar>
 
  404                                                              const ArrayScalar &    inputPoints,
 
  405                                                              const ArrayScalar &    cellVertices,
 
  406                                                              const EOperator        operatorType)
 const {
 
  407   TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
 
  408                       ">>> ERROR (Basis_HGRAD_HEX_Cn_FEM): FEM Basis calling an FVD member function");
 
  411   template<
class Scalar,
class ArrayScalar>
 
  415     for (
int k=0;k<ptsz_.dimension(0);k++)
 
  417         for (
int j=0;j<ptsy_.dimension(0);j++)
 
  419             for (
int i=0;i<ptsx_.dimension(0);i++)
 
  421                 dofCoords(cur,0) = ptsx_(i,0);
 
  422                 dofCoords(cur,1) = ptsy_(j,0);
 
  423                 dofCoords(cur,2) = ptsz_(k,0);
 
virtual int getCardinality() const 
Returns cardinality of the basis. 
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1...
virtual void getDofCoords(ArrayScalar &DofCoords) const 
implement the DofCoordsInterface interface 
EBasis basisType_
Type of the basis. 
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags. 
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized 
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined. 
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays. 
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const 
Evaluation of a FEM basis on a reference Hexahedron cell. 
Basis_HGRAD_HEX_Cn_FEM(const int orderx, const int ordery, const int orderz, const ArrayScalar &pts_x, const ArrayScalar &pts_y, const ArrayScalar &pts_z)
Constructor. 
static void lineProduct3d(const int dim0, const int entity0, const int dim1, const int entity1, const int dim2, const int entity2, int &resultdim, int &resultentity)
virtual void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const =0
Evaluation of a FEM basis on a reference cell. 
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis. 
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom. ...