52 template<
class Scalar,
class ArrayScalar>
56 coeffs_( (n+1)*(n+2)*(n+3)/2 , n*(n+1)*(n+3)/2 )
58 const int N = n*(n+1)*(n+3)/2;
62 = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >() );
68 const int littleN = n*(n+1)*(n+2)/2;
69 const int bigN = (n+1)*(n+2)*(n+3)/2;
70 const int start_PkH = (n-1)*n*(n+1)/6;
71 const int dim_PkH = n*(n+1)*(n+2)/6 - start_PkH;
72 const int scalarLittleN = littleN/3;
73 const int scalarBigN = bigN/3;
79 Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, N);
91 for (
int i=0;i<scalarLittleN;i++) {
92 for (
int k=0;k<3;k++) {
93 V1(i+k*scalarBigN,i+k*scalarLittleN) = 1.0;
106 Phis_.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
110 for (
int i=0;i<dim_PkH;i++) {
111 for (
int j=0;j<scalarBigN;j++) {
112 V1(j,littleN+i) = 0.0;
113 for (
int d=0;d<3;d++) {
115 V1(j+d*scalarBigN,littleN+i) +=
116 cubWeights(k) * cubPoints(k,d)
117 * phisAtCubPoints(start_PkH+i,k)
118 * phisAtCubPoints(j,k);
126 Teuchos::SerialDenseMatrix<int,Scalar> V2(N , bigN);
128 shards::CellTopology faceTop(shards::getCellTopologyData<shards::Triangle<3> >() );
134 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts ,
151 Scalar normal[][4] = { {0.0,0.5,-0.5,0.0},
153 {0.0,0.5,0.0,-0.5} };
155 for (
int i=0;i<4;i++) {
162 Phis_.getValues( phisAtFacePoints , facePts , OPERATOR_VALUE );
165 for (
int j=0;j<numPtsPerFace;j++) {
167 for (
int k=0;k<scalarBigN;k++) {
168 for (
int l=0;l<3;l++) {
169 V2(numPtsPerFace*i+j,k+l*scalarBigN) = normal[l][i] * phisAtFacePoints(k,j);
185 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
192 Phis_.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
195 for (
int i=0;i<numInternalPoints;i++) {
196 for (
int j=0;j<scalarBigN;j++) {
197 for (
int k=0;k<3;k++) {
198 V2(4*numPtsPerFace+k*numInternalPoints+i,k*scalarBigN+j) = phisAtInternalPoints(j,i);
204 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
207 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
213 Teuchos::SerialDenseSolver<int,Scalar> solver;
214 solver.setMatrix( rcp( &Vsdm ,
false ) );
218 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
219 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
223 for (
int i=0;i<bigN;i++) {
224 for (
int j=0;j<N;j++) {
230 template<
class Scalar,
class ArrayScalar>
241 int *tags =
new int[ tagSize * this->getCardinality() ];
243 const int deg = this->getDegree();
245 const int numPtsPerFace = deg*(deg+1)/2;
248 for (
int f=0;f<4;f++) {
249 for (
int i=0;i<numPtsPerFace;i++) {
250 tag_cur[0] = 2; tag_cur[1] = f; tag_cur[2] = i; tag_cur[3] = numPtsPerFace;
257 const int numInternalDof = this->getCardinality() - 4 * numPtsPerFace;
258 int internalDofCur=0;
259 for (
int i=4*numPtsPerFace;i<this->getCardinality();i++) {
260 tag_cur[0] = 3; tag_cur[1] = 0; tag_cur[2] = internalDofCur; tag_cur[3] = numInternalDof;
267 this -> ordinalToTag_,
269 this -> basisCardinality_,
281 template<
class Scalar,
class ArrayScalar>
283 const ArrayScalar & inputPoints,
287 #ifdef HAVE_INTREPID_DEBUG 288 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
291 this -> getBaseCellTopology(),
292 this -> getCardinality() );
294 const int numPts = inputPoints.dimension(0);
295 const int deg =
this -> getDegree();
296 const int scalarBigN = (deg+1)*(deg+2)*(deg+3)/6;
299 switch (operatorType) {
303 Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
305 for (
int i=0;i<outputValues.dimension(0);i++) {
306 for (
int j=0;j<outputValues.dimension(1);j++) {
307 for (
int l=0;l<3;l++) {
308 outputValues(i,j,l) = 0.0;
310 for (
int k=0;k<scalarBigN;k++) {
311 for (
int l=0;l<3;l++) {
312 outputValues(i,j,l) += coeffs_(k+l*scalarBigN,i) * phisCur(k,j);
322 Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
323 for (
int i=0;i<outputValues.dimension(0);i++) {
324 for (
int j=0;j<outputValues.dimension(1);j++) {
325 outputValues(i,j) = 0.0;
326 for (
int k=0;k<scalarBigN;k++) {
327 outputValues(i,j) += coeffs_(k,i) * phisCur(k,j,0);
328 outputValues(i,j) += coeffs_(k+scalarBigN,i) * phisCur(k,j,1);
329 outputValues(i,j) += coeffs_(k+2*scalarBigN,i) * phisCur(k,j,2);
336 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
337 ">>> ERROR (Basis_HDIV_TET_In_FEM): Operator type not implemented");
341 catch (std::invalid_argument &exception){
342 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
343 ">>> ERROR (Basis_HDIV_TET_In_FEM): Operator type not implemented");
350 template<
class Scalar,
class ArrayScalar>
352 const ArrayScalar & inputPoints,
353 const ArrayScalar & cellVertices,
355 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
356 ">>> ERROR (Basis_HDIV_TET_In_FEM): FEM Basis calling an FVD member function");
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
Evaluation of a FEM basis on a reference Tetrahedron cell.
virtual void initializeTags()
Initializes tagToOrdinal_ and ordinalToTag_ lookup arrays.
FieldContainer< Scalar > coeffs_
expansion coefficients of the nodal basis in terms of the orthgonal one
EBasis basisType_
Type of the basis.
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.
Defines direct integration rules on a tetrahedron.
bool basisTagsAreSet_
"true" if tagToOrdinal_ and ordinalToTag_ have been initialized
EOperator
Enumeration of primitive operators available in Intrepid. Primitive operators act on reconstructed fu...
ECoordinates basisCoordinates_
The coordinate system for which the basis is defined.
Basis_HDIV_TET_In_FEM(const int n, const EPointType pointType)
Constructor.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
shards::CellTopology basisCellTopology_
Base topology of the cells for which the basis is defined. See the Shards package http://trilinos...
Basis_HGRAD_TET_Cn_FEM_ORTH< Scalar, FieldContainer< Scalar > > Phis_
Orthogonal basis out of which the nodal basis is constructed.
virtual int getNumPoints() const
Returns the number of cubature points.
EPointType
Enumeration of types of point distributions in Intrepid.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...
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. ...