Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PolygonalMeshGenerationUtils.h"
      11             : 
      12             : Real
      13        2329 : PolygonalMeshGenerationUtils::dummyTRI6VolCalculator(const std::pair<Real, Real> & azi_pair)
      14             : {
      15             :   // The algorithm is copied from libMesh's face_tri6.C
      16             :   // Original license is LGPL so it can be used here.
      17             :   const Real & azi_1 = azi_pair.first;
      18             :   const Real & azi_2 = azi_pair.second;
      19             : 
      20             :   Point x0 = Point(0.0, 0.0, 0.0), x1 = Point(1.0, 0.0, 0.0),
      21        2329 :         x2 = Point(std::cos(azi_1 + azi_2), std::sin(azi_1 + azi_2), 0.0),
      22        2329 :         x3 = Point(0.5, 0.0, 0.0), x4 = Point(std::cos(azi_1), std::sin(azi_1), 0.0),
      23        2329 :         x5 = Point(std::cos(azi_1 + azi_2) / 2.0, std::sin(azi_1 + azi_2) / 2.0, 0.0);
      24             : 
      25             :   // Construct constant data vectors.
      26             :   // \vec{x}_{\xi}  = \vec{a1}*xi + \vec{b1}*eta + \vec{c1}
      27             :   // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}*eta + \vec{c2}
      28             :   Point a1 = 4 * x0 + 4 * x1 - 8 * x3, b1 = 4 * x0 - 4 * x3 + 4 * x4 - 4 * x5,
      29             :         c1 = -3 * x0 - 1 * x1 + 4 * x3, b2 = 4 * x0 + 4 * x2 - 8 * x5,
      30             :         c2 = -3 * x0 - 1 * x2 + 4 * x5;
      31             : 
      32             :   // 7-point rule, exact for quintics.
      33             :   const unsigned int N = 7;
      34             : 
      35             :   // Parameters of the quadrature rule
      36             :   const static Real w1 = Real(31) / 480 + Real(std::sqrt(15.0L) / 2400),
      37             :                     w2 = Real(31) / 480 - Real(std::sqrt(15.0L) / 2400),
      38             :                     q1 = Real(2) / 7 + Real(std::sqrt(15.0L) / 21),
      39             :                     q2 = Real(2) / 7 - Real(std::sqrt(15.0L) / 21);
      40             : 
      41             :   const static Real xi[N] = {Real(1) / 3, q1, q1, 1 - 2 * q1, q2, q2, 1 - 2 * q2};
      42             :   const static Real eta[N] = {Real(1) / 3, q1, 1 - 2 * q1, q1, q2, 1 - 2 * q2, q2};
      43             :   const static Real wts[N] = {Real(9) / 80, w1, w1, w1, w2, w2, w2};
      44             : 
      45             :   // Approximate the area with quadrature
      46             :   Real vol = 0.;
      47       18632 :   for (unsigned int q = 0; q < N; ++q)
      48       16303 :     vol += wts[q] * cross_norm(xi[q] * a1 + eta[q] * b1 + c1, xi[q] * b1 + eta[q] * b2 + c2);
      49             : 
      50        2329 :   return vol;
      51             : }
      52             : 
      53             : Real
      54        2471 : PolygonalMeshGenerationUtils::radiusCorrectionFactor(const std::vector<Real> & azimuthal_list,
      55             :                                                      const bool full_circle,
      56             :                                                      const unsigned int order,
      57             :                                                      const bool is_first_value_vertex)
      58             : {
      59             :   Real tmp_acc = 0.0;
      60             :   Real tmp_acc_azi = 0.0;
      61        2471 :   if (order == 1)
      62             :   {
      63             :     // A vector to collect the azimuthal intervals of all EDGE2 side elements
      64             :     std::vector<Real> azi_interval_list;
      65       53308 :     for (unsigned int i = 1; i < azimuthal_list.size(); i++)
      66       50986 :       azi_interval_list.push_back((azimuthal_list[i] - azimuthal_list[i - 1]) / 180.0 * M_PI);
      67        2322 :     if (full_circle)
      68        2313 :       azi_interval_list.push_back((azimuthal_list.front() + 360.0 - azimuthal_list.back()) / 180.0 *
      69             :                                   M_PI);
      70             :     // summation of triangles S = 0.5 * r * r * Sigma_i [sin (azi_i)]
      71             :     // Circle area S_c = pi * r_0 * r_0
      72             :     // r = sqrt{2 * pi / Sigma_i [sin (azi_i)]} * r_0
      73       55621 :     for (const auto i : index_range(azi_interval_list))
      74             :     {
      75       53299 :       tmp_acc += std::sin(azi_interval_list[i]);
      76       53299 :       tmp_acc_azi += azi_interval_list[i];
      77             :     }
      78        2322 :   }
      79             :   else // order == 2
      80             :   {
      81             :     // A vector to collect the pairs of azimuthal intervals of all EDGE3 side elements
      82             :     // The midpoint divides the interval into two parts, which are stored as first and second of
      83             :     // each pair; the sum of the two parts is the full interval of the EDGE3 element
      84             :     std::vector<std::pair<Real, Real>> azi_interval_list;
      85        2478 :     for (unsigned int i = is_first_value_vertex ? 2 : 3; i < azimuthal_list.size(); i += 2)
      86             :       azi_interval_list.push_back(
      87        2180 :           std::make_pair((azimuthal_list[i - 1] - azimuthal_list[i - 2]) / 180.0 * M_PI,
      88        2180 :                          (azimuthal_list[i] - azimuthal_list[i - 1]) / 180.0 * M_PI));
      89             :     // If it is not a full circle, the first value should always belong to a vertex
      90         149 :     if (full_circle)
      91             :     {
      92         149 :       if (is_first_value_vertex)
      93         132 :         azi_interval_list.push_back(std::make_pair(
      94         132 :             (azimuthal_list.back() - azimuthal_list[azimuthal_list.size() - 2]) / 180.0 * M_PI,
      95         132 :             (azimuthal_list.front() + 360.0 - azimuthal_list.back()) / 180.0 * M_PI));
      96             :       else
      97             :         azi_interval_list.push_back(
      98          17 :             std::make_pair((azimuthal_list.front() + 360.0 - azimuthal_list.back()) / 180.0 * M_PI,
      99          17 :                            (azimuthal_list[1] - azimuthal_list.front()) / 180.0 * M_PI));
     100             :     }
     101             :     // Use the libMesh TRI6 element volume algorithm
     102        2478 :     for (const auto i : index_range(azi_interval_list))
     103             :     {
     104        2329 :       tmp_acc += 2.0 * dummyTRI6VolCalculator(azi_interval_list[i]);
     105        2329 :       tmp_acc_azi += azi_interval_list[i].first + azi_interval_list[i].second;
     106             :     }
     107         149 :   }
     108        2471 :   return std::sqrt(tmp_acc_azi / tmp_acc);
     109             : }