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 "XFEMCutElem.h"
      11             : 
      12             : #include "libmesh/mesh.h"
      13             : #include "libmesh/elem.h"
      14             : #include "libmesh/node.h"
      15             : #include "libmesh/quadrature.h"
      16             : #include "libmesh/fe.h"
      17             : #include "libmesh/enum_fe_family.h"
      18             : 
      19             : #include "EFANode.h"
      20             : #include "EFAElement.h"
      21             : 
      22       11858 : XFEMCutElem::XFEMCutElem(Elem * elem, unsigned int n_qpoints, unsigned int n_sides)
      23       23716 :   : _n_nodes(elem->n_nodes()),
      24       11858 :     _n_qpoints(n_qpoints),
      25       11858 :     _n_sides(n_sides),
      26       11858 :     _nodes(_n_nodes, nullptr),
      27       11858 :     _elem_side_area(n_sides),
      28       11858 :     _physical_areafrac(n_sides),
      29       11858 :     _have_weights(false),
      30       11858 :     _have_face_weights(n_sides),
      31       23716 :     _new_face_weights(n_sides)
      32             : {
      33       79534 :   for (unsigned int i = 0; i < _n_nodes; ++i)
      34       67676 :     _nodes[i] = elem->node_ptr(i);
      35             : 
      36       60730 :   for (unsigned int i = 0; i < _n_sides; ++i)
      37             :   {
      38             :     _have_face_weights[i] = false;
      39       48872 :     _physical_areafrac[i] = 1.0;
      40             :     // TODO: this line does not support RZ/RSpherical geoms
      41       48872 :     _elem_side_area[i] = elem->side_ptr(i)->volume();
      42             :   }
      43             : 
      44             :   // TODO: this line does not support RZ/RSpherical geoms
      45       11858 :   _elem_volume = elem->volume();
      46       11858 : }
      47             : 
      48       23716 : XFEMCutElem::~XFEMCutElem() {}
      49             : 
      50             : Real
      51      135527 : XFEMCutElem::getPhysicalVolumeFraction() const
      52             : {
      53      135527 :   return _physical_volfrac;
      54             : }
      55             : 
      56             : Real
      57          78 : XFEMCutElem::getPhysicalFaceAreaFraction(unsigned int side) const
      58             : {
      59          78 :   return _physical_areafrac[side];
      60             : }
      61             : 
      62             : void
      63      626171 : XFEMCutElem::getWeightMultipliers(MooseArray<Real> & weights,
      64             :                                   QBase * qrule,
      65             :                                   Xfem::XFEM_QRULE xfem_qrule,
      66             :                                   const MooseArray<Point> & q_points)
      67             : {
      68      626171 :   if (!_have_weights)
      69        8021 :     computeXFEMWeights(qrule, xfem_qrule, q_points);
      70             : 
      71      626171 :   weights.resize(_new_weights.size());
      72     5063807 :   for (unsigned int qp = 0; qp < _new_weights.size(); ++qp)
      73     4437636 :     weights[qp] = _new_weights[qp];
      74      626171 : }
      75             : 
      76             : void
      77         290 : XFEMCutElem::getFaceWeightMultipliers(MooseArray<Real> & face_weights,
      78             :                                       QBase * qrule,
      79             :                                       Xfem::XFEM_QRULE xfem_qrule,
      80             :                                       const MooseArray<Point> & q_points,
      81             :                                       unsigned int side)
      82             : {
      83         290 :   if (!_have_face_weights[side])
      84          78 :     computeXFEMFaceWeights(qrule, xfem_qrule, q_points, side);
      85             : 
      86         290 :   face_weights.resize(_new_face_weights[side].size());
      87        1150 :   for (unsigned int qp = 0; qp < _new_face_weights[side].size(); ++qp)
      88         860 :     face_weights[qp] = _new_face_weights[side][qp];
      89         290 : }
      90             : 
      91             : void
      92          78 : XFEMCutElem::computeXFEMFaceWeights(QBase * qrule,
      93             :                                     Xfem::XFEM_QRULE /*xfem_qrule*/,
      94             :                                     const MooseArray<Point> & /*q_points*/,
      95             :                                     unsigned int side)
      96             : {
      97          78 :   _new_face_weights[side].clear();
      98          78 :   _new_face_weights[side].resize(qrule->n_points(), 1.0);
      99             : 
     100          78 :   computePhysicalFaceAreaFraction(side);
     101          78 :   Real surffrac = getPhysicalFaceAreaFraction(side);
     102         354 :   for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
     103         276 :     _new_face_weights[side][qp] = surffrac;
     104             : 
     105             :   _have_face_weights[side] = true;
     106          78 : }
     107             : 
     108             : void
     109        8021 : XFEMCutElem::computeXFEMWeights(QBase * qrule,
     110             :                                 Xfem::XFEM_QRULE xfem_qrule,
     111             :                                 const MooseArray<Point> & q_points)
     112             : {
     113        8021 :   _new_weights.clear();
     114             : 
     115        8021 :   _new_weights.resize(qrule->n_points(), 1.0);
     116             : 
     117        8021 :   switch (xfem_qrule)
     118             :   {
     119        7673 :     case Xfem::VOLFRAC:
     120             :     {
     121        7673 :       computePhysicalVolumeFraction();
     122        7673 :       Real volfrac = getPhysicalVolumeFraction();
     123       68719 :       for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
     124       61046 :         _new_weights[qp] = volfrac;
     125             :       break;
     126             :     }
     127             :     case Xfem::MOMENT_FITTING:
     128             :     {
     129             :       // These are the coordinates in parametric coordinates
     130         348 :       _qp_points = qrule->get_points();
     131         348 :       _qp_weights = qrule->get_weights();
     132             : 
     133         348 :       computeMomentFittingWeights();
     134             : 
     135             :       // Blend weights from moment fitting and volume fraction to avoid negative weights
     136             :       Real alpha = 1.0;
     137         348 :       if (*std::min_element(_new_weights.begin(), _new_weights.end()) < 0.0)
     138             :       {
     139             :         // One or more of the weights computed by moment fitting is negative.
     140             :         // Blend moment and volume fraction weights to keep them nonnegative.
     141             :         // Find the largest value of alpha that will keep all of the weights nonnegative.
     142         906 :         for (unsigned int i = 0; i < _n_qpoints; ++i)
     143             :         {
     144         732 :           const Real denominator = _physical_volfrac - _new_weights[i];
     145         732 :           if (denominator > 0.0) // Negative values would give a negative value for alpha, which
     146             :                                  // must be between 0 and 1.
     147             :           {
     148         420 :             const Real alpha_i = _physical_volfrac / denominator;
     149         420 :             if (alpha_i < alpha)
     150             :               alpha = alpha_i;
     151             :           }
     152             :         }
     153         906 :         for (unsigned int i = 0; i < _n_qpoints; ++i)
     154             :         {
     155         732 :           _new_weights[i] = alpha * _new_weights[i] + (1.0 - alpha) * _physical_volfrac;
     156         732 :           if (_new_weights[i] < 0.0) // We can end up with small (roundoff) negative weights
     157          96 :             _new_weights[i] = 0.0;
     158             :         }
     159             :       }
     160             :       break;
     161             :     }
     162             :     case Xfem::DIRECT: // remove q-points outside the partial element's physical domain
     163             :     {
     164             :       bool nonzero = false;
     165           0 :       for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
     166             :       {
     167             :         // q_points contains quadrature point locations in physical coordinates
     168           0 :         if (isPointPhysical(q_points[qp]))
     169             :         {
     170             :           nonzero = true;
     171           0 :           _new_weights[qp] = 1.0;
     172             :         }
     173             :         else
     174           0 :           _new_weights[qp] = 0.0;
     175             :       }
     176           0 :       if (!nonzero)
     177             :       {
     178             :         // Set the weights to a small value (1e-3) to avoid having DOFs
     179             :         // with zeros on the diagonal, which occurs for nodes that are
     180             :         // connected to no physical material.
     181           0 :         for (unsigned qp = 0; qp < qrule->n_points(); ++qp)
     182           0 :           _new_weights[qp] = 1e-3;
     183             :       }
     184             :       break;
     185             :     }
     186           0 :     default:
     187           0 :       mooseError("Undefined option for XFEM_QRULE");
     188             :   }
     189        8021 :   _have_weights = true;
     190        8021 : }
     191             : 
     192             : bool
     193      177040 : XFEMCutElem::isPointPhysical(const Point & p) const
     194             : {
     195             :   // determine whether the point is inside the physical domain of a partial element
     196             :   bool physical_flag = true;
     197      177040 :   unsigned int n_cut_planes = numCutPlanes();
     198      268386 :   for (unsigned int plane_id = 0; plane_id < n_cut_planes; ++plane_id)
     199             :   {
     200      177040 :     Point origin = getCutPlaneOrigin(plane_id);
     201      177040 :     Point normal = getCutPlaneNormal(plane_id);
     202             :     Point origin2qp = p - origin;
     203      177040 :     if (origin2qp * normal > 0.0)
     204             :     {
     205             :       physical_flag = false; // Point outside pysical domain
     206       85694 :       break;
     207             :     }
     208             :   }
     209      177040 :   return physical_flag;
     210             : }