54 const EPointType pointType ):
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;
99 CubatureDirectTetDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n );
100 FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 3 );
101 FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
102 myCub.getCubature( cubPoints , cubWeights );
105 FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
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++) {
114 for (
int k=0;k<myCub.getNumPoints();k++) {
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> >() );
133 FieldContainer<Scalar> twoDPts( numPtsPerFace , 2 );
134 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( twoDPts ,
141 FieldContainer<Scalar> facePts( numPtsPerFace , 3 );
142 FieldContainer<Scalar> phisAtFacePoints( scalarBigN ,
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);
184 FieldContainer<Scalar> internalPoints( numInternalPoints , 3 );
185 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
191 FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
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++) {
283 const ArrayScalar & inputPoints,
284 const EOperator operatorType)
const {
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) {
302 FieldContainer<Scalar> phisCur( scalarBigN , numPts );
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);
321 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 3 );
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");
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.