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> >() );
66 const int littleN = n*(n+1);
67 const int bigN = (n+1)*(n+2);
68 const int scalarSmallestN = (n-1)*n / 2;
69 const int scalarLittleN = littleN/2;
70 const int scalarBigN = bigN/2;
76 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,1)
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,0)
119 * phisAtCubPoints(scalarSmallestN+i,k)
120 * phisAtCubPoints(j,k);
129 Teuchos::SerialDenseMatrix<int,Scalar> V2(N , bigN);
133 FieldContainer<Scalar> linePts( n , 1 );
134 if (pointType == POINTTYPE_WARPBLEND) {
135 CubatureDirectLineGauss<Scalar> edgeRule( 2*n - 1 );
136 FieldContainer<Scalar> edgeCubWts( n );
137 edgeRule.getCubature( linePts , edgeCubWts );
139 else if (pointType == POINTTYPE_EQUISPACED ) {
140 shards::CellTopology linetop(shards::getCellTopologyData<shards::Line<2> >() );
142 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( linePts ,
145 POINTTYPE_EQUISPACED );
149 FieldContainer<Scalar> edgePts( n , 2 );
150 FieldContainer<Scalar> phisAtEdgePoints( scalarBigN , n );
151 FieldContainer<Scalar> edgeTan(2);
153 for (
int i=0;i<3;i++) {
158 for (
int j=0;j<2;j++) {
168 Phis_.getValues( phisAtEdgePoints , edgePts , OPERATOR_VALUE );
171 for (
int j=0;j<n;j++) {
173 for (
int k=0;k<scalarBigN;k++) {
174 V2(n*i+j,k) = edgeTan(0) * phisAtEdgePoints(k,j);
175 V2(n*i+j,k+scalarBigN) = edgeTan(1) * phisAtEdgePoints(k,j);
187 if (numInternalPoints > 0) {
188 FieldContainer<Scalar> internalPoints( numInternalPoints , 2 );
189 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( internalPoints ,
195 FieldContainer<Scalar> phisAtInternalPoints( scalarBigN , numInternalPoints );
196 Phis_.getValues( phisAtInternalPoints , internalPoints , OPERATOR_VALUE );
199 for (
int i=0;i<numInternalPoints;i++) {
200 for (
int j=0;j<scalarBigN;j++) {
202 V2(3*n+i,j) = phisAtInternalPoints(j,i);
204 V2(3*n+numInternalPoints+i,scalarBigN+j) = phisAtInternalPoints(j,i);
212 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm( N , N );
215 Vsdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V2 , V1 , 0.0 );
221 Teuchos::SerialDenseSolver<int,Scalar> solver;
222 solver.setMatrix( rcp( &Vsdm ,
false ) );
225 Teuchos::SerialDenseMatrix<int,Scalar> Csdm( bigN , N );
226 Csdm.multiply( Teuchos::NO_TRANS , Teuchos::NO_TRANS , 1.0 , V1 , Vsdm , 0.0 );
230 for (
int i=0;i<bigN;i++) {
231 for (
int j=0;j<N;j++) {
289 const ArrayScalar & inputPoints,
290 const EOperator operatorType)
const {
293#ifdef HAVE_INTREPID_DEBUG
294 Intrepid::getValues_HCURL_Args<Scalar, ArrayScalar>(outputValues,
297 this -> getBaseCellTopology(),
298 this -> getCardinality() );
300 const int numPts = inputPoints.dimension(0);
301 const int deg =
this -> getDegree();
302 const int scalarBigN = (deg+1)*(deg+2)/2;
305 switch (operatorType) {
308 FieldContainer<Scalar> phisCur( scalarBigN , numPts );
309 Phis_.getValues( phisCur , inputPoints , OPERATOR_VALUE );
311 for (
int i=0;i<outputValues.dimension(0);i++) {
312 for (
int j=0;j<outputValues.dimension(1);j++) {
313 outputValues(i,j,0) = 0.0;
314 outputValues(i,j,1) = 0.0;
315 for (
int k=0;k<scalarBigN;k++) {
316 outputValues(i,j,0) += coeffs_(k,i) * phisCur(k,j);
317 outputValues(i,j,1) += coeffs_(k+scalarBigN,i) * phisCur(k,j);
325 FieldContainer<Scalar> phisCur( scalarBigN , numPts , 2 );
326 Phis_.getValues( phisCur , inputPoints , OPERATOR_GRAD );
327 for (
int i=0;i<outputValues.dimension(0);i++) {
328 for (
int j=0;j<outputValues.dimension(1);j++) {
330 outputValues(i,j) = 0.0;
331 for (
int k=0;k<scalarBigN;k++) {
332 outputValues(i,j) -= coeffs_(k,i) * phisCur(k,j,1);
335 for (
int k=0;k<scalarBigN;k++) {
336 outputValues(i,j) += coeffs_(k+scalarBigN,i) * phisCur(k,j,0);
343 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
344 ">>> ERROR (Basis_HCURL_TRI_In_FEM): Operator type not implemented");
348 catch (std::invalid_argument &exception){
349 TEUCHOS_TEST_FOR_EXCEPTION(
true , std::invalid_argument,
350 ">>> ERROR (Basis_HCURL_TRI_In_FEM): Operator type not implemented");