1#ifndef INTREPID_HGRAD_QUAD_CN_FEMDEF_HPP
2#define INTREPID_HGRAD_QUAD_CN_FEMDEF_HPP
55 template<
class Scalar,
class ArrayScalar>
57 const ArrayScalar &pts_x ,
58 const ArrayScalar &pts_y ):
59 ptsx_( pts_x.dimension(0) , 1 ) ,
60 ptsy_( pts_y.dimension(0) , 1 )
62 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
66 this->setBases( bases );
69 if (orderx > ordery) {
75 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
80 for (
int i=0;i<pts_x.dimension(0);i++)
82 ptsx_(i,0) = pts_x(i,0);
85 for (
int i=0;i<pts_y.dimension(0);i++)
87 ptsy_(i,0) = pts_y(i,0);
92 template<
class Scalar,
class ArrayScalar>
94 const EPointType &pointType ):
95 ptsx_( order+1 , 1 ) ,
98 Array<Array<RCP<Basis<Scalar,ArrayScalar> > > > bases(1);
102 this->setBases( bases );
106 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
112 EPointType pt = (pointType==POINTTYPE_EQUISPACED)?pointType:POINTTYPE_WARPBLEND;
113 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( ptsx_ ,
114 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
119 for (
int i=0;i<order+1;i++)
121 ptsy_(i,0) = ptsx_(i,0);
125 template<
class Scalar,
class ArrayScalar>
136 std::vector<int> tags( tagSize * this->getCardinality() );
139 for (
int i=0;i<this->getCardinality();i++) {
141 tags[tagSize*i+1] = 0;
142 tags[tagSize*i+2] = i;
143 tags[tagSize*i+3] = this->getCardinality();
151 const std::vector<std::vector<int> >& xdoftags = xBasis_.
getAllDofTags();
152 const std::vector<std::vector<int> >& ydoftags = yBasis_.
getAllDofTags();
154 std::map<int,std::map<int,int> > total_dof_per_entity;
155 std::map<int,std::map<int,int> > current_dof_per_entity;
157 for (
int i=0;i<4;i++) {
158 total_dof_per_entity[0][i] = 0;
159 current_dof_per_entity[0][i] = 0;
161 for (
int i=0;i<4;i++) {
162 total_dof_per_entity[1][i] = 0;
163 current_dof_per_entity[1][i] = 0;
165 total_dof_per_entity[2][0] = 0;
166 current_dof_per_entity[2][0] = 0;
170 const int ydim = ydoftags[j][0];
171 const int yent = ydoftags[j][1];
173 const int xdim = xdoftags[i][0];
174 const int xent = xdoftags[i][1];
178 total_dof_per_entity[dofdim][dofent] += 1;
184 const int ydim = ydoftags[j][0];
185 const int yent = ydoftags[j][1];
187 const int xdim = xdoftags[i][0];
188 const int xent = xdoftags[i][1];
192 tags[4*tagcur] = dofdim;
193 tags[4*tagcur+1] = dofent;
194 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
195 current_dof_per_entity[dofdim][dofent]++;
196 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
202 this -> ordinalToTag_,
204 this -> basisCardinality_,
212 template<
class Scalar,
class ArrayScalar>
214 const ArrayScalar &inputPoints ,
215 const EOperator operatorType )
const
217#ifdef HAVE_INTREPID_DEBUG
218 getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues,
221 this -> getBaseCellTopology(),
222 this -> getCardinality() );
231 for (
int i=0;i<inputPoints.dimension(0);i++) {
232 xInputPoints(i,0) = inputPoints(i,0);
233 yInputPoints(i,0) = inputPoints(i,1);
236 switch (operatorType) {
242 xBasis_.
getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
243 yBasis_.
getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
248 for (
int k=0;k<inputPoints.dimension(0);k++) {
249 outputValues(bfcur,k) = xBasisValues(i,k) * yBasisValues(j,k);
264 xBasis_.
getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
265 yBasis_.
getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
266 xBasis_.
getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
267 yBasis_.
getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);
274 for (
int k=0;k<inputPoints.dimension(0);k++) {
275 outputValues(bfcur,k,0) = xBasisDerivs(i,k,0) * yBasisValues(j,k);
276 outputValues(bfcur,k,1) = xBasisValues(i,k) * yBasisDerivs(j,k,0);
296 Teuchos::Array<int> partialMult;
300 if (partialMult[0] == 0) {
302 xBasis_.
getValues( xBasisValues , xInputPoints, OPERATOR_VALUE );
306 EOperator xop = (EOperator) ( (
int) OPERATOR_D1 + partialMult[0] - 1 );
307 xBasis_.
getValues( xBasisValues , xInputPoints, xop );
310 if (partialMult[1] == 0) {
312 yBasis_.
getValues( yBasisValues , yInputPoints, OPERATOR_VALUE );
316 EOperator yop = (EOperator) ( (
int) OPERATOR_D1 + partialMult[1] - 1 );
317 yBasis_.
getValues( yBasisValues , yInputPoints, yop );
325 for (
int k=0;k<inputPoints.dimension(0);k++) {
326 outputValues(bfcur,k,d) = xBasisValues(i,k) * yBasisValues(j,k);
341 xBasis_.
getValues(xBasisValues,xInputPoints,OPERATOR_VALUE);
342 yBasis_.
getValues(yBasisValues,yInputPoints,OPERATOR_VALUE);
343 xBasis_.
getValues(xBasisDerivs,xInputPoints,OPERATOR_D1);
344 yBasis_.
getValues(yBasisDerivs,yInputPoints,OPERATOR_D1);
351 for (
int k=0;k<inputPoints.dimension(0);k++) {
352 outputValues(bfcur,k,0) = xBasisValues(i,k) * yBasisDerivs(j,k,0);
353 outputValues(bfcur,k,1) = -xBasisDerivs(i,k,0) * yBasisValues(j,k);
361 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
362 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): Operator type not implemented");
367 template<
class Scalar,
class ArrayScalar>
369 const ArrayScalar & inputPoints,
370 const ArrayScalar & cellVertices,
371 const EOperator operatorType)
const {
372 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
373 ">>> ERROR (Basis_HGRAD_QUAD_Cn_FEM): FEM Basis calling an FVD member function");
376 template<
class Scalar,
class ArrayScalar>
380 for (
int j=0;j<ptsy_.dimension(0);j++)
382 for (
int i=0;i<ptsx_.dimension(0);i++)
384 dofCoords(cur,0) = ptsx_(i,0);
385 dofCoords(cur,1) = ptsy_(j,0);
void setOrdinalTagData(std::vector< std::vector< std::vector< int > > > &tagToOrdinal, std::vector< std::vector< int > > &ordinalToTag, const int *tags, const int basisCard, const int tagSize, const int posScDim, const int posScOrd, const int posDfOrd)
Fills ordinalToTag_ and tagToOrdinal_ by basis-specific tag data.
int getDkCardinality(const EOperator operatorType, const int spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
void getDkMultiplicities(Teuchos::Array< int > &partialMult, const int derivativeEnum, const EOperator operatorType, const int spaceDim)
Returns multiplicities of dx, dy, and dz based on the enumeration of the partial derivative,...
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
virtual void getDofCoords(ArrayScalar &dofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Quadrilateral cell.
Basis_HGRAD_QUAD_Cn_FEM(const int orderx, const int ordery, const ArrayScalar &pts_x, const ArrayScalar &pts_y)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
virtual const std::vector< std::vector< int > > & getAllDofTags()
Retrieves all DoF tags.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const =0
Evaluation of a FEM basis on a reference cell.
int basisCardinality_
Cardinality of the basis, i.e., the number of basis functions/degrees-of-freedom.
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
EBasis basisType_
Type of the basis.
int basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos....
virtual int getCardinality() const
Returns cardinality of the basis.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0.
int dimension(const int whichDim) const
Returns the specified dimension.
static void lineProduct2d(const int dim0, const int entity0, const int dim1, const int entity1, int &resultdim, int &resultentity)