quadrature_monomial.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 #include "libmesh/quadrature_monomial.h" 00020 00021 namespace libMesh 00022 { 00023 00024 00025 // See the files: 00026 00027 // quadrature_monomial_2D.C 00028 // quadrature_monomial_3D.C 00029 00030 // for implementation of specific element types. 00031 00032 // Constructor 00033 QMonomial::QMonomial(const unsigned int d, 00034 const Order o) : QBase(d,o) 00035 { 00036 } 00037 00038 00039 // Destructor 00040 QMonomial::~QMonomial() 00041 { 00042 } 00043 00044 void QMonomial::wissmann_rule(const Real rule_data[][3], 00045 const unsigned int n_pts) 00046 { 00047 for (unsigned int i=0, c=0; i<n_pts; ++i) 00048 { 00049 _points[c] = Point( rule_data[i][0], rule_data[i][1] ); 00050 _weights[c++] = rule_data[i][2]; 00051 00052 // This may be an (x1,x2) -> (-x1,x2) point, in which case 00053 // we will also generate the mirror point using the same weight. 00054 if (rule_data[i][0] != static_cast<Real>(0.0)) 00055 { 00056 _points[c] = Point( -rule_data[i][0], rule_data[i][1] ); 00057 _weights[c++] = rule_data[i][2]; 00058 } 00059 } 00060 } 00061 00062 00063 00064 void QMonomial::stroud_rule(const Real rule_data[][3], 00065 const unsigned int* rule_symmetry, 00066 const unsigned int n_pts) 00067 { 00068 for (unsigned int i=0, c=0; i<n_pts; ++i) 00069 { 00070 const Real 00071 x=rule_data[i][0], 00072 y=rule_data[i][1], 00073 wt=rule_data[i][2]; 00074 00075 switch(rule_symmetry[i]) 00076 { 00077 case 0: // Single point (no symmetry) 00078 { 00079 _points[c] = Point( x, y); 00080 _weights[c++] = wt; 00081 00082 break; 00083 } 00084 case 1: // Fully-symmetric (x,y) 00085 { 00086 _points[c] = Point( x, y); 00087 _weights[c++] = wt; 00088 00089 _points[c] = Point(-x, y); 00090 _weights[c++] = wt; 00091 00092 _points[c] = Point( x,-y); 00093 _weights[c++] = wt; 00094 00095 _points[c] = Point(-x,-y); 00096 _weights[c++] = wt; 00097 00098 _points[c] = Point( y, x); 00099 _weights[c++] = wt; 00100 00101 _points[c] = Point(-y, x); 00102 _weights[c++] = wt; 00103 00104 _points[c] = Point( y,-x); 00105 _weights[c++] = wt; 00106 00107 _points[c] = Point(-y,-x); 00108 _weights[c++] = wt; 00109 00110 break; 00111 } 00112 case 2: // Fully-symmetric (x,x) 00113 { 00114 _points[c] = Point( x, x); 00115 _weights[c++] = wt; 00116 00117 _points[c] = Point(-x, x); 00118 _weights[c++] = wt; 00119 00120 _points[c] = Point( x,-x); 00121 _weights[c++] = wt; 00122 00123 _points[c] = Point(-x,-x); 00124 _weights[c++] = wt; 00125 00126 break; 00127 } 00128 case 3: // Fully-symmetric (x,0) 00129 { 00130 libmesh_assert_equal_to (y, 0.0); 00131 00132 _points[c] = Point( x,0.); 00133 _weights[c++] = wt; 00134 00135 _points[c] = Point(-x,0.); 00136 _weights[c++] = wt; 00137 00138 _points[c] = Point(0., x); 00139 _weights[c++] = wt; 00140 00141 _points[c] = Point(0.,-x); 00142 _weights[c++] = wt; 00143 00144 break; 00145 } 00146 case 4: // Rotational invariant 00147 { 00148 _points[c] = Point( x, y); 00149 _weights[c++] = wt; 00150 00151 _points[c] = Point(-x,-y); 00152 _weights[c++] = wt; 00153 00154 _points[c] = Point(-y, x); 00155 _weights[c++] = wt; 00156 00157 _points[c] = Point( y,-x); 00158 _weights[c++] = wt; 00159 00160 break; 00161 } 00162 case 5: // Partial symmetry (Wissman's rules) 00163 { 00164 libmesh_assert_not_equal_to (x, 0.0); 00165 00166 _points[c] = Point( x, y); 00167 _weights[c++] = wt; 00168 00169 _points[c] = Point(-x, y); 00170 _weights[c++] = wt; 00171 00172 break; 00173 } 00174 case 6: // Rectangular symmetry 00175 { 00176 _points[c] = Point( x, y); 00177 _weights[c++] = wt; 00178 00179 _points[c] = Point(-x, y); 00180 _weights[c++] = wt; 00181 00182 _points[c] = Point(-x,-y); 00183 _weights[c++] = wt; 00184 00185 _points[c] = Point( x,-y); 00186 _weights[c++] = wt; 00187 00188 break; 00189 } 00190 case 7: // Central symmetry 00191 { 00192 libmesh_assert_equal_to (x, 0.0); 00193 libmesh_assert_not_equal_to (y, 0.0); 00194 00195 _points[c] = Point(0., y); 00196 _weights[c++] = wt; 00197 00198 _points[c] = Point(0.,-y); 00199 _weights[c++] = wt; 00200 00201 break; 00202 } 00203 default: 00204 libmesh_error_msg("Unknown symmetry!"); 00205 } // end switch(rule_symmetry[i]) 00206 } 00207 } 00208 00209 00210 00211 00212 void QMonomial::kim_rule(const Real rule_data[][4], 00213 const unsigned int* rule_id, 00214 const unsigned int n_pts) 00215 { 00216 for (unsigned int i=0, c=0; i<n_pts; ++i) 00217 { 00218 const Real 00219 x=rule_data[i][0], 00220 y=rule_data[i][1], 00221 z=rule_data[i][2], 00222 wt=rule_data[i][3]; 00223 00224 switch(rule_id[i]) 00225 { 00226 case 0: // (0,0,0) 1 permutation 00227 { 00228 _points[c] = Point( x, y, z); _weights[c++] = wt; 00229 00230 break; 00231 } 00232 case 1: // (x,0,0) 6 permutations 00233 { 00234 _points[c] = Point( x, 0., 0.); _weights[c++] = wt; 00235 _points[c] = Point(0., -x, 0.); _weights[c++] = wt; 00236 _points[c] = Point(-x, 0., 0.); _weights[c++] = wt; 00237 _points[c] = Point(0., x, 0.); _weights[c++] = wt; 00238 _points[c] = Point(0., 0., -x); _weights[c++] = wt; 00239 _points[c] = Point(0., 0., x); _weights[c++] = wt; 00240 00241 break; 00242 } 00243 case 2: // (x,x,0) 12 permutations 00244 { 00245 _points[c] = Point( x, x, 0.); _weights[c++] = wt; 00246 _points[c] = Point( x, -x, 0.); _weights[c++] = wt; 00247 _points[c] = Point(-x, -x, 0.); _weights[c++] = wt; 00248 _points[c] = Point(-x, x, 0.); _weights[c++] = wt; 00249 _points[c] = Point( x, 0., -x); _weights[c++] = wt; 00250 _points[c] = Point( x, 0., x); _weights[c++] = wt; 00251 _points[c] = Point(0., x, -x); _weights[c++] = wt; 00252 _points[c] = Point(0., x, x); _weights[c++] = wt; 00253 _points[c] = Point(0., -x, -x); _weights[c++] = wt; 00254 _points[c] = Point(-x, 0., -x); _weights[c++] = wt; 00255 _points[c] = Point(0., -x, x); _weights[c++] = wt; 00256 _points[c] = Point(-x, 0., x); _weights[c++] = wt; 00257 00258 break; 00259 } 00260 case 3: // (x,y,0) 24 permutations 00261 { 00262 _points[c] = Point( x, y, 0.); _weights[c++] = wt; 00263 _points[c] = Point( y, -x, 0.); _weights[c++] = wt; 00264 _points[c] = Point(-x, -y, 0.); _weights[c++] = wt; 00265 _points[c] = Point(-y, x, 0.); _weights[c++] = wt; 00266 _points[c] = Point( x, 0., -y); _weights[c++] = wt; 00267 _points[c] = Point( x, -y, 0.); _weights[c++] = wt; 00268 _points[c] = Point( x, 0., y); _weights[c++] = wt; 00269 _points[c] = Point(0., y, -x); _weights[c++] = wt; 00270 _points[c] = Point(-x, y, 0.); _weights[c++] = wt; 00271 _points[c] = Point(0., y, x); _weights[c++] = wt; 00272 _points[c] = Point( y, 0., -x); _weights[c++] = wt; 00273 _points[c] = Point(0., -y, -x); _weights[c++] = wt; 00274 _points[c] = Point(-y, 0., -x); _weights[c++] = wt; 00275 _points[c] = Point( y, x, 0.); _weights[c++] = wt; 00276 _points[c] = Point(-y, -x, 0.); _weights[c++] = wt; 00277 _points[c] = Point( y, 0., x); _weights[c++] = wt; 00278 _points[c] = Point(0., -y, x); _weights[c++] = wt; 00279 _points[c] = Point(-y, 0., x); _weights[c++] = wt; 00280 _points[c] = Point(-x, 0., y); _weights[c++] = wt; 00281 _points[c] = Point(0., -x, -y); _weights[c++] = wt; 00282 _points[c] = Point(0., -x, y); _weights[c++] = wt; 00283 _points[c] = Point(-x, 0., -y); _weights[c++] = wt; 00284 _points[c] = Point(0., x, y); _weights[c++] = wt; 00285 _points[c] = Point(0., x, -y); _weights[c++] = wt; 00286 00287 break; 00288 } 00289 case 4: // (x,x,x) 8 permutations 00290 { 00291 _points[c] = Point( x, x, x); _weights[c++] = wt; 00292 _points[c] = Point( x, -x, x); _weights[c++] = wt; 00293 _points[c] = Point(-x, -x, x); _weights[c++] = wt; 00294 _points[c] = Point(-x, x, x); _weights[c++] = wt; 00295 _points[c] = Point( x, x, -x); _weights[c++] = wt; 00296 _points[c] = Point( x, -x, -x); _weights[c++] = wt; 00297 _points[c] = Point(-x, x, -x); _weights[c++] = wt; 00298 _points[c] = Point(-x, -x, -x); _weights[c++] = wt; 00299 00300 break; 00301 } 00302 case 5: // (x,x,z) 24 permutations 00303 { 00304 _points[c] = Point( x, x, z); _weights[c++] = wt; 00305 _points[c] = Point( x, -x, z); _weights[c++] = wt; 00306 _points[c] = Point(-x, -x, z); _weights[c++] = wt; 00307 _points[c] = Point(-x, x, z); _weights[c++] = wt; 00308 _points[c] = Point( x, z, -x); _weights[c++] = wt; 00309 _points[c] = Point( x, -x, -z); _weights[c++] = wt; 00310 _points[c] = Point( x, -z, x); _weights[c++] = wt; 00311 _points[c] = Point( z, x, -x); _weights[c++] = wt; 00312 _points[c] = Point(-x, x, -z); _weights[c++] = wt; 00313 _points[c] = Point(-z, x, x); _weights[c++] = wt; 00314 _points[c] = Point( x, -z, -x); _weights[c++] = wt; 00315 _points[c] = Point(-z, -x, -x); _weights[c++] = wt; 00316 _points[c] = Point(-x, z, -x); _weights[c++] = wt; 00317 _points[c] = Point( x, x, -z); _weights[c++] = wt; 00318 _points[c] = Point(-x, -x, -z); _weights[c++] = wt; 00319 _points[c] = Point( x, z, x); _weights[c++] = wt; 00320 _points[c] = Point( z, -x, x); _weights[c++] = wt; 00321 _points[c] = Point(-x, -z, x); _weights[c++] = wt; 00322 _points[c] = Point(-x, z, x); _weights[c++] = wt; 00323 _points[c] = Point( z, -x, -x); _weights[c++] = wt; 00324 _points[c] = Point(-z, -x, x); _weights[c++] = wt; 00325 _points[c] = Point(-x, -z, -x); _weights[c++] = wt; 00326 _points[c] = Point( z, x, x); _weights[c++] = wt; 00327 _points[c] = Point(-z, x, -x); _weights[c++] = wt; 00328 00329 break; 00330 } 00331 case 6: // (x,y,z) 24 permutations 00332 { 00333 _points[c] = Point( x, y, z); _weights[c++] = wt; 00334 _points[c] = Point( y, -x, z); _weights[c++] = wt; 00335 _points[c] = Point(-x, -y, z); _weights[c++] = wt; 00336 _points[c] = Point(-y, x, z); _weights[c++] = wt; 00337 _points[c] = Point( x, z, -y); _weights[c++] = wt; 00338 _points[c] = Point( x, -y, -z); _weights[c++] = wt; 00339 _points[c] = Point( x, -z, y); _weights[c++] = wt; 00340 _points[c] = Point( z, y, -x); _weights[c++] = wt; 00341 _points[c] = Point(-x, y, -z); _weights[c++] = wt; 00342 _points[c] = Point(-z, y, x); _weights[c++] = wt; 00343 _points[c] = Point( y, -z, -x); _weights[c++] = wt; 00344 _points[c] = Point(-z, -y, -x); _weights[c++] = wt; 00345 _points[c] = Point(-y, z, -x); _weights[c++] = wt; 00346 _points[c] = Point( y, x, -z); _weights[c++] = wt; 00347 _points[c] = Point(-y, -x, -z); _weights[c++] = wt; 00348 _points[c] = Point( y, z, x); _weights[c++] = wt; 00349 _points[c] = Point( z, -y, x); _weights[c++] = wt; 00350 _points[c] = Point(-y, -z, x); _weights[c++] = wt; 00351 _points[c] = Point(-x, z, y); _weights[c++] = wt; 00352 _points[c] = Point( z, -x, -y); _weights[c++] = wt; 00353 _points[c] = Point(-z, -x, y); _weights[c++] = wt; 00354 _points[c] = Point(-x, -z, -y); _weights[c++] = wt; 00355 _points[c] = Point( z, x, y); _weights[c++] = wt; 00356 _points[c] = Point(-z, x, -y); _weights[c++] = wt; 00357 00358 break; 00359 } 00360 default: 00361 libmesh_error_msg("Unknown rule ID: " << rule_id[i] << "!"); 00362 } // end switch(rule_id[i]) 00363 } 00364 } 00365 00366 } // namespace libMesh
generated by
