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