54 const EPointType pointType ):
56 coeffs( (n+1)*(n+2) , n*(n+2) )
58 const int N = n*(n+2);
61 this ->
basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() );
67 const int littleN = n*(n+1);
68 const int bigN = (n+1)*(n+2);
69 const int scalarSmallestN = (n-1)*n / 2;
70 const int scalarLittleN = littleN/2;
71 const int scalarBigN = bigN/2;
77 Teuchos::SerialDenseMatrix<int,Scalar> V1(bigN, N);
87 for (
int i=0;i<scalarLittleN;i++) {
89 V1(scalarBigN+i,scalarLittleN+i) = 1.0;
94 CubatureDirectTriDefault<Scalar,FieldContainer<Scalar> > myCub( 2 * n );
95 FieldContainer<Scalar> cubPoints( myCub.getNumPoints() , 2 );
96 FieldContainer<Scalar> cubWeights( myCub.getNumPoints() );
97 myCub.getCubature( cubPoints , cubWeights );
100 FieldContainer<Scalar> phisAtCubPoints( scalarBigN , myCub.getNumPoints() );
101 Phis.getValues( phisAtCubPoints , cubPoints , OPERATOR_VALUE );
104 for (
int i=0;i<n;i++) {
105 for (
int j=0;j<scalarBigN;j++) {
106 V1(j,littleN+i) = 0.0;
107 for (
int k=0;k<myCub.getNumPoints();k++) {
109 cubWeights(k) * cubPoints(k,0)
110 * phisAtCubPoints(scalarSmallestN+i,k)
111 * phisAtCubPoints(j,k);
114 for (
int j=0;j<scalarBigN;j++) {
115 V1(j+scalarBigN,littleN+i) = 0.0;
116 for (
int k=0;k<myCub.getNumPoints();k++) {
117 V1(j+scalarBigN,littleN+i) +=
118 cubWeights(k) * cubPoints(k,1)
119 * phisAtCubPoints(scalarSmallestN+i,k)
120 * phisAtCubPoints(j,k);
129 Teuchos::SerialDenseMatrix<int,Scalar> V2(N , bigN);
133 FieldContainer<Scalar> linePts( n , 1 );
138 shards::CellTopology linetop(shards::getCellTopologyData<shards::Line<2> >() );
140 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( linePts ,
145 FieldContainer<Scalar> edgePts( n , 2 );
146 FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , n );
149 const Scalar nx[] = {0.0,1.0,-1.0};
150 const Scalar ny[] = {-1.0,1.0,0.0};
152 for (
int i=0;i<3;i++) {
159 Phis.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
162 for (
int j=0;j<n;j++) {
164 for (
int k=0;k<scalarBigN;k++) {
165 V2(n*i+j,k) = nx[i] * phisAtEdgePoints(k,j);
166 V2(n*i+j,k+scalarBigN) = ny[i] * phisAtEdgePoints(k,j);
184 if (numInternalPoints > 0) {
185 FieldContainer<Scalar> internalPoints( numInternalPoints , 2 );
186 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
192 FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
193 Phis.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
196 for (
int i=0;i<numInternalPoints;i++) {
197 for (
int j=0;j<scalarBigN;j++) {
199 V2(3*n+i,j) = phisAtInternalPoints(j,i);
201 V2(3*n+numInternalPoints+i,scalarBigN+j) = phisAtInternalPoints(j,i);
209 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
212 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
218 Teuchos::SerialDenseSolver<int,Scalar> solver;
219 solver.setMatrix( rcp( &Vsdm ,
false ) );
222 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
223 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
227 for (
int i=0;i<bigN;i++) {
228 for (
int j=0;j<N;j++) {
286 const ArrayScalar & inputPoints,
287 const EOperator operatorType)
const {
290#ifdef HAVE_INTREPID_DEBUG
291 Intrepid::getValues_HDIV_Args<Scalar, ArrayScalar>(outputValues,
294 this -> getBaseCellTopology(),
295 this -> getCardinality() );
297 const int numPts = inputPoints.dimension(0);
298 const int deg =
this -> getDegree();
299 const int scalarBigN = (deg+1)*(deg+2)/2;
302 switch (operatorType) {
305 FieldContainer<Scalar> phisCur( scalarBigN , numPts );
306 Phis.getValues( phisCur , inputPoints , OPERATOR_VALUE );
308 for (
int i=0;i<outputValues.dimension(0);i++) {
309 for (
int j=0;j<outputValues.dimension(1);j++) {
310 outputValues(i,j,0) = 0.0;
311 outputValues(i,j,1) = 0.0;
312 for (
int k=0;k<scalarBigN;k++) {
313 outputValues(i,j,0) += coeffs(k,i) * phisCur(k,j);
314 outputValues(i,j,1) += coeffs(k+scalarBigN,i) * phisCur(k,j);
322 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 2 );
323 Phis.getValues( phisCur , inputPoints , OPERATOR_GRAD );
324 for (
int i=0;i<outputValues.dimension(0);i++) {
325 for (
int j=0;j<outputValues.dimension(1);j++) {
327 outputValues(i,j) = 0.0;
328 for (
int k=0;k<scalarBigN;k++) {
329 outputValues(i,j) += coeffs(k,i) * phisCur(k,j,0);
332 for (
int k=0;k<scalarBigN;k++) {
333 outputValues(i,j) += coeffs(k+scalarBigN,i) * phisCur(k,j,1);
340 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
341 ">>> ERROR (Basis_HDIV_TRI_In_FEM): Operator type not implemented");
345 catch (std::invalid_argument &exception){
346 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
347 ">>> ERROR (Basis_HDIV_TRI_In_FEM): Operator type not implemented");