51 template<
class Scalar,
class ArrayScalar>
53 const ArrayScalar & ptsClosed ,
54 const ArrayScalar & ptsOpen):
55 closedBasis_( order , ptsClosed ),
56 openBasis_( order-1 , ptsOpen ),
57 closedPts_( ptsClosed ),
61 this ->
basisCardinality_ = 3 * closedBasis_.getCardinality() * openBasis_.getCardinality() * openBasis_.getCardinality();
62 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
67 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3);
68 bases[0].resize(3); bases[1].resize(3); bases[2].resize(3);
69 bases[0][0] = rcp( &closedBasis_ ,
false );
70 bases[0][1] = rcp( &openBasis_ ,
false );
71 bases[0][2] = rcp( &openBasis_ ,
false );
72 bases[1][0] = rcp( &openBasis_ ,
false );
73 bases[1][1] = rcp( &closedBasis_ ,
false );
74 bases[1][2] = rcp( &openBasis_ ,
false );
75 bases[2][0] = rcp( &openBasis_ ,
false );
76 bases[2][1] = rcp( &openBasis_ ,
false );
77 bases[2][2] = rcp( &closedBasis_ ,
false );
78 this->setBases( bases );
82 template<
class Scalar,
class ArrayScalar>
84 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
85 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ),
86 closedPts_( order+1 , 1 ),
90 this ->
basisCardinality_ = 3 * closedBasis_.getCardinality() * openBasis_.getCardinality() * openBasis_.getCardinality();
91 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
96 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
97 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
100 pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
102 if (pointType == POINTTYPE_SPECTRAL)
104 PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
109 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
110 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
113 POINTTYPE_EQUISPACED );
118 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3);
119 bases[0].resize(3); bases[1].resize(3); bases[2].resize(3);
120 bases[0][0] = rcp( &closedBasis_ ,
false );
121 bases[0][1] = rcp( &openBasis_ ,
false );
122 bases[0][2] = rcp( &openBasis_ ,
false );
123 bases[1][0] = rcp( &openBasis_ ,
false );
124 bases[1][1] = rcp( &closedBasis_ ,
false );
125 bases[1][2] = rcp( &openBasis_ ,
false );
126 bases[2][0] = rcp( &openBasis_ ,
false );
127 bases[2][1] = rcp( &openBasis_ ,
false );
128 bases[2][2] = rcp( &closedBasis_ ,
false );
129 this->setBases( bases );
135 template<
class Scalar,
class ArrayScalar>
144 std::vector<int> tags( tagSize * this->getCardinality() );
146 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
147 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
149 std::map<int,std::map<int,int> > total_dof_per_entity;
150 std::map<int,std::map<int,int> > current_dof_per_entity;
153 for (
int i=0;i<4;i++) {
154 total_dof_per_entity[0][i] = 0;
155 current_dof_per_entity[0][i] = 0;
158 for (
int i=0;i<12;i++) {
159 total_dof_per_entity[1][i] = 0;
160 current_dof_per_entity[1][i] = 0;
163 for (
int i=0;i<6;i++) {
164 total_dof_per_entity[2][i] = 0;
165 current_dof_per_entity[2][i] = 0;
167 total_dof_per_entity[3][0] = 0;
168 current_dof_per_entity[3][0] = 0;
171 for (
int i=0;i<6;i++) {
172 total_dof_per_entity[2][i] = openBasis_.getCardinality() * openBasis_.getCardinality();
175 total_dof_per_entity[2][0] = this->getCardinality() - 6 * openBasis_.getCardinality()* openBasis_.getCardinality();
180 for (
int k=0;k<openBasis_.getCardinality();k++) {
181 const int odimk = openDofTags[k][0];
182 const int oentk = openDofTags[k][1];
183 for (
int j=0;j<openBasis_.getCardinality();j++) {
184 const int odimj = openDofTags[j][0];
185 const int oentj = openDofTags[j][1];
186 for (
int i=0;i<closedBasis_.getCardinality();i++) {
187 const int cdim = closedDofTags[i][0];
188 const int cent = closedDofTags[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];
204 for (
int k=0;k<openBasis_.getCardinality();k++) {
205 const int odimk = openDofTags[k][0];
206 const int oentk = openDofTags[k][1];
207 for (
int j=0;j<closedBasis_.getCardinality();j++) {
208 const int cdim = closedDofTags[j][0];
209 const int cent = closedDofTags[j][1];
210 for (
int i=0;i<openBasis_.getCardinality();i++) {
211 const int odimi = openDofTags[i][0];
212 const int oenti = openDofTags[i][1];
216 tags[4*tagcur] = dofdim;
217 tags[4*tagcur+1] = dofent;
218 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
219 current_dof_per_entity[dofdim][dofent]++;
220 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
228 for (
int k=0;k<closedBasis_.getCardinality();k++) {
229 const int cdim = closedDofTags[k][0];
230 const int cent = closedDofTags[k][1];
231 for (
int j=0;j<openBasis_.getCardinality();j++) {
232 const int odimj = openDofTags[j][0];
233 const int oentj = openDofTags[j][1];
234 for (
int i=0;i<openBasis_.getCardinality();i++) {
235 const int odimi = openDofTags[i][0];
236 const int oenti = openDofTags[i][1];
240 tags[4*tagcur] = dofdim;
241 tags[4*tagcur+1] = dofent;
242 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
243 current_dof_per_entity[dofdim][dofent]++;
244 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
258 Intrepid::setOrdinalTagData(
this -> tagToOrdinal_,
259 this -> ordinalToTag_,
261 this -> basisCardinality_,
269 template<
class Scalar,
class ArrayScalar>
271 const ArrayScalar & inputPoints,
272 const EOperator operatorType)
const {
275 #ifdef HAVE_INTREPID_DEBUG
276 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
279 this -> getBaseCellTopology(),
280 this -> getCardinality() );
284 int dim0 = inputPoints.dimension(0);
292 for (
int i=0;i<dim0;i++) {
293 xPoints(i,0) = inputPoints(i,0);
294 yPoints(i,0) = inputPoints(i,1);
295 zPoints(i,0) = inputPoints(i,2);
298 switch (operatorType) {
308 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
309 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
310 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE );
311 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
312 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
313 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE );
317 for (
int k=0;k<openBasis_.getCardinality();k++) {
318 for (
int j=0;j<openBasis_.getCardinality();j++) {
319 for (
int i=0;i<closedBasis_.getCardinality();i++) {
320 for (
int l=0;l<dim0;l++) {
321 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * openBasisValsZPts(k,l);
322 outputValues(bfcur,l,1) = 0.0;
323 outputValues(bfcur,l,2) = 0.0;
331 for (
int k=0;k<openBasis_.getCardinality();k++) {
332 for (
int j=0;j<closedBasis_.getCardinality();j++) {
333 for (
int i=0;i<openBasis_.getCardinality();i++) {
334 for (
int l=0;l<dim0;l++) {
335 outputValues(bfcur,l,0) = 0.0;
336 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l);
337 outputValues(bfcur,l,2) = 0.0;
345 for (
int k=0;k<closedBasis_.getCardinality();k++) {
346 for (
int j=0;j<openBasis_.getCardinality();j++) {
347 for (
int i=0;i<openBasis_.getCardinality();i++) {
348 for (
int l=0;l<dim0;l++) {
349 outputValues(bfcur,l,0) = 0.0;
350 outputValues(bfcur,l,1) = 0.0;
351 outputValues(bfcur,l,2) = openBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l);
370 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
371 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
372 closedBasis_.getValues( closedBasisDerivsZPts , zPoints , OPERATOR_D1 );
373 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
374 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
375 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE );
380 for (
int k=0;k<openBasis_.getCardinality();k++) {
381 for (
int j=0;j<openBasis_.getCardinality();j++) {
382 for (
int i=0;i<closedBasis_.getCardinality();i++) {
383 for (
int l=0;l<dim0;l++) {
384 outputValues(bfcur,l) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l) * openBasisValsZPts(k,l);
392 for (
int k=0;k<openBasis_.getCardinality();k++) {
393 for (
int j=0;j<closedBasis_.getCardinality();j++) {
394 for (
int i=0;i<openBasis_.getCardinality();i++) {
395 for (
int l=0;l<dim0;l++) {
396 outputValues(bfcur,l) = openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * openBasisValsZPts(k,l);
404 for (
int k=0;k<closedBasis_.getCardinality();k++) {
405 for (
int j=0;j<openBasis_.getCardinality();j++) {
406 for (
int i=0;i<openBasis_.getCardinality();i++) {
407 for (
int l=0;l<dim0;l++) {
408 outputValues(bfcur,l) = openBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0);
417 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument,
418 ">>> ERROR (Basis_HDIV_HEX_In_FEM): CURL is invalid operator for HDIV Basis Functions");
422 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
423 ">>> ERROR (Basis_HDIV_HEX_In_FEM): GRAD is invalid operator for HDIV Basis Functions");
436 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
437 (operatorType == OPERATOR_D2) ||
438 (operatorType == OPERATOR_D3) ||
439 (operatorType == OPERATOR_D4) ||
440 (operatorType == OPERATOR_D5) ||
441 (operatorType == OPERATOR_D6) ||
442 (operatorType == OPERATOR_D7) ||
443 (operatorType == OPERATOR_D8) ||
444 (operatorType == OPERATOR_D9) ||
445 (operatorType == OPERATOR_D10) ),
446 std::invalid_argument,
447 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Invalid operator type");
451 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
452 (operatorType != OPERATOR_GRAD) &&
453 (operatorType != OPERATOR_CURL) &&
454 (operatorType != OPERATOR_DIV) &&
455 (operatorType != OPERATOR_D1) &&
456 (operatorType != OPERATOR_D2) &&
457 (operatorType != OPERATOR_D3) &&
458 (operatorType != OPERATOR_D4) &&
459 (operatorType != OPERATOR_D5) &&
460 (operatorType != OPERATOR_D6) &&
461 (operatorType != OPERATOR_D7) &&
462 (operatorType != OPERATOR_D8) &&
463 (operatorType != OPERATOR_D9) &&
464 (operatorType != OPERATOR_D10) ),
465 std::invalid_argument,
466 ">>> ERROR (Basis_HDIV_HEX_In_FEM): Invalid operator type");
472 template<
class Scalar,
class ArrayScalar>
474 const ArrayScalar & inputPoints,
475 const ArrayScalar & cellVertices,
476 const EOperator operatorType)
const {
477 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
478 ">>> ERROR (Basis_HDIV_HEX_In_FEM): FEM Basis calling an FVD member function");
481 template<
class Scalar,
class ArrayScalar>
487 for (
int k=0;k<openPts_.dimension(0);k++)
489 for (
int j=0;j<openPts_.dimension(0);j++)
491 for (
int i=0;i<closedPts_.dimension(0);i++)
493 DofCoords(cur,0) = closedPts_(i,0);
494 DofCoords(cur,1) = openPts_(j,0);
495 DofCoords(cur,2) = openPts_(k,0);
502 for (
int k=0;k<openPts_.dimension(0);k++)
504 for (
int j=0;j<closedPts_.dimension(0);j++)
506 for (
int i=0;i<openPts_.dimension(0);i++)
508 DofCoords(cur,0) = openPts_(i,0);
509 DofCoords(cur,1) = closedPts_(j,0);
510 DofCoords(cur,2) = openPts_(k,0);
517 for (
int k=0;k<closedPts_.dimension(0);k++)
519 for (
int j=0;j<openPts_.dimension(0);j++)
521 for (
int i=0;i<openPts_.dimension(0);i++)
523 DofCoords(cur,0) = openPts_(i,0);
524 DofCoords(cur,1) = openPts_(j,0);
525 DofCoords(cur,2) = closedPts_(k,0);
EBasis basisType_
Type of the basis.
Basis_HDIV_HEX_In_FEM(int order, const ArrayScalar &ptsClosed, const ArrayScalar &ptsOpen)
Constructor.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
virtual void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
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)
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedral 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. ...