1 #ifndef INTREPID_HGRAD_TRI_CN_FEMDEF_HPP 2 #define INTREPID_HGRAD_TRI_CN_FEMDEF_HPP 53 template<
class Scalar,
class ArrayScalar>
57 V((n+1)*(n+2)/2,(n+1)*(n+2)/2),
58 Vinv((n+1)*(n+2)/2,(n+1)*(n+2)/2),
59 latticePts( (n+1)*(n+2)/2 , 2 )
61 TEUCHOS_TEST_FOR_EXCEPTION( n <= 0, std::invalid_argument, "polynomial order must be >= 1
"); 63 const int N = (n+1)*(n+2)/2; 64 this -> basisCardinality_ = N; 65 this -> basisDegree_ = n; 66 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() ); 67 this -> basisType_ = BASIS_FEM_FIAT; 68 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 69 this -> basisTagsAreSet_ = false; 73 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() ); 75 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts , 82 // form Vandermonde matrix. Actually, this is the transpose of the VDM, 83 // so we transpose on copy below. 85 Phis.getValues( V , latticePts , OPERATOR_VALUE ); 87 // now I need to copy V into a Teuchos array to do the inversion 88 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N); 89 for (int i=0;i<N;i++) { 90 for (int j=0;j<N;j++) { 96 Teuchos::SerialDenseSolver<int,Scalar> solver; 97 solver.setMatrix( rcp( &Vsdm , false ) ); 100 // now I need to copy the inverse into Vinv 101 for (int i=0;i<N;i++) { 102 for (int j=0;j<N;j++) { 103 Vinv(i,j) = Vsdm(j,i); 110 template<class Scalar, class ArrayScalar> 111 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::initializeTags() { 113 // Basis-dependent initializations 114 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 115 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 116 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 117 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 119 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 121 int *tags = new int[ tagSize * this->getCardinality() ]; 123 const int degree = this->getDegree(); 125 // BEGIN DOF ALONG BOTTOM EDGE 127 // the first dof is on vertex 0 128 tag_cur[0] = 0; tag_cur[1] = 0; tag_cur[2] = 0; tag_cur[3] = 1; 131 // next degree-1 dof are on edge 0 132 for (int i=1;i<degree;i++) { 133 tag_cur[0] = 1; tag_cur[1] = 0; tag_cur[2] = i-1; tag_cur[3] = degree-1; 137 // last dof is on vertex 1 138 tag_cur[0] = 0; tag_cur[1] = 1; tag_cur[2] = 0; tag_cur[3] = 1; 141 // END DOF ALONG BOTTOM EDGE 143 int num_internal_dof = PointTools::getLatticeSize( this->getBaseCellTopology() , 147 int internal_dof_cur = 0; 149 // BEGIN DOF ALONG INTERNAL HORIZONTAL LINES 150 for (int i=1;i<degree;i++) { 151 // first dof along line is on edge #2 152 tag_cur[0] = 1; tag_cur[1] = 2; tag_cur[2] = i-1; tag_cur[3] = degree-1; 155 // next dof are internal 156 for (int j=1;j<degree-i;j++) { 157 tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = internal_dof_cur; tag_cur[3] = num_internal_dof; 162 // last dof along line is on edge 1 163 tag_cur[0] = 1; tag_cur[1] = 1; tag_cur[2] = i-1; tag_cur[3] = degree-1; 167 // END DOF ALONG INTERNAL HORIZONTAL LINES 169 // LAST DOF IS ON VERTEX 2 170 tag_cur[0] = 0; tag_cur[1] = 2; tag_cur[2] = 0; tag_cur[3] = 1; 174 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 175 this -> ordinalToTag_, 177 this -> basisCardinality_, 189 template<class Scalar, class ArrayScalar> 190 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 191 const ArrayScalar & inputPoints, 192 const EOperator operatorType) const { 195 #ifdef HAVE_INTREPID_DEBUG 196 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 199 this -> getBaseCellTopology(), 200 this -> getCardinality() ); 202 const int numPts = inputPoints.dimension(0); 203 const int numBf = this->getCardinality(); 206 switch (operatorType) { 209 FieldContainer<Scalar> phisCur( numBf , numPts ); 210 Phis.getValues( phisCur , inputPoints , operatorType ); 211 for (int i=0;i<outputValues.dimension(0);i++) { 212 for (int j=0;j<outputValues.dimension(1);j++) { 213 outputValues(i,j) = 0.0; 214 for (int k=0;k<this->getCardinality();k++) { 215 outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j); 234 (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,2): getDkCardinality(operatorType,2); 236 FieldContainer<Scalar> phisCur( numBf , numPts , dkcard ); 237 Phis.getValues( phisCur , inputPoints , operatorType ); 239 for (int i=0;i<outputValues.dimension(0);i++) { 240 for (int j=0;j<outputValues.dimension(1);j++) { 241 for (int k=0;k<outputValues.dimension(2);k++) { 242 outputValues(i,j,k) = 0.0; 243 for (int l=0;l<this->getCardinality();l++) { 244 outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k); 251 case OPERATOR_CURL: // only works in 2d. first component is -d/dy, second is d/dx 253 FieldContainer<Scalar> phisCur( numBf , numPts , getDkCardinality( OPERATOR_D1 , 2 ) ); 254 Phis.getValues( phisCur , inputPoints , OPERATOR_D1 ); 256 for (int i=0;i<outputValues.dimension(0);i++) { 257 for (int j=0;j<outputValues.dimension(1);j++) { 258 outputValues(i,j,0) = 0.0; 259 outputValues(i,j,1) = 0.0; 260 for (int k=0;k<this->getCardinality();k++) { 261 outputValues(i,j,0) += this->Vinv(k,i) * phisCur(k,j,1); 263 for (int k=0;k<this->getCardinality();k++) { 264 outputValues(i,j,1) -= this->Vinv(k,i) * phisCur(k,j,0); 271 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument, 276 catch (std::invalid_argument &exception){ 277 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument, 285 template<class Scalar, class ArrayScalar> 286 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 287 const ArrayScalar & inputPoints, 288 const ArrayScalar & cellVertices, 289 const EOperator operatorType) const { 290 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 295 }// namespace Intrepid Basis_HGRAD_TRI_Cn_FEM(const int n, const EPointType pointType)
Constructor.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell...
EPointType
Enumeration of types of point distributions in Intrepid.