|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00035 namespace Intrepid { 00036 00037 template<class Scalar, class ArrayScalar> 00038 Basis_HGRAD_TRI_Cn_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TRI_Cn_FEM( const int n , 00039 const EPointType pointType ): 00040 Phis( n ), 00041 V((n+1)*(n+2)/2,(n+1)*(n+2)/2), 00042 Vinv((n+1)*(n+2)/2,(n+1)*(n+2)/2), 00043 latticePts( (n+1)*(n+2)/2 , 2 ) 00044 { 00045 const int N = (n+1)*(n+2)/2; 00046 this -> basisCardinality_ = N; 00047 this -> basisDegree_ = n; 00048 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() ); 00049 this -> basisType_ = BASIS_FEM_FIAT; 00050 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00051 this -> basisTagsAreSet_ = false; 00052 00053 // construct lattice 00054 00055 shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() ); 00056 PointTools::getLattice<Scalar,FieldContainer<Scalar> >( latticePts , 00057 myTri_3 , 00058 n , 00059 0 , 00060 pointType ); 00061 00062 00063 // form Vandermonde matrix. Actually, this is the transpose of the VDM, 00064 // so we transpose on copy below. 00065 00066 Phis.getValues( V , latticePts , OPERATOR_VALUE ); 00067 00068 // now I need to copy V into a Teuchos array to do the inversion 00069 Teuchos::SerialDenseMatrix<int,Scalar> Vsdm(N,N); 00070 for (int i=0;i<N;i++) { 00071 for (int j=0;j<N;j++) { 00072 Vsdm(i,j) = V(i,j); 00073 } 00074 } 00075 00076 // invert the matrix 00077 Teuchos::SerialDenseSolver<int,Scalar> solver; 00078 solver.setMatrix( rcp( &Vsdm , false ) ); 00079 solver.invert( ); 00080 00081 // now I need to copy the inverse into Vinv 00082 for (int i=0;i<N;i++) { 00083 for (int j=0;j<N;j++) { 00084 Vinv(i,j) = Vsdm(j,i); 00085 } 00086 } 00087 00088 } 00089 00090 00091 template<class Scalar, class ArrayScalar> 00092 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::initializeTags() { 00093 00094 // Basis-dependent initializations 00095 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00096 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00097 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00098 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00099 00100 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00101 00102 int *tags = new int[ tagSize * this->getCardinality() ]; 00103 int *tag_cur = tags; 00104 const int degree = this->getDegree(); 00105 00106 // BEGIN DOF ALONG BOTTOM EDGE 00107 00108 // the first dof is on vertex 0 00109 tag_cur[0] = 0; tag_cur[1] = 0; tag_cur[2] = 0; tag_cur[3] = 1; 00110 tag_cur += tagSize; 00111 00112 // next degree-1 dof are on edge 0 00113 for (int i=1;i<degree;i++) { 00114 tag_cur[0] = 1; tag_cur[1] = 0; tag_cur[2] = i-1; tag_cur[3] = degree-1; 00115 tag_cur += tagSize; 00116 } 00117 00118 // last dof is on vertex 1 00119 tag_cur[0] = 0; tag_cur[1] = 1; tag_cur[2] = 0; tag_cur[3] = 1; 00120 tag_cur += tagSize; 00121 00122 // END DOF ALONG BOTTOM EDGE 00123 00124 int num_internal_dof = PointTools::getLatticeSize( this->getBaseCellTopology() , 00125 this->getDegree() , 00126 1 ); 00127 00128 int internal_dof_cur = 0; 00129 00130 // BEGIN DOF ALONG INTERNAL HORIZONTAL LINES 00131 for (int i=1;i<degree;i++) { 00132 // first dof along line is on edge #2 00133 tag_cur[0] = 1; tag_cur[1] = 2; tag_cur[2] = i-1; tag_cur[3] = degree-1; 00134 tag_cur += tagSize; 00135 00136 // next dof are internal 00137 for (int j=1;j<degree-i;j++) { 00138 tag_cur[0] = 2; tag_cur[1] = 0; tag_cur[2] = internal_dof_cur; tag_cur[3] = num_internal_dof; 00139 internal_dof_cur++; 00140 tag_cur += tagSize; 00141 } 00142 00143 // last dof along line is on edge 1 00144 tag_cur[0] = 1; tag_cur[1] = 1; tag_cur[2] = i-1; tag_cur[3] = degree-1; 00145 tag_cur += tagSize; 00146 00147 } 00148 // END DOF ALONG INTERNAL HORIZONTAL LINES 00149 00150 // LAST DOF IS ON VERTEX 2 00151 tag_cur[0] = 0; tag_cur[1] = 2; tag_cur[2] = 0; tag_cur[3] = 1; 00152 // END LAST DOF 00153 00154 00155 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00156 this -> ordinalToTag_, 00157 tags, 00158 this -> basisCardinality_, 00159 tagSize, 00160 posScDim, 00161 posScOrd, 00162 posDfOrd); 00163 00164 delete []tags; 00165 00166 } 00167 00168 00169 00170 template<class Scalar, class ArrayScalar> 00171 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00172 const ArrayScalar & inputPoints, 00173 const EOperator operatorType) const { 00174 00175 // Verify arguments 00176 #ifdef HAVE_INTREPID_DEBUG 00177 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00178 inputPoints, 00179 operatorType, 00180 this -> getBaseCellTopology(), 00181 this -> getCardinality() ); 00182 #endif 00183 const int numPts = inputPoints.dimension(0); 00184 const int numBf = this->getCardinality(); 00185 00186 try { 00187 switch (operatorType) { 00188 case OPERATOR_VALUE: 00189 { 00190 FieldContainer<Scalar> phisCur( numBf , numPts ); 00191 Phis.getValues( phisCur , inputPoints , operatorType ); 00192 for (int i=0;i<outputValues.dimension(0);i++) { 00193 for (int j=0;j<outputValues.dimension(1);j++) { 00194 outputValues(i,j) = 0.0; 00195 for (int k=0;k<this->getCardinality();k++) { 00196 outputValues(i,j) += this->Vinv(k,i) * phisCur(k,j); 00197 } 00198 } 00199 } 00200 } 00201 break; 00202 case OPERATOR_GRAD: 00203 case OPERATOR_D1: 00204 case OPERATOR_D2: 00205 case OPERATOR_D3: 00206 case OPERATOR_D4: 00207 case OPERATOR_D5: 00208 case OPERATOR_D6: 00209 case OPERATOR_D7: 00210 case OPERATOR_D8: 00211 case OPERATOR_D9: 00212 case OPERATOR_D10: 00213 { 00214 const int dkcard = 00215 (operatorType == OPERATOR_GRAD)? getDkCardinality(OPERATOR_D1,2): getDkCardinality(operatorType,2); 00216 00217 FieldContainer<Scalar> phisCur( numBf , numPts , dkcard ); 00218 Phis.getValues( phisCur , inputPoints , operatorType ); 00219 00220 for (int i=0;i<outputValues.dimension(0);i++) { 00221 for (int j=0;j<outputValues.dimension(1);j++) { 00222 for (int k=0;k<outputValues.dimension(2);k++) { 00223 outputValues(i,j,k) = 0.0; 00224 for (int l=0;l<this->getCardinality();l++) { 00225 outputValues(i,j,k) += this->Vinv(l,i) * phisCur(l,j,k); 00226 } 00227 } 00228 } 00229 } 00230 } 00231 break; 00232 case OPERATOR_CURL: // only works in 2d. first component is -d/dy, second is d/dx 00233 { 00234 FieldContainer<Scalar> phisCur( numBf , numPts , getDkCardinality( OPERATOR_D1 , 2 ) ); 00235 Phis.getValues( phisCur , inputPoints , OPERATOR_D1 ); 00236 00237 for (int i=0;i<outputValues.dimension(0);i++) { 00238 for (int j=0;j<outputValues.dimension(1);j++) { 00239 outputValues(i,j,0) = 0.0; 00240 outputValues(i,j,1) = 0.0; 00241 for (int k=0;k<this->getCardinality();k++) { 00242 outputValues(i,j,0) = this->Vinv(k,i) * phisCur(k,j,1); 00243 } 00244 for (int k=0;k<this->getCardinality();k++) { 00245 outputValues(i,j,1) -= this->Vinv(k,i) * phisCur(k,j,0); 00246 } 00247 } 00248 } 00249 } 00250 break; 00251 default: 00252 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00253 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented"); 00254 break; 00255 } 00256 } 00257 catch (std::invalid_argument &exception){ 00258 TEST_FOR_EXCEPTION( true , std::invalid_argument, 00259 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): Operator type not implemented"); 00260 } 00261 00262 } 00263 00264 00265 00266 template<class Scalar, class ArrayScalar> 00267 void Basis_HGRAD_TRI_Cn_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00268 const ArrayScalar & inputPoints, 00269 const ArrayScalar & cellVertices, 00270 const EOperator operatorType) const { 00271 TEST_FOR_EXCEPTION( (true), std::logic_error, 00272 ">>> ERROR (Basis_HGRAD_TRI_Cn_FEM): FEM Basis calling an FVD member function"); 00273 } 00274 00275 00276 }// namespace Intrepid
1.7.4