Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Integrate/Quadrature.cpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Quadrature coordinates and weights for numerical integration
       9                 :            :   \details   This file contains functions that provide Gauss quadrature
      10                 :            :     coordinates and weights for numerical integration.
      11                 :            : */
      12                 :            : // *****************************************************************************
      13                 :            : 
      14                 :            : #include "Quadrature.hpp"
      15                 :            : 
      16                 :            : void
      17                 :   20429716 : tk::GaussQuadratureTet( const std::size_t NG,
      18                 :            :                         std::array< std::vector< real >, 3>& coordgp,
      19                 :            :                         std::vector< real >& wgp )
      20                 :            : // *****************************************************************************
      21                 :            : //! Initialize Gaussian quadrature points locations and weights for a tetrahedron
      22                 :            : //! \param[in] NG number of quadrature points
      23                 :            : //! \param[in,out] coordgp 3 spatial coordinates of quadrature points
      24                 :            : //! \param[in,out] wgp Weights of quadrature points
      25                 :            : // *****************************************************************************
      26                 :            : {
      27 [ -  + ][ -  - ]:   20429716 :   Assert( coordgp[0].size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
      28 [ -  + ][ -  - ]:   20429716 :   Assert( coordgp[1].size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
      29 [ -  + ][ -  - ]:   20429716 :   Assert( coordgp[1].size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
      30 [ -  + ][ -  - ]:   20429716 :   Assert( wgp.size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
      31                 :            : 
      32 [ +  + ][ +  + ]:   20429716 :   switch( NG )
                 [ +  - ]
      33                 :            :   {
      34                 :   10311825 :     case 1:
      35                 :            :     {
      36                 :   10311825 :       coordgp[0][0] = 0.25;
      37                 :   10311825 :       coordgp[1][0] = 0.25;
      38                 :   10311825 :       coordgp[2][0] = 0.25;
      39                 :   10311825 :       wgp[0]        = 1.0;
      40                 :            :     }
      41                 :   10311825 :     break;
      42                 :            : 
      43                 :     500647 :     case 4:
      44                 :            :     {
      45                 :     500647 :       const real a1 = 0.5854101966249685;
      46                 :     500647 :       const real a2 = 0.1381966011250105;
      47                 :            : 
      48                 :     500647 :       coordgp[0][0] = a2;
      49                 :     500647 :       coordgp[1][0] = a2;
      50                 :     500647 :       coordgp[2][0] = a2;
      51                 :     500647 :       wgp[0]        = 0.25;
      52                 :            : 
      53                 :     500647 :       coordgp[0][1] = a1;
      54                 :     500647 :       coordgp[1][1] = a2;
      55                 :     500647 :       coordgp[2][1] = a2;
      56                 :     500647 :       wgp[1]        = 0.25;
      57                 :            : 
      58                 :     500647 :       coordgp[0][2] = a2;
      59                 :     500647 :       coordgp[1][2] = a1;
      60                 :     500647 :       coordgp[2][2] = a2;
      61                 :     500647 :       wgp[2]        = 0.25;
      62                 :            : 
      63                 :     500647 :       coordgp[0][3] = a2;
      64                 :     500647 :       coordgp[1][3] = a2;
      65                 :     500647 :       coordgp[2][3] = a1;
      66                 :     500647 :       wgp[3]        = 0.25;
      67                 :            :     }
      68                 :     500647 :     break;
      69                 :            : 
      70                 :    7301898 :     case 5:
      71                 :            :     {
      72                 :    7301898 :       coordgp[0][0] = 0.25;
      73                 :    7301898 :       coordgp[1][0] = 0.25;
      74                 :    7301898 :       coordgp[2][0] = 0.25;
      75                 :    7301898 :       wgp[0]        = -12.0/15.0;
      76                 :            : 
      77                 :    7301898 :       coordgp[0][1] = 1.0/6.0;
      78                 :    7301898 :       coordgp[1][1] = 1.0/6.0;
      79                 :    7301898 :       coordgp[2][1] = 1.0/6.0;
      80                 :    7301898 :       wgp[1]        = 9.0/20.0;
      81                 :            : 
      82                 :    7301898 :       coordgp[0][2] = 0.5;
      83                 :    7301898 :       coordgp[1][2] = 1.0/6.0;
      84                 :    7301898 :       coordgp[2][2] = 1.0/6.0;
      85                 :    7301898 :       wgp[2]        = 9.0/20.0;
      86                 :            : 
      87                 :    7301898 :       coordgp[0][3] = 1.0/6.0;
      88                 :    7301898 :       coordgp[1][3] = 0.5;
      89                 :    7301898 :       coordgp[2][3] = 1.0/6.0;
      90                 :    7301898 :       wgp[3]        = 9.0/20.0;
      91                 :            : 
      92                 :    7301898 :       coordgp[0][4] = 1.0/6.0;
      93                 :    7301898 :       coordgp[1][4] = 1.0/6.0;
      94                 :    7301898 :       coordgp[2][4] = 0.5;
      95                 :    7301898 :       wgp[4]        = 9.0/20.0;
      96                 :            :     }
      97                 :    7301898 :     break;
      98                 :            : 
      99                 :    2211205 :     case 11:
     100                 :            :     {
     101                 :    2211205 :       const tk::real c1 = 0.3994035761667992;
     102                 :    2211205 :       const tk::real c2 = 0.1005964238332008;
     103                 :    2211205 :       const tk::real c3 = 343.0 / 7500.0;
     104                 :    2211205 :       const tk::real c4 = 56.0 / 375.0;
     105                 :            : 
     106                 :    2211205 :       coordgp[0][0] = 0.25;
     107                 :    2211205 :       coordgp[1][0] = 0.25;
     108                 :    2211205 :       coordgp[2][0] = 0.25;
     109                 :    2211205 :       wgp[0]        = -148.0/1875.0;
     110                 :            : 
     111                 :    2211205 :       coordgp[0][1] = 11.0/14.0;
     112                 :    2211205 :       coordgp[1][1] = 1.0/14.0;
     113                 :    2211205 :       coordgp[2][1] = 1.0/14.0;
     114                 :    2211205 :       wgp[1]        = c3;
     115                 :            : 
     116                 :    2211205 :       coordgp[0][2] = 1.0/14.0;
     117                 :    2211205 :       coordgp[1][2] = 11.0/14.0;
     118                 :    2211205 :       coordgp[2][2] = 1.0/14.0;
     119                 :    2211205 :       wgp[2]        = c3;
     120                 :            : 
     121                 :    2211205 :       coordgp[0][3] = 1.0/14.0;
     122                 :    2211205 :       coordgp[1][3] = 1.0/14.0;
     123                 :    2211205 :       coordgp[2][3] = 11.0/14.0;
     124                 :    2211205 :       wgp[3]        = c3;
     125                 :            : 
     126                 :    2211205 :       coordgp[0][4] = 1.0/14.0;
     127                 :    2211205 :       coordgp[1][4] = 1.0/14.0;
     128                 :    2211205 :       coordgp[2][4] = 1.0/14.0;
     129                 :    2211205 :       wgp[4]        = c3;
     130                 :            : 
     131                 :    2211205 :       coordgp[0][5] = c1;
     132                 :    2211205 :       coordgp[1][5] = c1;
     133                 :    2211205 :       coordgp[2][5] = c2;
     134                 :    2211205 :       wgp[5]        = c4;
     135                 :            : 
     136                 :    2211205 :       coordgp[0][6] = c1;
     137                 :    2211205 :       coordgp[1][6] = c2;
     138                 :    2211205 :       coordgp[2][6] = c1;
     139                 :    2211205 :       wgp[6]        = c4;
     140                 :            : 
     141                 :    2211205 :       coordgp[0][7] = c1;
     142                 :    2211205 :       coordgp[1][7] = c2;
     143                 :    2211205 :       coordgp[2][7] = c2;
     144                 :    2211205 :       wgp[7]        = c4;
     145                 :            : 
     146                 :    2211205 :       coordgp[0][8] = c2;
     147                 :    2211205 :       coordgp[1][8] = c1;
     148                 :    2211205 :       coordgp[2][8] = c1;
     149                 :    2211205 :       wgp[8]        = c4;
     150                 :            : 
     151                 :    2211205 :       coordgp[0][9] = c2;
     152                 :    2211205 :       coordgp[1][9] = c1;
     153                 :    2211205 :       coordgp[2][9] = c2;
     154                 :    2211205 :       wgp[9]        = c4;
     155                 :            : 
     156                 :    2211205 :       coordgp[0][10] = c2;
     157                 :    2211205 :       coordgp[1][10] = c2;
     158                 :    2211205 :       coordgp[2][10] = c1;
     159                 :    2211205 :       wgp[10]        = c4;
     160                 :            :     }
     161                 :    2211205 :     break;
     162                 :            : 
     163                 :     104141 :     case 14:
     164                 :            :     {
     165                 :     104141 :       const tk::real a = 0.0673422422100983;
     166                 :     104141 :       const tk::real b = 0.3108859192633005;
     167                 :     104141 :       const tk::real c = 0.7217942490673264;
     168                 :     104141 :       const tk::real d = 0.0927352503108912;
     169                 :     104141 :       const tk::real e = 0.4544962958743506;
     170                 :     104141 :       const tk::real f = 0.0455037041256494;
     171                 :     104141 :       const tk::real p = 0.1126879257180162;
     172                 :     104141 :       const tk::real q = 0.0734930431163619;
     173                 :     104141 :       const tk::real r = 0.0425460207770812;
     174                 :            : 
     175                 :     104141 :       coordgp[0][0] = a;
     176                 :     104141 :       coordgp[1][0] = b;
     177                 :     104141 :       coordgp[2][0] = b;
     178                 :     104141 :       wgp[0]        = p;
     179                 :            : 
     180                 :     104141 :       coordgp[0][1] = b;
     181                 :     104141 :       coordgp[1][1] = a;
     182                 :     104141 :       coordgp[2][1] = b;
     183                 :     104141 :       wgp[1]        = p;
     184                 :            : 
     185                 :     104141 :       coordgp[0][2] = b;
     186                 :     104141 :       coordgp[1][2] = b;
     187                 :     104141 :       coordgp[2][2] = a;
     188                 :     104141 :       wgp[2]        = p;
     189                 :            : 
     190                 :     104141 :       coordgp[0][3] = b;
     191                 :     104141 :       coordgp[1][3] = b;
     192                 :     104141 :       coordgp[2][3] = b;
     193                 :     104141 :       wgp[3]        = p;
     194                 :            : 
     195                 :     104141 :       coordgp[0][4] = c;
     196                 :     104141 :       coordgp[1][4] = d;
     197                 :     104141 :       coordgp[2][4] = d;
     198                 :     104141 :       wgp[4]        = q;
     199                 :            : 
     200                 :     104141 :       coordgp[0][5] = d;
     201                 :     104141 :       coordgp[1][5] = c;
     202                 :     104141 :       coordgp[2][5] = d;
     203                 :     104141 :       wgp[5]        = q;
     204                 :            : 
     205                 :     104141 :       coordgp[0][6] = d;
     206                 :     104141 :       coordgp[1][6] = d;
     207                 :     104141 :       coordgp[2][6] = c;
     208                 :     104141 :       wgp[6]        = q;
     209                 :            : 
     210                 :     104141 :       coordgp[0][7] = d;
     211                 :     104141 :       coordgp[1][7] = d;
     212                 :     104141 :       coordgp[2][7] = d;
     213                 :     104141 :       wgp[7]        = q;
     214                 :            : 
     215                 :     104141 :       coordgp[0][8] = e;
     216                 :     104141 :       coordgp[1][8] = e;
     217                 :     104141 :       coordgp[2][8] = f;
     218                 :     104141 :       wgp[8]        = r;
     219                 :            : 
     220                 :     104141 :       coordgp[0][9] = e;
     221                 :     104141 :       coordgp[1][9] = f;
     222                 :     104141 :       coordgp[2][9] = e;
     223                 :     104141 :       wgp[9]        = r;
     224                 :            : 
     225                 :     104141 :       coordgp[0][10] = e;
     226                 :     104141 :       coordgp[1][10] = f;
     227                 :     104141 :       coordgp[2][10] = f;
     228                 :     104141 :       wgp[10]        = r;
     229                 :            : 
     230                 :     104141 :       coordgp[0][11] = f;
     231                 :     104141 :       coordgp[1][11] = e;
     232                 :     104141 :       coordgp[2][11] = e;
     233                 :     104141 :       wgp[11]        = r;
     234                 :            : 
     235                 :     104141 :       coordgp[0][12] = f;
     236                 :     104141 :       coordgp[1][12] = e;
     237                 :     104141 :       coordgp[2][12] = f;
     238                 :     104141 :       wgp[12]        = r;
     239                 :            : 
     240                 :     104141 :       coordgp[0][13] = f;
     241                 :     104141 :       coordgp[1][13] = f;
     242                 :     104141 :       coordgp[2][13] = e;
     243                 :     104141 :       wgp[13]        = r;
     244                 :            :     }
     245                 :     104141 :     break;
     246                 :            :   }
     247                 :   20429716 : }
     248                 :            : 
     249                 :            : void
     250                 :   50815992 : tk::GaussQuadratureTri( const std::size_t NG,
     251                 :            :                         std::array< std::vector< real >, 2>& coordgp,
     252                 :            :                         std::vector< real >& wgp )
     253                 :            : // *****************************************************************************
     254                 :            : //! Initialize Gaussian quadrature points locations and weights for a triangle
     255                 :            : //! \param[in] NG number of quadrature points
     256                 :            : //! \param[in,out] coordgp 2 spatial coordinates of quadrature points
     257                 :            : //! \param[in,out] wgp Weights of quadrature points
     258                 :            : // *****************************************************************************
     259                 :            : {
     260 [ -  + ][ -  - ]:   50815992 :   Assert( coordgp[0].size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
     261 [ -  + ][ -  - ]:   50815992 :   Assert( coordgp[1].size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
     262 [ -  + ][ -  - ]:   50815992 :   Assert( wgp.size() == NG, "Size mismatch" );
         [ -  - ][ -  - ]
     263                 :            :   
     264 [ +  + ][ -  + ]:   50815992 :   switch( NG )
                    [ - ]
     265                 :            :   {
     266                 :   28734708 :     case 1:
     267                 :   28734708 :       coordgp[0][0] = 1.0/3.0;
     268                 :   28734708 :       coordgp[1][0] = 1.0/3.0;
     269                 :   28734708 :       wgp[0]        = 1.0;
     270                 :   28734708 :       break;
     271                 :            : 
     272                 :   18430054 :     case 3:
     273                 :   18430054 :       coordgp[0][0] = 2.0/3.0;
     274                 :   18430054 :       coordgp[1][0] = 1.0/6.0;
     275                 :   18430054 :       wgp[0]        = 1.0/3.0;
     276                 :            : 
     277                 :   18430054 :       coordgp[0][1] = 1.0/6.0;
     278                 :   18430054 :       coordgp[1][1] = 2.0/3.0;
     279                 :   18430054 :       wgp[1]        = 1.0/3.0;
     280                 :            : 
     281                 :   18430054 :       coordgp[0][2] = 1.0/6.0;
     282                 :   18430054 :       coordgp[1][2] = 1.0/6.0;
     283                 :   18430054 :       wgp[2]        = 1.0/3.0;
     284                 :   18430054 :       break;
     285                 :            : 
     286                 :          0 :     case 4:
     287                 :          0 :       coordgp[0][0] = 1.0/3.0;
     288                 :          0 :       coordgp[1][0] = 1.0/3.0;
     289                 :          0 :       wgp[0]        = -27.0/48.0;
     290                 :            : 
     291                 :          0 :       coordgp[0][1] = 1.0/5.0;
     292                 :          0 :       coordgp[1][1] = 1.0/5.0;
     293                 :          0 :       wgp[1]        = 25.0/48.0;
     294                 :            : 
     295                 :          0 :       coordgp[0][2] = 3.0/5.0;
     296                 :          0 :       coordgp[1][2] = 1.0/5.0;
     297                 :          0 :       wgp[2]        = 25.0/48.0;
     298                 :            : 
     299                 :          0 :       coordgp[0][3] = 1.0/5.0;
     300                 :          0 :       coordgp[1][3] = 3.0/5.0;
     301                 :          0 :       wgp[3]        = 25.0/48.0;
     302                 :          0 :       break;
     303                 :            : 
     304                 :    3651230 :     case 6:
     305                 :    3651230 :       const tk::real c1 = 0.816847572980459;
     306                 :    3651230 :       const tk::real c2 = 0.091576213509771;
     307                 :    3651230 :       const tk::real c3 = 0.091576213509771;
     308                 :    3651230 :       const tk::real c4 = 0.108103018168070;
     309                 :    3651230 :       const tk::real c5 = 0.445948490915965;
     310                 :    3651230 :       const tk::real c6 = 0.445948490915965;
     311                 :    3651230 :       const tk::real w1 = 0.109951743655322;
     312                 :    3651230 :       const tk::real w2 = 0.223381589678011;
     313                 :            : 
     314                 :    3651230 :       coordgp[0][0] = c1;
     315                 :    3651230 :       coordgp[1][0] = c2;
     316                 :    3651230 :       wgp[0]        = w1;
     317                 :            : 
     318                 :    3651230 :       coordgp[0][1] = c2;
     319                 :    3651230 :       coordgp[1][1] = c3;
     320                 :    3651230 :       wgp[1]        = w1;
     321                 :            : 
     322                 :    3651230 :       coordgp[0][2] = c3;
     323                 :    3651230 :       coordgp[1][2] = c1;
     324                 :    3651230 :       wgp[2]        = w1;
     325                 :            : 
     326                 :    3651230 :       coordgp[0][3] = c4;
     327                 :    3651230 :       coordgp[1][3] = c5;
     328                 :    3651230 :       wgp[3]        = w2;
     329                 :            : 
     330                 :    3651230 :       coordgp[0][4] = c5;
     331                 :    3651230 :       coordgp[1][4] = c6;
     332                 :    3651230 :       wgp[4]        = w2;
     333                 :            : 
     334                 :    3651230 :       coordgp[0][5] = c6;
     335                 :    3651230 :       coordgp[1][5] = c4;
     336                 :    3651230 :       wgp[5]        = w2;
     337                 :    3651230 :       break;
     338                 :            :   }
     339                 :   50815992 : }