quadrature_gm_3D.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 00020 // Local includes 00021 #include "libmesh/quadrature_gm.h" 00022 00023 namespace libMesh 00024 { 00025 00026 00027 00028 void QGrundmann_Moller::init_3D(const ElemType type_in, 00029 unsigned int p) 00030 { 00031 // Nearly all GM rules contain negative weights, so if you are not 00032 // allowing rules with negative weights, we cannot continue! 00033 if (!allow_rules_with_negative_weights) 00034 libmesh_error_msg("You requested a Grundmann-Moller rule but\n" \ 00035 << "are not allowing rules with negative weights!\n" \ 00036 << "Either select a different quadrature class or\n" \ 00037 << "set allow_rules_with_negative_weights==true."); 00038 00039 switch (type_in) 00040 { 00041 case TET4: 00042 case TET10: 00043 { 00044 switch(_order + 2*p) 00045 { 00046 // We hard-code the first few orders based on output from 00047 // the mp-quadrature library: 00048 // 00049 // https://code.google.com/p/mp-quadrature 00050 // 00051 // The points are ordered in such a way that the amount of 00052 // round-off error in the quadrature calculations is 00053 // (hopefully) minimized. These orderings were obtained 00054 // via a simple permutation optimization strategy designed 00055 // to produce the smallest overall average error while 00056 // integrating all polynomials of degree <= d. 00057 case CONSTANT: 00058 case FIRST: 00059 { 00060 _points.resize(1); 00061 _weights.resize(1); 00062 00063 _points[0](0) = 1.L/4.L; 00064 _points[0](1) = 1.L/4.L; 00065 _points[0](2) = 1.L/4.L; 00066 00067 _weights[0] = 1.L/6.L; 00068 return; 00069 } 00070 00071 case SECOND: 00072 case THIRD: 00073 { 00074 _points.resize(5); 00075 _weights.resize(5); 00076 00077 _points[0](0) = 1.L/6.L; _points[0](1) = 1.L/6.L; _points[0](2) = 1.L/2.L; _weights[0] = 3.L/40.L; 00078 _points[1](0) = 1.L/6.L; _points[1](1) = 1.L/6.L; _points[1](2) = 1.L/6.L; _weights[1] = 3.L/40.L; 00079 _points[2](0) = 1.L/4.L; _points[2](1) = 1.L/4.L; _points[2](2) = 1.L/4.L; _weights[2] = -2.L/15.L; 00080 _points[3](0) = 1.L/2.L; _points[3](1) = 1.L/6.L; _points[3](2) = 1.L/6.L; _weights[3] = 3.L/40.L; 00081 _points[4](0) = 1.L/6.L; _points[4](1) = 1.L/2.L; _points[4](2) = 1.L/6.L; _weights[4] = 3.L/40.L; 00082 return; 00083 } 00084 00085 case FOURTH: 00086 case FIFTH: 00087 { 00088 _points.resize(15); 00089 _weights.resize(15); 00090 00091 _points[ 0](0) = 1.L/8.L; _points[ 0](1) = 3.L/8.L; _points[ 0](2) = 1.L/8.L; _weights[ 0] = 16.L/315.L; 00092 _points[ 1](0) = 1.L/8.L; _points[ 1](1) = 1.L/8.L; _points[ 1](2) = 5.L/8.L; _weights[ 1] = 16.L/315.L; 00093 _points[ 2](0) = 3.L/8.L; _points[ 2](1) = 1.L/8.L; _points[ 2](2) = 1.L/8.L; _weights[ 2] = 16.L/315.L; 00094 _points[ 3](0) = 1.L/6.L; _points[ 3](1) = 1.L/2.L; _points[ 3](2) = 1.L/6.L; _weights[ 3] = -27.L/280.L; 00095 _points[ 4](0) = 3.L/8.L; _points[ 4](1) = 1.L/8.L; _points[ 4](2) = 3.L/8.L; _weights[ 4] = 16.L/315.L; 00096 _points[ 5](0) = 1.L/8.L; _points[ 5](1) = 3.L/8.L; _points[ 5](2) = 3.L/8.L; _weights[ 5] = 16.L/315.L; 00097 _points[ 6](0) = 1.L/6.L; _points[ 6](1) = 1.L/6.L; _points[ 6](2) = 1.L/6.L; _weights[ 6] = -27.L/280.L; 00098 _points[ 7](0) = 1.L/6.L; _points[ 7](1) = 1.L/6.L; _points[ 7](2) = 1.L/2.L; _weights[ 7] = -27.L/280.L; 00099 _points[ 8](0) = 1.L/8.L; _points[ 8](1) = 1.L/8.L; _points[ 8](2) = 1.L/8.L; _weights[ 8] = 16.L/315.L; 00100 _points[ 9](0) = 1.L/4.L; _points[ 9](1) = 1.L/4.L; _points[ 9](2) = 1.L/4.L; _weights[ 9] = 2.L/45.L; 00101 _points[10](0) = 1.L/8.L; _points[10](1) = 5.L/8.L; _points[10](2) = 1.L/8.L; _weights[10] = 16.L/315.L; 00102 _points[11](0) = 1.L/2.L; _points[11](1) = 1.L/6.L; _points[11](2) = 1.L/6.L; _weights[11] = -27.L/280.L; 00103 _points[12](0) = 1.L/8.L; _points[12](1) = 1.L/8.L; _points[12](2) = 3.L/8.L; _weights[12] = 16.L/315.L; 00104 _points[13](0) = 5.L/8.L; _points[13](1) = 1.L/8.L; _points[13](2) = 1.L/8.L; _weights[13] = 16.L/315.L; 00105 _points[14](0) = 3.L/8.L; _points[14](1) = 3.L/8.L; _points[14](2) = 1.L/8.L; _weights[14] = 16.L/315.L; 00106 return; 00107 } 00108 00109 case SIXTH: 00110 case SEVENTH: 00111 { 00112 _points.resize(35); 00113 _weights.resize(35); 00114 00115 _points[ 0](0) = 3.L/10.L; _points[ 0](1) = 1.L/10.L; _points[ 0](2) = 3.L/10.L; _weights[ 0] = 3125.L/72576.L; 00116 _points[ 1](0) = 1.L/6.L; _points[ 1](1) = 1.L/2.L; _points[ 1](2) = 1.L/6.L; _weights[ 1] = 243.L/4480.L; 00117 _points[ 2](0) = 1.L/6.L; _points[ 2](1) = 1.L/6.L; _points[ 2](2) = 1.L/2.L; _weights[ 2] = 243.L/4480.L; 00118 _points[ 3](0) = 1.L/2.L; _points[ 3](1) = 1.L/10.L; _points[ 3](2) = 1.L/10.L; _weights[ 3] = 3125.L/72576.L; 00119 _points[ 4](0) = 3.L/10.L; _points[ 4](1) = 1.L/10.L; _points[ 4](2) = 1.L/10.L; _weights[ 4] = 3125.L/72576.L; 00120 _points[ 5](0) = 3.L/10.L; _points[ 5](1) = 3.L/10.L; _points[ 5](2) = 1.L/10.L; _weights[ 5] = 3125.L/72576.L; 00121 _points[ 6](0) = 1.L/2.L; _points[ 6](1) = 1.L/6.L; _points[ 6](2) = 1.L/6.L; _weights[ 6] = 243.L/4480.L; 00122 _points[ 7](0) = 3.L/10.L; _points[ 7](1) = 1.L/10.L; _points[ 7](2) = 1.L/2.L; _weights[ 7] = 3125.L/72576.L; 00123 _points[ 8](0) = 7.L/10.L; _points[ 8](1) = 1.L/10.L; _points[ 8](2) = 1.L/10.L; _weights[ 8] = 3125.L/72576.L; 00124 _points[ 9](0) = 3.L/8.L; _points[ 9](1) = 1.L/8.L; _points[ 9](2) = 1.L/8.L; _weights[ 9] = -256.L/2835.L; 00125 _points[10](0) = 3.L/10.L; _points[10](1) = 3.L/10.L; _points[10](2) = 3.L/10.L; _weights[10] = 3125.L/72576.L; 00126 _points[11](0) = 1.L/10.L; _points[11](1) = 1.L/2.L; _points[11](2) = 3.L/10.L; _weights[11] = 3125.L/72576.L; 00127 _points[12](0) = 1.L/10.L; _points[12](1) = 1.L/10.L; _points[12](2) = 7.L/10.L; _weights[12] = 3125.L/72576.L; 00128 _points[13](0) = 1.L/2.L; _points[13](1) = 1.L/10.L; _points[13](2) = 3.L/10.L; _weights[13] = 3125.L/72576.L; 00129 _points[14](0) = 1.L/10.L; _points[14](1) = 7.L/10.L; _points[14](2) = 1.L/10.L; _weights[14] = 3125.L/72576.L; 00130 _points[15](0) = 1.L/10.L; _points[15](1) = 1.L/2.L; _points[15](2) = 1.L/10.L; _weights[15] = 3125.L/72576.L; 00131 _points[16](0) = 1.L/6.L; _points[16](1) = 1.L/6.L; _points[16](2) = 1.L/6.L; _weights[16] = 243.L/4480.L; 00132 _points[17](0) = 3.L/8.L; _points[17](1) = 1.L/8.L; _points[17](2) = 3.L/8.L; _weights[17] = -256.L/2835.L; 00133 _points[18](0) = 1.L/8.L; _points[18](1) = 1.L/8.L; _points[18](2) = 5.L/8.L; _weights[18] = -256.L/2835.L; 00134 _points[19](0) = 1.L/10.L; _points[19](1) = 1.L/10.L; _points[19](2) = 3.L/10.L; _weights[19] = 3125.L/72576.L; 00135 _points[20](0) = 1.L/8.L; _points[20](1) = 3.L/8.L; _points[20](2) = 3.L/8.L; _weights[20] = -256.L/2835.L; 00136 _points[21](0) = 5.L/8.L; _points[21](1) = 1.L/8.L; _points[21](2) = 1.L/8.L; _weights[21] = -256.L/2835.L; 00137 _points[22](0) = 1.L/8.L; _points[22](1) = 5.L/8.L; _points[22](2) = 1.L/8.L; _weights[22] = -256.L/2835.L; 00138 _points[23](0) = 1.L/10.L; _points[23](1) = 3.L/10.L; _points[23](2) = 1.L/10.L; _weights[23] = 3125.L/72576.L; 00139 _points[24](0) = 1.L/4.L; _points[24](1) = 1.L/4.L; _points[24](2) = 1.L/4.L; _weights[24] = -8.L/945.L; 00140 _points[25](0) = 1.L/8.L; _points[25](1) = 1.L/8.L; _points[25](2) = 3.L/8.L; _weights[25] = -256.L/2835.L; 00141 _points[26](0) = 3.L/8.L; _points[26](1) = 3.L/8.L; _points[26](2) = 1.L/8.L; _weights[26] = -256.L/2835.L; 00142 _points[27](0) = 1.L/8.L; _points[27](1) = 3.L/8.L; _points[27](2) = 1.L/8.L; _weights[27] = -256.L/2835.L; 00143 _points[28](0) = 1.L/10.L; _points[28](1) = 3.L/10.L; _points[28](2) = 1.L/2.L; _weights[28] = 3125.L/72576.L; 00144 _points[29](0) = 3.L/10.L; _points[29](1) = 1.L/2.L; _points[29](2) = 1.L/10.L; _weights[29] = 3125.L/72576.L; 00145 _points[30](0) = 1.L/10.L; _points[30](1) = 1.L/10.L; _points[30](2) = 1.L/2.L; _weights[30] = 3125.L/72576.L; 00146 _points[31](0) = 1.L/2.L; _points[31](1) = 3.L/10.L; _points[31](2) = 1.L/10.L; _weights[31] = 3125.L/72576.L; 00147 _points[32](0) = 1.L/8.L; _points[32](1) = 1.L/8.L; _points[32](2) = 1.L/8.L; _weights[32] = -256.L/2835.L; 00148 _points[33](0) = 1.L/10.L; _points[33](1) = 3.L/10.L; _points[33](2) = 3.L/10.L; _weights[33] = 3125.L/72576.L; 00149 _points[34](0) = 1.L/10.L; _points[34](1) = 1.L/10.L; _points[34](2) = 1.L/10.L; _weights[34] = 3125.L/72576.L; 00150 return; 00151 } 00152 00153 default: 00154 { 00155 // Untested above _order=23 but should work... 00156 gm_rule((_order + 2*p)/2, /*dim=*/3); 00157 return; 00158 } 00159 } // end switch (order) 00160 } // end case TET4, TET10 00161 00162 default: 00163 libmesh_error_msg("ERROR: Unsupported element type: " << type_in); 00164 } // end switch (type_in) 00165 } 00166 00167 } // namespace libMesh
generated by
