Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : // Local includes
      21             : #include "libmesh/quadrature_gm.h"
      22             : #include "libmesh/enum_to_string.h"
      23             : 
      24             : namespace libMesh
      25             : {
      26             : 
      27        6607 : void QGrundmann_Moller::init_3D()
      28             : {
      29             :   // Nearly all GM rules contain negative weights, so if you are not
      30             :   // allowing rules with negative weights, we cannot continue!
      31        6607 :   libmesh_error_msg_if(!allow_rules_with_negative_weights,
      32             :                        "You requested a Grundmann-Moller rule but\n"
      33             :                        "are not allowing rules with negative weights!\n"
      34             :                        "Either select a different quadrature class or\n"
      35             :                        "set allow_rules_with_negative_weights==true.");
      36             : 
      37        6607 :   switch (_type)
      38             :     {
      39        6607 :     case TET4:
      40             :     case TET10:
      41             :     case TET14:
      42             :       {
      43        6881 :         switch(get_order())
      44             :           {
      45             :             // We hard-code the first few orders based on output from
      46             :             // the mp-quadrature library:
      47             :             //
      48             :             // https://code.google.com/p/mp-quadrature
      49             :             //
      50             :             // The points are ordered in such a way that the amount of
      51             :             // round-off error in the quadrature calculations is
      52             :             // (hopefully) minimized.  These orderings were obtained
      53             :             // via a simple permutation optimization strategy designed
      54             :             // to produce the smallest overall average error while
      55             :             // integrating all polynomials of degree <= d.
      56         142 :           case CONSTANT:
      57             :           case FIRST:
      58             :             {
      59         142 :               _points.resize(1);
      60         142 :               _weights.resize(1);
      61             : 
      62         142 :               _points[0](0) = Real(1)/4;
      63         142 :               _points[0](1) = Real(1)/4;
      64         142 :               _points[0](2) = Real(1)/4;
      65             : 
      66         142 :               _weights[0] = Real(1)/6;
      67         142 :               return;
      68             :             }
      69             : 
      70         142 :           case SECOND:
      71             :           case THIRD:
      72             :             {
      73         142 :               _points.resize(5);
      74         142 :               _weights.resize(5);
      75             : 
      76         142 :               _points[0](0) = Real(1)/6;  _points[0](1) = Real(1)/6;  _points[0](2) = Real(1)/2;  _weights[0] =  Real(3)/40;
      77         142 :               _points[1](0) = Real(1)/6;  _points[1](1) = Real(1)/6;  _points[1](2) = Real(1)/6;  _weights[1] =  Real(3)/40;
      78         142 :               _points[2](0) = Real(1)/4;  _points[2](1) = Real(1)/4;  _points[2](2) = Real(1)/4;  _weights[2] = -Real(2)/15;
      79         142 :               _points[3](0) = Real(1)/2;  _points[3](1) = Real(1)/6;  _points[3](2) = Real(1)/6;  _weights[3] =  Real(3)/40;
      80         142 :               _points[4](0) = Real(1)/6;  _points[4](1) = Real(1)/2;  _points[4](2) = Real(1)/6;  _weights[4] =  Real(3)/40;
      81         142 :               return;
      82             :             }
      83             : 
      84         142 :           case FOURTH:
      85             :           case FIFTH:
      86             :             {
      87         142 :               _points.resize(15);
      88         142 :               _weights.resize(15);
      89             : 
      90         142 :               _points[ 0](0) = Real(1)/8;  _points[ 0](1) = Real(3)/8;  _points[ 0](2) = Real(1)/8;  _weights[ 0] =  Real(16)/315;
      91         142 :               _points[ 1](0) = Real(1)/8;  _points[ 1](1) = Real(1)/8;  _points[ 1](2) = Real(5)/8;  _weights[ 1] =  Real(16)/315;
      92         142 :               _points[ 2](0) = Real(3)/8;  _points[ 2](1) = Real(1)/8;  _points[ 2](2) = Real(1)/8;  _weights[ 2] =  Real(16)/315;
      93         142 :               _points[ 3](0) = Real(1)/6;  _points[ 3](1) = Real(1)/2;  _points[ 3](2) = Real(1)/6;  _weights[ 3] = -Real(27)/280;
      94         142 :               _points[ 4](0) = Real(3)/8;  _points[ 4](1) = Real(1)/8;  _points[ 4](2) = Real(3)/8;  _weights[ 4] =  Real(16)/315;
      95         142 :               _points[ 5](0) = Real(1)/8;  _points[ 5](1) = Real(3)/8;  _points[ 5](2) = Real(3)/8;  _weights[ 5] =  Real(16)/315;
      96         142 :               _points[ 6](0) = Real(1)/6;  _points[ 6](1) = Real(1)/6;  _points[ 6](2) = Real(1)/6;  _weights[ 6] = -Real(27)/280;
      97         142 :               _points[ 7](0) = Real(1)/6;  _points[ 7](1) = Real(1)/6;  _points[ 7](2) = Real(1)/2;  _weights[ 7] = -Real(27)/280;
      98         142 :               _points[ 8](0) = Real(1)/8;  _points[ 8](1) = Real(1)/8;  _points[ 8](2) = Real(1)/8;  _weights[ 8] =  Real(16)/315;
      99         142 :               _points[ 9](0) = Real(1)/4;  _points[ 9](1) = Real(1)/4;  _points[ 9](2) = Real(1)/4;  _weights[ 9] =    Real(2)/45;
     100         142 :               _points[10](0) = Real(1)/8;  _points[10](1) = Real(5)/8;  _points[10](2) = Real(1)/8;  _weights[10] =  Real(16)/315;
     101         142 :               _points[11](0) = Real(1)/2;  _points[11](1) = Real(1)/6;  _points[11](2) = Real(1)/6;  _weights[11] = -Real(27)/280;
     102         142 :               _points[12](0) = Real(1)/8;  _points[12](1) = Real(1)/8;  _points[12](2) = Real(3)/8;  _weights[12] =  Real(16)/315;
     103         142 :               _points[13](0) = Real(5)/8;  _points[13](1) = Real(1)/8;  _points[13](2) = Real(1)/8;  _weights[13] =  Real(16)/315;
     104         142 :               _points[14](0) = Real(3)/8;  _points[14](1) = Real(3)/8;  _points[14](2) = Real(1)/8;  _weights[14] =  Real(16)/315;
     105         142 :               return;
     106             :             }
     107             : 
     108          71 :           case SIXTH:
     109             :           case SEVENTH:
     110             :             {
     111          71 :               _points.resize(35);
     112          71 :               _weights.resize(35);
     113             : 
     114          71 :               _points[ 0](0) = Real(3)/10;  _points[ 0](1) = Real(1)/10;  _points[ 0](2) = Real(3)/10;  _weights[ 0] = Real(3125)/72576;
     115          71 :               _points[ 1](0) =  Real(1)/6;  _points[ 1](1) =  Real(1)/2;  _points[ 1](2) =  Real(1)/6;  _weights[ 1] =   Real(243)/4480;
     116          71 :               _points[ 2](0) =  Real(1)/6;  _points[ 2](1) =  Real(1)/6;  _points[ 2](2) =  Real(1)/2;  _weights[ 2] =   Real(243)/4480;
     117          71 :               _points[ 3](0) =  Real(1)/2;  _points[ 3](1) = Real(1)/10;  _points[ 3](2) = Real(1)/10;  _weights[ 3] = Real(3125)/72576;
     118          71 :               _points[ 4](0) = Real(3)/10;  _points[ 4](1) = Real(1)/10;  _points[ 4](2) = Real(1)/10;  _weights[ 4] = Real(3125)/72576;
     119          71 :               _points[ 5](0) = Real(3)/10;  _points[ 5](1) = Real(3)/10;  _points[ 5](2) = Real(1)/10;  _weights[ 5] = Real(3125)/72576;
     120          71 :               _points[ 6](0) =  Real(1)/2;  _points[ 6](1) =  Real(1)/6;  _points[ 6](2) =  Real(1)/6;  _weights[ 6] =   Real(243)/4480;
     121          71 :               _points[ 7](0) = Real(3)/10;  _points[ 7](1) = Real(1)/10;  _points[ 7](2) =  Real(1)/2;  _weights[ 7] = Real(3125)/72576;
     122          71 :               _points[ 8](0) = Real(7)/10;  _points[ 8](1) = Real(1)/10;  _points[ 8](2) = Real(1)/10;  _weights[ 8] = Real(3125)/72576;
     123          71 :               _points[ 9](0) =  Real(3)/8;  _points[ 9](1) =  Real(1)/8;  _points[ 9](2) =  Real(1)/8;  _weights[ 9] =  -Real(256)/2835;
     124          71 :               _points[10](0) = Real(3)/10;  _points[10](1) = Real(3)/10;  _points[10](2) = Real(3)/10;  _weights[10] = Real(3125)/72576;
     125          71 :               _points[11](0) = Real(1)/10;  _points[11](1) =  Real(1)/2;  _points[11](2) = Real(3)/10;  _weights[11] = Real(3125)/72576;
     126          71 :               _points[12](0) = Real(1)/10;  _points[12](1) = Real(1)/10;  _points[12](2) = Real(7)/10;  _weights[12] = Real(3125)/72576;
     127          71 :               _points[13](0) =  Real(1)/2;  _points[13](1) = Real(1)/10;  _points[13](2) = Real(3)/10;  _weights[13] = Real(3125)/72576;
     128          71 :               _points[14](0) = Real(1)/10;  _points[14](1) = Real(7)/10;  _points[14](2) = Real(1)/10;  _weights[14] = Real(3125)/72576;
     129          71 :               _points[15](0) = Real(1)/10;  _points[15](1) =  Real(1)/2;  _points[15](2) = Real(1)/10;  _weights[15] = Real(3125)/72576;
     130          71 :               _points[16](0) =  Real(1)/6;  _points[16](1) =  Real(1)/6;  _points[16](2) =  Real(1)/6;  _weights[16] =   Real(243)/4480;
     131          71 :               _points[17](0) =  Real(3)/8;  _points[17](1) =  Real(1)/8;  _points[17](2) =  Real(3)/8;  _weights[17] =  -Real(256)/2835;
     132          71 :               _points[18](0) =  Real(1)/8;  _points[18](1) =  Real(1)/8;  _points[18](2) =  Real(5)/8;  _weights[18] =  -Real(256)/2835;
     133          71 :               _points[19](0) = Real(1)/10;  _points[19](1) = Real(1)/10;  _points[19](2) = Real(3)/10;  _weights[19] = Real(3125)/72576;
     134          71 :               _points[20](0) =  Real(1)/8;  _points[20](1) =  Real(3)/8;  _points[20](2) =  Real(3)/8;  _weights[20] =  -Real(256)/2835;
     135          71 :               _points[21](0) =  Real(5)/8;  _points[21](1) =  Real(1)/8;  _points[21](2) =  Real(1)/8;  _weights[21] =  -Real(256)/2835;
     136          71 :               _points[22](0) =  Real(1)/8;  _points[22](1) =  Real(5)/8;  _points[22](2) =  Real(1)/8;  _weights[22] =  -Real(256)/2835;
     137          71 :               _points[23](0) = Real(1)/10;  _points[23](1) = Real(3)/10;  _points[23](2) = Real(1)/10;  _weights[23] = Real(3125)/72576;
     138          71 :               _points[24](0) =  Real(1)/4;  _points[24](1) =  Real(1)/4;  _points[24](2) =  Real(1)/4;  _weights[24] =     -Real(8)/945;
     139          71 :               _points[25](0) =  Real(1)/8;  _points[25](1) =  Real(1)/8;  _points[25](2) =  Real(3)/8;  _weights[25] =  -Real(256)/2835;
     140          71 :               _points[26](0) =  Real(3)/8;  _points[26](1) =  Real(3)/8;  _points[26](2) =  Real(1)/8;  _weights[26] =  -Real(256)/2835;
     141          71 :               _points[27](0) =  Real(1)/8;  _points[27](1) =  Real(3)/8;  _points[27](2) =  Real(1)/8;  _weights[27] =  -Real(256)/2835;
     142          71 :               _points[28](0) = Real(1)/10;  _points[28](1) = Real(3)/10;  _points[28](2) =  Real(1)/2;  _weights[28] = Real(3125)/72576;
     143          71 :               _points[29](0) = Real(3)/10;  _points[29](1) =  Real(1)/2;  _points[29](2) = Real(1)/10;  _weights[29] = Real(3125)/72576;
     144          71 :               _points[30](0) = Real(1)/10;  _points[30](1) = Real(1)/10;  _points[30](2) =  Real(1)/2;  _weights[30] = Real(3125)/72576;
     145          71 :               _points[31](0) =  Real(1)/2;  _points[31](1) = Real(3)/10;  _points[31](2) = Real(1)/10;  _weights[31] = Real(3125)/72576;
     146          71 :               _points[32](0) =  Real(1)/8;  _points[32](1) =  Real(1)/8;  _points[32](2) =  Real(1)/8;  _weights[32] =  -Real(256)/2835;
     147          71 :               _points[33](0) = Real(1)/10;  _points[33](1) = Real(3)/10;  _points[33](2) = Real(3)/10;  _weights[33] = Real(3125)/72576;
     148          71 :               _points[34](0) = Real(1)/10;  _points[34](1) = Real(1)/10;  _points[34](2) = Real(1)/10;  _weights[34] = Real(3125)/72576;
     149          71 :               return;
     150             :             }
     151             : 
     152        6110 :           default:
     153             :             {
     154             :               // Untested above _order=23 but should work...
     155        6110 :               gm_rule(get_order()/2, /*dim=*/3);
     156        6110 :               return;
     157             :             }
     158             :           } // end switch (order)
     159             :       } // end case TET
     160             : 
     161           0 :     default:
     162           0 :       libmesh_error_msg("ERROR: Unsupported element type: " << Utility::enum_to_string(_type));
     163             :     } // end switch (_type)
     164             : }
     165             : 
     166             : } // namespace libMesh