|
Intrepid
|
00001 #ifndef INTREPID_HGRAD_PYR_C1_FEMDEF_HPP 00002 #define INTREPID_HGRAD_PYR_C1_FEMDEF_HPP 00003 00004 #include <limits> 00005 00006 // @HEADER 00007 // ************************************************************************ 00008 // 00009 // Intrepid Package 00010 // Copyright (2007) Sandia Corporation 00011 // 00012 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00013 // license for use of this work by or on behalf of the U.S. Government. 00014 // 00015 // Redistribution and use in source and binary forms, with or without 00016 // modification, are permitted provided that the following conditions are 00017 // met: 00018 // 00019 // 1. Redistributions of source code must retain the above copyright 00020 // notice, this list of conditions and the following disclaimer. 00021 // 00022 // 2. Redistributions in binary form must reproduce the above copyright 00023 // notice, this list of conditions and the following disclaimer in the 00024 // documentation and/or other materials provided with the distribution. 00025 // 00026 // 3. Neither the name of the Corporation nor the names of the 00027 // contributors may be used to endorse or promote products derived from 00028 // this software without specific prior written permission. 00029 // 00030 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00031 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00032 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00033 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00034 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00035 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00036 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00037 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00038 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00039 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00040 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00041 // 00042 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00043 // Denis Ridzal (dridzal@sandia.gov), or 00044 // Kara Peterson (kjpeter@sandia.gov) 00045 // 00046 // ************************************************************************ 00047 // @HEADER 00048 00054 namespace Intrepid { 00055 00056 template<class Scalar, class ArrayScalar> 00057 Basis_HGRAD_PYR_C1_FEM<Scalar, ArrayScalar>::Basis_HGRAD_PYR_C1_FEM() 00058 { 00059 this -> basisCardinality_ = 5; 00060 this -> basisDegree_ = 1; 00061 this -> basisCellTopology_ = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >() ); 00062 this -> basisType_ = BASIS_FEM_DEFAULT; 00063 this -> basisCoordinates_ = COORDINATES_CARTESIAN; 00064 this -> basisTagsAreSet_ = false; 00065 } 00066 00067 00068 template<class Scalar, class ArrayScalar> 00069 void Basis_HGRAD_PYR_C1_FEM<Scalar, ArrayScalar>::initializeTags() { 00070 00071 // Basis-dependent intializations 00072 int tagSize = 4; // size of DoF tag 00073 int posScDim = 0; // position in the tag, counting from 0, of the subcell dim 00074 int posScOrd = 1; // position in the tag, counting from 0, of the subcell ordinal 00075 int posDfOrd = 2; // position in the tag, counting from 0, of DoF ordinal relative to the subcell 00076 00077 // An array with local DoF tags assigned to basis functions, in the order of their local enumeration 00078 int tags[] = { 0, 0, 0, 1, 00079 0, 1, 0, 1, 00080 0, 2, 0, 1, 00081 0, 3, 0, 1, 00082 0, 4, 0, 1 }; 00083 00084 // Basis-independent function sets tag and enum data in tagToOrdinal_ and ordinalToTag_ arrays: 00085 Intrepid::setOrdinalTagData(this -> tagToOrdinal_, 00086 this -> ordinalToTag_, 00087 tags, 00088 this -> basisCardinality_, 00089 tagSize, 00090 posScDim, 00091 posScOrd, 00092 posDfOrd); 00093 } 00094 00095 00096 00097 template<class Scalar, class ArrayScalar> 00098 void Basis_HGRAD_PYR_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar & outputValues, 00099 const ArrayScalar & inputPoints, 00100 const EOperator operatorType) const { 00101 00102 // Verify arguments 00103 #ifdef HAVE_INTREPID_DEBUG 00104 Intrepid::getValues_HGRAD_Args<Scalar, ArrayScalar>(outputValues, 00105 inputPoints, 00106 operatorType, 00107 this -> getBaseCellTopology(), 00108 this -> getCardinality() ); 00109 #endif 00110 00111 // Number of evaluation points = dim 0 of inputPoints 00112 int dim0 = inputPoints.dimension(0); 00113 00114 // Temporaries: (x,y,z) coordinates of the evaluation point 00115 Scalar x = 0.0; 00116 Scalar y = 0.0; 00117 Scalar z = 0.0; 00118 const Scalar eps = std::numeric_limits<Scalar>::epsilon( ); 00119 00120 switch (operatorType) { 00121 00122 case OPERATOR_VALUE: 00123 for (int i0 = 0; i0 < dim0; i0++) { 00124 x = inputPoints(i0, 0); 00125 y = inputPoints(i0, 1); 00126 z = inputPoints(i0, 2); 00127 00128 //be sure that the basis functions are defined when z is very close to 1. 00129 if(fabs(z-1.0) < eps) { 00130 if(z <= 1.0) z = 1.0-eps; 00131 else z = 1.0+eps; 00132 } 00133 00134 00135 Scalar zTerm = 0.25/(1.0 - z); 00136 00137 // outputValues is a rank-2 array with dimensions (basisCardinality_, dim0) 00138 outputValues(0, i0) = (1.0 - x - z) * (1.0 - y - z) * zTerm; 00139 outputValues(1, i0) = (1.0 + x - z) * (1.0 - y - z) * zTerm; 00140 outputValues(2, i0) = (1.0 + x - z) * (1.0 + y - z) * zTerm; 00141 outputValues(3, i0) = (1.0 - x - z) * (1.0 + y - z) * zTerm; 00142 outputValues(4, i0) = z; 00143 } 00144 break; 00145 00146 case OPERATOR_GRAD: 00147 case OPERATOR_D1: 00148 for (int i0 = 0; i0 < dim0; i0++) { 00149 00150 x = inputPoints(i0, 0); 00151 y = inputPoints(i0, 1); 00152 z = inputPoints(i0, 2); 00153 00154 00155 //be sure that the basis functions are defined when z is very close to 1. 00156 //warning, the derivatives are discontinuous in (0, 0, 1) 00157 if(fabs(z-1.0) < eps) { 00158 if(z <= 1.0) z = 1.0-eps; 00159 else z = 1.0+eps; 00160 } 00161 00162 00163 Scalar zTerm = 0.25/(1.0 - z); 00164 Scalar zTerm2 = 4.0 * zTerm * zTerm; 00165 00166 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, spaceDim) 00167 outputValues(0, i0, 0) = (y + z - 1.0) * zTerm; 00168 outputValues(0, i0, 1) = (x + z - 1.0) * zTerm; 00169 outputValues(0, i0, 2) = x * y * zTerm2 - 0.25; 00170 00171 outputValues(1, i0, 0) = (1.0 - y - z) * zTerm; 00172 outputValues(1, i0, 1) = (z - x - 1.0) * zTerm; 00173 outputValues(1, i0, 2) = - x*y * zTerm2 - 0.25; 00174 00175 outputValues(2, i0, 0) = (1.0 + y - z) * zTerm; 00176 outputValues(2, i0, 1) = (1.0 + x - z) * zTerm; 00177 outputValues(2, i0, 2) = x * y * zTerm2 - 0.25; 00178 00179 outputValues(3, i0, 0) = (z - y - 1.0) * zTerm; 00180 outputValues(3, i0, 1) = (1.0 - x - z) * zTerm; 00181 outputValues(3, i0, 2) = - x*y * zTerm2 - 0.25; 00182 00183 outputValues(4, i0, 0) = 0.0; 00184 outputValues(4, i0, 1) = 0.0; 00185 outputValues(4, i0, 2) = 1; 00186 } 00187 break; 00188 00189 case OPERATOR_CURL: 00190 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_CURL), std::invalid_argument, 00191 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): CURL is invalid operator for rank-0 (scalar) functions in 3D"); 00192 break; 00193 00194 case OPERATOR_DIV: 00195 TEUCHOS_TEST_FOR_EXCEPTION( (operatorType == OPERATOR_DIV), std::invalid_argument, 00196 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): DIV is invalid operator for rank-0 (scalar) functions in 3D"); 00197 break; 00198 00199 case OPERATOR_D2: 00200 for (int i0 = 0; i0 < dim0; i0++) { 00201 x = inputPoints(i0,0); 00202 y = inputPoints(i0,1); 00203 z = inputPoints(i0,2); 00204 00205 //be sure that the basis functions are defined when z is very close to 1. 00206 //warning, the derivatives are discontinuous in (0, 0, 1) 00207 if(fabs(z-1.0) < eps) { 00208 if(z <= 1.0) z = 1.0-eps; 00209 else z = 1.0+eps; 00210 } 00211 00212 00213 Scalar zTerm = 0.25/(1.0 - z); 00214 Scalar zTerm2 = 4.0 * zTerm * zTerm; 00215 Scalar zTerm3 = 8.0 * zTerm * zTerm2; 00216 00217 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, D2Cardinality = 6) 00218 outputValues(0, i0, 0) = 0.0; // {2, 0, 0} 00219 outputValues(0, i0, 1) = zTerm; // {1, 1, 0} 00220 outputValues(0, i0, 2) = y*zTerm2; // {1, 0, 1} 00221 outputValues(0, i0, 3) = 0.0; // {0, 2, 0} 00222 outputValues(0, i0, 4) = x*zTerm2; // {0, 1, 1} 00223 outputValues(0, i0, 5) = x*y*zTerm3; // {0, 0, 2} 00224 00225 outputValues(1, i0, 0) = 0.0; // {2, 0, 0} 00226 outputValues(1, i0, 1) = -zTerm; // {1, 1, 0} 00227 outputValues(1, i0, 2) = -y*zTerm2; // {1, 0, 1} 00228 outputValues(1, i0, 3) = 0.0; // {0, 2, 0} 00229 outputValues(1, i0, 4) = -x*zTerm2; // {0, 1, 1} 00230 outputValues(1, i0, 5) = -x*y*zTerm3; // {0, 0, 2} 00231 00232 outputValues(2, i0, 0) = 0.0; // {2, 0, 0} 00233 outputValues(2, i0, 1) = zTerm; // {1, 1, 0} 00234 outputValues(2, i0, 2) = y*zTerm2; // {1, 0, 1} 00235 outputValues(2, i0, 3) = 0.0; // {0, 2, 0} 00236 outputValues(2, i0, 4) = x*zTerm2; // {0, 1, 1} 00237 outputValues(2, i0, 5) = x*y*zTerm3; // {0, 0, 2} 00238 00239 outputValues(3, i0, 0) = 0.0; // {2, 0, 0} 00240 outputValues(3, i0, 1) = -zTerm; // {1, 1, 0} 00241 outputValues(3, i0, 2) = -y*zTerm2; // {1, 0, 1} 00242 outputValues(3, i0, 3) = 0.0; // {0, 2, 0} 00243 outputValues(3, i0, 4) = -x*zTerm2; // {0, 1, 1} 00244 outputValues(3, i0, 5) = -x*y*zTerm3; // {0, 0, 2} 00245 00246 outputValues(4, i0, 0) = 0.0; // {2, 0, 0} 00247 outputValues(4, i0, 1) = 0.0; // {1, 1, 0} 00248 outputValues(4, i0, 2) = 0.0; // {1, 0, 1} 00249 outputValues(4, i0, 3) = 0.0; // {0, 2, 0} 00250 outputValues(4, i0, 4) = 0.0; // {0, 1, 1} 00251 outputValues(4, i0, 5) = 0.0; // {0, 0, 2} 00252 } 00253 break; 00254 00255 case OPERATOR_D3: 00256 case OPERATOR_D4: 00257 case OPERATOR_D5: 00258 case OPERATOR_D6: 00259 case OPERATOR_D7: 00260 case OPERATOR_D8: 00261 case OPERATOR_D9: 00262 case OPERATOR_D10: 00263 { 00264 // outputValues is a rank-3 array with dimensions (basisCardinality_, dim0, DkCardinality) 00265 int DkCardinality = Intrepid::getDkCardinality(operatorType, 00266 this -> basisCellTopology_.getDimension() ); 00267 for(int dofOrd = 0; dofOrd < this -> basisCardinality_; dofOrd++) { 00268 for (int i0 = 0; i0 < dim0; i0++) { 00269 for(int dkOrd = 0; dkOrd < DkCardinality; dkOrd++){ 00270 outputValues(dofOrd, i0, dkOrd) = 0.0; 00271 } 00272 } 00273 } 00274 } 00275 break; 00276 default: 00277 TEUCHOS_TEST_FOR_EXCEPTION( !( Intrepid::isValidOperator(operatorType) ), std::invalid_argument, 00278 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): Invalid operator type"); 00279 } 00280 } 00281 00282 00283 00284 template<class Scalar, class ArrayScalar> 00285 void Basis_HGRAD_PYR_C1_FEM<Scalar, ArrayScalar>::getValues(ArrayScalar& outputValues, 00286 const ArrayScalar & inputPoints, 00287 const ArrayScalar & cellVertices, 00288 const EOperator operatorType) const { 00289 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error, 00290 ">>> ERROR (Basis_HGRAD_PYR_C1_FEM): FEM Basis calling an FVD member function"); 00291 } 00292 }// namespace Intrepid 00293 #endif
1.7.6.1