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             : #ifndef LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
      19             : #define LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H
      20             : 
      21             : // libMesh includes
      22             : #include "libmesh/enum_fe_family.h"
      23             : #include "libmesh/fem_function_base.h"
      24             : #include "libmesh/fem_system.h"
      25             : #include "libmesh/libmesh_common.h"
      26             : 
      27             : // C++ includes
      28             : #include <map>
      29             : #include <memory>
      30             : 
      31             : namespace libMesh
      32             : {
      33             : 
      34             : // FEMSystem, TimeSolver and  NewtonSolver will handle most tasks,
      35             : // but we must specify element residuals
      36          96 : class VariationalSmootherSystem : public libMesh::FEMSystem
      37             : {
      38             : /**
      39             :  * This is an FEMSystem to solve the optimization probelem posed by the
      40             :  * VariationalMeshSmoother class.
      41             :  *
      42             :  * The residual is coded as the gradient of the distortion-dilation metric, and
      43             :  * the jacobian as analytically coded as the Hessian of the metric.
      44             :  *
      45             :  * The nodes of the system mesh are updated during the solve.
      46             :  */
      47             : public:
      48         852 :   VariationalSmootherSystem(libMesh::EquationSystems & es,
      49             :                 const std::string & name,
      50             :                 const unsigned int number)
      51         852 :   : libMesh::FEMSystem(es, name, number),
      52         804 :     _epsilon_squared(1e-10),
      53         804 :     _ref_vol(1.),
      54         876 :     _dilation_weight(0.5)
      55         852 :   {}
      56             : 
      57             :   // Default destructor
      58             :   ~VariationalSmootherSystem() override;
      59             : 
      60             :   /**
      61             :    * Assembly method to update the mesh based on the smoother solve.
      62             :    */
      63             :   virtual void assembly (bool get_residual,
      64             :                          bool get_jacobian,
      65             :                          bool apply_heterogeneous_constraints = false,
      66             :                          bool apply_no_constraints = false) override;
      67             : 
      68          24 :   Real & get_dilation_weight() { return _dilation_weight; }
      69             : 
      70             :   /**
      71             :    * Get the target element for a given element type.
      72             :    * @param type Element type
      73             :    * @return a std::pair containing the target element for type and the
      74             :    * corresponding nodes that must be kept in scope while the target element is
      75             :    * used.
      76             :    */
      77             :   static std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
      78             :   get_target_elem(const ElemType & type);
      79             : 
      80             :   /**
      81             :    * Get the jacobians (and determinants) of the target-to-reference element mapping.
      82             :    * @param target_elem Target element.
      83             :    * @param femcontext Context used to build mapping.
      84             :    * @param jacobian Vector in which to store the jacobians for each quadrature point.
      85             :    * @param jacobian_dets Vector in which to store the determinant of the jacobians
      86             :    * for each quadrature point.
      87             :    */
      88             :   static void get_target_to_reference_jacobian(const Elem * const target_elem,
      89             :                                                const FEMContext & femcontext,
      90             :                                                std::vector<RealTensor> & jacobians,
      91             :                                                std::vector<Real> & jacobian_dets);
      92             : 
      93             : protected:
      94             : 
      95             :   // System initialization
      96             :   virtual void init_data () override;
      97             : 
      98             :   // Context initialization
      99             :   virtual void init_context (libMesh::DiffContext & context) override;
     100             : 
     101             :   // Element residual and jacobian calculations
     102             :   // Time dependent parts
     103             :   virtual bool element_time_derivative (bool request_jacobian,
     104             :                                         libMesh::DiffContext & context) override;
     105             : 
     106             :   /* Computes the element reference volume used in the dilation metric
     107             :    * The reference value is set to the averaged value of all elements' average
     108             :    * |J|. Also computes any applicable target element inverse Jacobians. Target
     109             :    * elements are relavant when the reference element does not minimize the
     110             :    * distortion metric.
     111             :    */
     112             :   void prepare_for_smoothing();
     113             : 
     114             :   /// The small nonzero constant to prevent zero denominators (degenerate elements only)
     115             :   const Real _epsilon_squared;
     116             : 
     117             :   /// The reference volume for each element
     118             :   Real _ref_vol;
     119             : 
     120             :   /// The relative weight to give the dilation metric. The distortion metric is given weight 1 - _dilation_weight.
     121             :   Real _dilation_weight;
     122             : 
     123             :   /* Map to hold target qp-dependent element target-to-reference mapping
     124             :    * Jacobians, if any
     125             :    */
     126             :   std::map<ElemType, std::vector<RealTensor>> _target_jacobians;
     127             : };
     128             : 
     129             : } // namespace libMesh
     130             : 
     131             : #endif // LIBMESH_VARIATIONAL_SMOOTHER_SYSTEM_H