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. ...