|
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 00038 template<class Scalar, class ArrayScalar> 00039 Basis_HGRAD_TRI_C2_FEM<Scalar,ArrayScalar>::Basis_HGRAD_TRI_C2_FEM() 00040 { 00041 this -> basisCardinality_ = 6; 00042 this -> basisDegree_ = 2; 00043 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >() ); 00044 this -> basisType_ = BASIS_FEM_DEFAULT; 00045 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00046 this -> basisTagsAreSet_ = false; 00047 } 00048 00049 00050 00051 template<class Scalar, class ArrayScalar> 00052 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::initializeTags() { 00053 00054 // Basis-dependent initializations 00055 int tagSize = 4; // size of DoF tag, i.e., number of fields in the tag 00056 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00057 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00058 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00059 00060 // An array with local DoF tags assigned to the basis functions, in the order of their local enumeration 00061 int tags[] = { 0, 0, 0, 1, 00062 0, 1, 0, 1, 00063 0, 2, 0, 1, 00064 1, 0, 0, 1, 00065 1, 1, 0, 1, 00066 1, 2, 0, 1}; 00067 00068 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00069 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00070 this -> ordinalToTag_, 00071 tags, 00072 this -> basisCardinality_, 00073 tagSize, 00074 posScDim, 00075 posScOrd, 00076 posDfOrd); 00077 } 00078 00079 00080 00081 template<class Scalar, class ArrayScalar> 00082 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00083 const ArrayScalar & inputPoints, 00084 const EOperator operatorType) const { 00085 00086 // Verify arguments 00087 #ifdef HAVE_INTREPID_DEBUG 00088 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00089 inputPoints, 00090 operatorType, 00091 this -> getBaseCellTopology(), 00092 this -> getCardinality() ); 00093 #endif 00094 00095 // Number of evaluation points = dim 0 of inputPoints 00096 int dim0 = inputPoints.dimension(0); 00097 00098 // Temporaries: (x,y) coordinates of the evaluation point 00099 Scalar x = 0.0; 00100 Scalar y = 0.0; 00101 00102 switch (operatorType) { 00103 00104 case OPERATOR_VALUE: 00105 for (int i0 = 0; i0 < dim0; i0++) { 00106 x = inputPoints(i0, 0); 00107 y = inputPoints(i0, 1); 00108 00109 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00110 outputValues(0, i0) = (x + y - 1.0)*(2.0*x + 2.0*y - 1.0); 00111 outputValues(1, i0) = x*(2.0*x - 1.0); 00112 outputValues(2, i0) = y*(2.0*y - 1.0); 00113 outputValues(3, i0) = -4.0*x*(x + y - 1.0); 00114 outputValues(4, i0) = 4.0*x*y; 00115 outputValues(5, i0) = -4.0*y*(x + y - 1.0); 00116 00117 } 00118 break; 00119 00120 case OPERATOR_GRAD: 00121 case OPERATOR_D1: 00122 for (int i0 = 0; i0 < dim0; i0++) { 00123 x = inputPoints(i0, 0); 00124 y = inputPoints(i0, 1); 00125 00126 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00127 outputValues(0, i0, 0) = 4.0*x + 4.0*y - 3.0; 00128 outputValues(0, i0, 1) = 4.0*x + 4.0*y - 3.0; 00129 00130 outputValues(1, i0, 0) = 4.0*x - 1.0; 00131 outputValues(1, i0, 1) = 0.0; 00132 00133 outputValues(2, i0, 0) = 0.0; 00134 outputValues(2, i0, 1) = 4.0*y - 1.0; 00135 00136 outputValues(3, i0, 0) = -4.0*(2.0*x + y - 1.0); 00137 outputValues(3, i0, 1) = -4.0*x; 00138 00139 outputValues(4, i0, 0) = 4.0*y; 00140 outputValues(4, i0, 1) = 4.0*x; 00141 00142 outputValues(5, i0, 0) = -4.0*y; 00143 outputValues(5, i0, 1) = -4.0*(x + 2.0*y - 1.0); 00144 } 00145 break; 00146 00147 case OPERATOR_CURL: 00148 for (int i0 = 0; i0 < dim0; i0++) { 00149 x = inputPoints(i0, 0); 00150 y = inputPoints(i0, 1); 00151 00152 // CURL(u) = (u_y, -u_x), is rotated GRAD 00153 outputValues(0, i0, 1) =-(4.0*x + 4.0*y - 3.0); 00154 outputValues(0, i0, 0) = 4.0*x + 4.0*y - 3.0; 00155 00156 outputValues(1, i0, 1) =-(4.0*x - 1.0); 00157 outputValues(1, i0, 0) = 0.0; 00158 00159 outputValues(2, i0, 1) = 0.0; 00160 outputValues(2, i0, 0) = 4.0*y - 1.0; 00161 00162 outputValues(3, i0, 1) = 4.0*(2.0*x + y - 1.0); 00163 outputValues(3, i0, 0) = -4.0*x; 00164 00165 outputValues(4, i0, 1) = -4.0*y; 00166 outputValues(4, i0, 0) = 4.0*x; 00167 00168 outputValues(5, i0, 1) = 4.0*y; 00169 outputValues(5, i0, 0) = -4.0*(x + 2.0*y - 1.0); 00170 } 00171 break; 00172 00173 case OPERATOR_DIV: 00174 TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00175 ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): DIV is invalid operator for rank-0 (scalar) fields in 2D."); 00176 break; 00177 00178 case OPERATOR_D2: 00179 for (int i0 = 0; i0 < dim0; i0++) { 00180 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00181 // D2 -> (2,0) -> dx^2. 00182 outputValues(0, i0, 0) = 4.0; 00183 outputValues(1, i0, 0) = 4.0; 00184 outputValues(2, i0, 0) = 0.0; 00185 outputValues(3, i0, 0) =-8.0; 00186 outputValues(4, i0, 0) = 0.0; 00187 outputValues(5, i0, 0) = 0.0; 00188 00189 // D2 -> (1,1) -> dx dy 00190 outputValues(0, i0, 1) = 4.0; 00191 outputValues(1, i0, 1) = 0.0; 00192 outputValues(2, i0, 1) = 0.0; 00193 outputValues(3, i0, 1) =-4.0; 00194 outputValues(4, i0, 1) = 4.0; 00195 outputValues(5, i0, 1) =-4.0; 00196 00197 // D2 -> (0,2) -> dy^2 00198 outputValues(0, i0, 2) = 4.0; 00199 outputValues(1, i0, 2) = 0.0; 00200 outputValues(2, i0, 2) = 4.0; 00201 outputValues(3, i0, 2) = 0.0; 00202 outputValues(4, i0, 2) = 0.0; 00203 outputValues(5, i0, 2) =-8.0; 00204 }// for i0 00205 break; 00206 00207 case OPERATOR_D3: 00208 case OPERATOR_D4: 00209 case OPERATOR_D5: 00210 case OPERATOR_D6: 00211 case OPERATOR_D7: 00212 case OPERATOR_D8: 00213 case OPERATOR_D9: 00214 case OPERATOR_D10: 00215 { 00216 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00217 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00218 this -> basisCellTopology_.getDimension() ); 00219 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00220 for (int i0 = 0; i0 < dim0; i0++) { 00221 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00222 outputValues(dofOrd, i0, dkOrd) = 0.0; 00223 } 00224 } 00225 } 00226 } 00227 break; 00228 00229 default: 00230 TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00231 ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): Invalid operator type"); 00232 } 00233 } 00234 00235 00236 00237 template<class Scalar, class ArrayScalar> 00238 void Basis_HGRAD_TRI_C2_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00239 const ArrayScalar & inputPoints, 00240 const ArrayScalar & cellVertices, 00241 const EOperator operatorType) const { 00242 TEST_FOR_EXCEPTION( (true), std::logic_error, 00243 ">>> ERROR (Basis_HGRAD_TRI_C2_FEM): FEM Basis calling an FVD member function"); 00244 } 00245 00246 00247 }// namespace Intrepid
1.7.4