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 ),
62 * closedBasis_.getCardinality() * openBasis_.getCardinality();
63 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
68 Array<Array<RCP<Basis<Scalar,ArrayScalar > > > > bases(3);
69 bases[0].resize(3); bases[1].resize(3); bases[2].resize(3);
70 bases[0][0] = rcp( &openBasis_ ,
false );
71 bases[0][1] = rcp( &closedBasis_ ,
false );
72 bases[0][2] = rcp( &closedBasis_ ,
false );
73 bases[1][0] = rcp( &closedBasis_ ,
false );
74 bases[1][1] = rcp( &openBasis_ ,
false );
75 bases[1][2] = rcp( &closedBasis_ ,
false );
76 bases[2][0] = rcp( &closedBasis_ ,
false );
77 bases[2][1] = rcp( &closedBasis_ ,
false );
78 bases[2][2] = rcp( &openBasis_ ,
false );
79 this->setBases( bases );
83 template<
class Scalar,
class ArrayScalar>
85 closedBasis_( order , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL:POINTTYPE_EQUISPACED ),
86 openBasis_( order-1 , pointType==POINTTYPE_SPECTRAL?POINTTYPE_SPECTRAL_OPEN:POINTTYPE_EQUISPACED ),
87 closedPts_( order+1 , 1 ),
92 * closedBasis_.getCardinality() * openBasis_.getCardinality();
93 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >() );
98 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( closedPts_ ,
99 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
102 pointType==POINTTYPE_SPECTRAL?POINTTYPE_WARPBLEND:POINTTYPE_EQUISPACED );
104 if (pointType == POINTTYPE_SPECTRAL)
106 PointTools::getGaussPoints<Scalar,FieldContainer<Scalar> >( openPts_ ,
111 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( openPts_ ,
112 shards::CellTopology(shards::getCellTopologyData<shards::Line<2> >()) ,
115 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( &openBasis_ ,
false );
121 bases[0][1] = rcp( &closedBasis_ ,
false );
122 bases[0][2] = rcp( &closedBasis_ ,
false );
123 bases[1][0] = rcp( &closedBasis_ ,
false );
124 bases[1][1] = rcp( &openBasis_ ,
false );
125 bases[1][2] = rcp( &closedBasis_ ,
false );
126 bases[2][0] = rcp( &closedBasis_ ,
false );
127 bases[2][1] = rcp( &closedBasis_ ,
false );
128 bases[2][2] = rcp( &openBasis_ ,
false );
129 this->setBases( bases );
133 template<
class Scalar,
class ArrayScalar>
142 std::vector<int> tags( tagSize * this->getCardinality() );
144 const std::vector<std::vector<int> >& closedDofTags = closedBasis_.getAllDofTags();
145 const std::vector<std::vector<int> >& openDofTags = openBasis_.getAllDofTags();
147 std::map<int,std::map<int,int> > total_dof_per_entity;
148 std::map<int,std::map<int,int> > current_dof_per_entity;
150 const int order = this->basisDegree_;
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;
172 for (
int i=0;i<12;i++) {
173 total_dof_per_entity[1][i] = order;
176 for (
int i=0;i<6;i++) {
177 total_dof_per_entity[2][i] = 2 * (order - 1) * order;
180 total_dof_per_entity[3][0] = this->basisCardinality_ - 12 * order - 6 * total_dof_per_entity[2][0];
183 for (
int k=0;k<closedBasis_.getCardinality();k++) {
184 const int dimk = closedDofTags[k][0];
185 const int entk = closedDofTags[k][1];
186 for (
int j=0;j<closedBasis_.getCardinality();j++) {
187 const int dimj = closedDofTags[j][0];
188 const int entj = closedDofTags[j][1];
189 for (
int i=0;i<openBasis_.getCardinality();i++) {
190 const int dimi = openDofTags[i][0];
191 const int enti = openDofTags[i][1];
195 tags[4*tagcur] = dofdim;
196 tags[4*tagcur+1] = dofent;
197 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
198 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
199 current_dof_per_entity[dofdim][dofent]++;
204 for (
int k=0;k<closedBasis_.getCardinality();k++) {
205 const int dimk = closedDofTags[k][0];
206 const int entk = closedDofTags[k][1];
207 for (
int j=0;j<openBasis_.getCardinality();j++) {
208 const int dimj = openDofTags[j][0];
209 const int entj = openDofTags[j][1];
210 for (
int i=0;i<closedBasis_.getCardinality();i++) {
211 const int dimi = closedDofTags[i][0];
212 const int enti = closedDofTags[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 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
220 current_dof_per_entity[dofdim][dofent]++;
225 for (
int k=0;k<openBasis_.getCardinality();k++) {
226 const int dimk = openDofTags[k][0];
227 const int entk = openDofTags[k][1];
228 for (
int j=0;j<closedBasis_.getCardinality();j++) {
229 const int dimj = closedDofTags[j][0];
230 const int entj = closedDofTags[j][1];
231 for (
int i=0;i<closedBasis_.getCardinality();i++) {
232 const int dimi = closedDofTags[i][0];
233 const int enti = closedDofTags[i][1];
237 tags[4*tagcur] = dofdim;
238 tags[4*tagcur+1] = dofent;
239 tags[4*tagcur+2] = current_dof_per_entity[dofdim][dofent];
240 tags[4*tagcur+3] = total_dof_per_entity[dofdim][dofent];
241 current_dof_per_entity[dofdim][dofent]++;
248 Intrepid::setOrdinalTagData(
this -> tagToOrdinal_,
249 this -> ordinalToTag_,
251 this -> basisCardinality_,
259 template<
class Scalar,
class ArrayScalar>
261 const ArrayScalar & inputPoints,
262 const EOperator operatorType)
const {
265 #ifdef HAVE_INTREPID_DEBUG
266 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
269 this -> getBaseCellTopology(),
270 this -> getCardinality() );
274 int dim0 = inputPoints.dimension(0);
281 for (
int i=0;i<dim0;i++) {
282 xPoints(i,0) = inputPoints(i,0);
283 yPoints(i,0) = inputPoints(i,1);
284 zPoints(i,0) = inputPoints(i,2);
287 switch (operatorType) {
297 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
298 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
299 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE );
300 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
301 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
302 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE );
306 for (
int k=0;k<closedBasis_.getCardinality();k++) {
307 for (
int j=0;j<closedBasis_.getCardinality();j++) {
308 for (
int i=0;i<openBasis_.getCardinality();i++) {
309 for (
int l=0;l<dim0;l++) {
310 outputValues(bfcur,l,0) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * closedBasisValsZPts(k,l);
311 outputValues(bfcur,l,1) = 0.0;
312 outputValues(bfcur,l,2) = 0.0;
320 for (
int k=0;k<closedBasis_.getCardinality();k++) {
321 for (
int j=0;j<openBasis_.getCardinality();j++) {
322 for (
int i=0;i<closedBasis_.getCardinality();i++) {
323 for (
int l=0;l<dim0;l++) {
324 outputValues(bfcur,l,0) = 0.0;
325 outputValues(bfcur,l,1) = closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l);
326 outputValues(bfcur,l,2) = 0.0;
334 for (
int k=0;k<openBasis_.getCardinality();k++) {
335 for (
int j=0;j<closedBasis_.getCardinality();j++) {
336 for (
int i=0;i<closedBasis_.getCardinality();i++) {
337 for (
int l=0;l<dim0;l++) {
338 outputValues(bfcur,l,0) = 0.0;
339 outputValues(bfcur,l,1) = 0.0;
340 outputValues(bfcur,l,2) = closedBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l);
360 closedBasis_.getValues( closedBasisValsXPts , xPoints , OPERATOR_VALUE );
361 closedBasis_.getValues( closedBasisValsYPts , yPoints , OPERATOR_VALUE );
362 closedBasis_.getValues( closedBasisValsZPts , zPoints , OPERATOR_VALUE );
363 closedBasis_.getValues( closedBasisDerivsXPts , xPoints , OPERATOR_D1 );
364 closedBasis_.getValues( closedBasisDerivsYPts , yPoints , OPERATOR_D1 );
365 closedBasis_.getValues( closedBasisDerivsZPts , zPoints , OPERATOR_D1 );
366 openBasis_.getValues( openBasisValsXPts , xPoints , OPERATOR_VALUE );
367 openBasis_.getValues( openBasisValsYPts , yPoints , OPERATOR_VALUE );
368 openBasis_.getValues( openBasisValsZPts , zPoints , OPERATOR_VALUE );
373 for (
int k=0;k<closedBasis_.getCardinality();k++) {
374 for (
int j=0;j<closedBasis_.getCardinality();j++) {
375 for (
int i=0;i<openBasis_.getCardinality();i++) {
376 for (
int l=0;l<dim0;l++) {
377 outputValues(bfcur,l,0) = 0.0;
378 outputValues(bfcur,l,1) = openBasisValsXPts(i,l) * closedBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0);
379 outputValues(bfcur,l,2) = -openBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * closedBasisValsZPts(k,l);
387 for (
int k=0;k<closedBasis_.getCardinality();k++) {
388 for (
int j=0;j<openBasis_.getCardinality();j++) {
389 for (
int i=0;i<closedBasis_.getCardinality();i++) {
390 for (
int l=0;l<dim0;l++) {
391 outputValues(bfcur,l,0) = -closedBasisValsXPts(i,l) * openBasisValsYPts(j,l) * closedBasisDerivsZPts(k,l,0);
392 outputValues(bfcur,l,1) = 0.0;
393 outputValues(bfcur,l,2) = closedBasisDerivsXPts(i,l,0) * openBasisValsYPts(j,l) * closedBasisValsZPts(k,l);
401 for (
int k=0;k<openBasis_.getCardinality();k++) {
402 for (
int j=0;j<closedBasis_.getCardinality();j++) {
403 for (
int i=0;i<closedBasis_.getCardinality();i++) {
404 for (
int l=0;l<dim0;l++) {
405 outputValues(bfcur,l,0) = closedBasisValsXPts(i,l) * closedBasisDerivsYPts(j,l,0) * openBasisValsZPts(k,l);
406 outputValues(bfcur,l,1) = -closedBasisDerivsXPts(i,l,0) * closedBasisValsYPts(j,l) * openBasisValsZPts(k,l);
407 outputValues(bfcur,l,2) = 0.0;
416 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument,
417 ">>> ERROR (Basis_HCURL_HEX_In_FEM): DIV is invalid operator for HCURL Basis Functions");
421 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_GRAD), std::invalid_argument,
422 ">>> ERROR (Basis_HCURL_HEX_In_FEM): GRAD is invalid operator for HCURL Basis Functions");
435 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType == OPERATOR_D1) ||
436 (operatorType == OPERATOR_D2) ||
437 (operatorType == OPERATOR_D3) ||
438 (operatorType == OPERATOR_D4) ||
439 (operatorType == OPERATOR_D5) ||
440 (operatorType == OPERATOR_D6) ||
441 (operatorType == OPERATOR_D7) ||
442 (operatorType == OPERATOR_D8) ||
443 (operatorType == OPERATOR_D9) ||
444 (operatorType == OPERATOR_D10) ),
445 std::invalid_argument,
446 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Invalid operator type");
450 TEUCHOS_TEST_FOR_EXCEPTION( ( (operatorType != OPERATOR_VALUE) &&
451 (operatorType != OPERATOR_GRAD) &&
452 (operatorType != OPERATOR_CURL) &&
453 (operatorType != OPERATOR_CURL) &&
454 (operatorType != OPERATOR_D1) &&
455 (operatorType != OPERATOR_D2) &&
456 (operatorType != OPERATOR_D3) &&
457 (operatorType != OPERATOR_D4) &&
458 (operatorType != OPERATOR_D5) &&
459 (operatorType != OPERATOR_D6) &&
460 (operatorType != OPERATOR_D7) &&
461 (operatorType != OPERATOR_D8) &&
462 (operatorType != OPERATOR_D9) &&
463 (operatorType != OPERATOR_D10) ),
464 std::invalid_argument,
465 ">>> ERROR (Basis_HCURL_HEX_In_FEM): Invalid operator type");
471 template<
class Scalar,
class ArrayScalar>
473 const ArrayScalar & inputPoints,
474 const ArrayScalar & cellVertices,
475 const EOperator operatorType)
const {
476 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
477 ">>> ERROR (Basis_HCURL_HEX_In_FEM): FEM Basis calling an FVD member function");
480 template<
class Scalar,
class ArrayScalar>
486 for (
int k=0;k<closedPts_.dimension(0);k++)
488 for (
int j=0;j<closedPts_.dimension(0);j++)
490 for (
int i=0;i<openPts_.dimension(0);i++)
492 DofCoords(cur,0) = openPts_(i,0);
493 DofCoords(cur,1) = closedPts_(j,0);
494 DofCoords(cur,2) = closedPts_(k,0);
500 for (
int k=0;k<closedPts_.dimension(0);k++)
502 for (
int j=0;j<openPts_.dimension(0);j++)
504 for (
int i=0;i<closedPts_.dimension(0);i++)
506 DofCoords(cur,0) = closedPts_(i,0);
507 DofCoords(cur,1) = openPts_(j,0);
508 DofCoords(cur,2) = closedPts_(k,0);
515 for (
int k=0;k<openPts_.dimension(0);k++)
517 for (
int j=0;j<closedPts_.dimension(0);j++)
519 for (
int i=0;i<closedPts_.dimension(0);i++)
521 DofCoords(cur,0) = closedPts_(i,0);
522 DofCoords(cur,1) = closedPts_(j,0);
523 DofCoords(cur,2) = openPts_(k,0);
virtual void getDofCoords(ArrayScalar &DofCoords) const
Returns spatial locations (coordinates) of degrees of freedom on a reference cell; defined for interp...
EBasis basisType_
Type of the basis.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Hexahedral cell.
Basis_HCURL_HEX_In_FEM(int order, const ArrayScalar &ptsClosed, const ArrayScalar &ptsOpen)
Constructor.
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)
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. ...
void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.