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             : #include "libmesh/first_order_unsteady_solver.h"
      19             : #include "libmesh/diff_system.h"
      20             : #include "libmesh/quadrature.h"
      21             : 
      22             : namespace libMesh
      23             : {
      24             : 
      25      653824 : void FirstOrderUnsteadySolver::prepare_accel(DiffContext & context)
      26             : {
      27       60104 :   context.get_elem_solution_accel() = context.get_elem_solution_rate();
      28             : 
      29      653824 :   context.elem_solution_accel_derivative = context.get_elem_solution_rate_derivative();
      30      653824 : }
      31             : 
      32      980736 : bool FirstOrderUnsteadySolver::compute_second_order_eqns(bool compute_jacobian, DiffContext & c)
      33             : {
      34       90156 :   FEMContext & context = cast_ref<FEMContext &>(c);
      35             : 
      36       90156 :   unsigned int n_qpoints = context.get_element_qrule().n_points();
      37             : 
      38     6180192 :   for (auto var : make_range(context.n_vars()))
      39             :     {
      40     5670552 :       if (!this->_system.is_second_order_var(var))
      41     2364180 :         continue;
      42             : 
      43     2599728 :       unsigned int dot_var = this->_system.get_second_order_dot_var(var);
      44             : 
      45             :       // We're assuming that the FE space for var and dot_var are the same
      46      235548 :       libmesh_assert( context.get_system().variable(var).type() ==
      47             :                       context.get_system().variable(dot_var).type() );
      48             : 
      49     2364180 :       FEBase * elem_fe = nullptr;
      50     2364180 :       context.get_element_fe( var, elem_fe );
      51             : 
      52     2599728 :       const std::vector<Real> & JxW = elem_fe->get_JxW();
      53             : 
      54      235548 :       const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
      55             : 
      56             :       const unsigned int n_dofs = cast_int<unsigned int>
      57      471096 :         (context.get_dof_indices(dot_var).size());
      58             : 
      59      235548 :       DenseSubVector<Number> & Fu = context.get_elem_residual(var);
      60      235548 :       DenseSubMatrix<Number> & Kuu = context.get_elem_jacobian( var, var );
      61      235548 :       DenseSubMatrix<Number> & Kuv = context.get_elem_jacobian( var, dot_var );
      62             : 
      63    21026400 :       for (unsigned int qp = 0; qp != n_qpoints; ++qp)
      64             :         {
      65      895020 :           Number udot, v;
      66    18426672 :           context.interior_rate(var, qp, udot);
      67    18426672 :           context.interior_value(dot_var, qp, v);
      68             : 
      69   165877344 :           for (unsigned int i = 0; i < n_dofs; i++)
      70             :             {
      71   186882948 :               Fu(i) += JxW[qp]*(udot-v)*phi[i][qp];
      72             : 
      73   147450672 :               if (compute_jacobian)
      74             :                 {
      75    73725336 :                   Number rate_factor = JxW[qp]*context.get_elem_solution_rate_derivative()*phi[i][qp];
      76    73725336 :                   Number soln_factor = JxW[qp]*context.get_elem_solution_derivative()*phi[i][qp];
      77             : 
      78    70144386 :                   Kuu(i,i) += rate_factor*phi[i][qp];
      79    70144386 :                   Kuv(i,i) -= soln_factor*phi[i][qp];
      80             : 
      81   337192344 :                   for (unsigned int j = i+1; j < n_dofs; j++)
      82             :                     {
      83   310579200 :                       Kuu(i,j) += rate_factor*phi[j][qp];
      84   250656768 :                       Kuu(j,i) += rate_factor*phi[j][qp];
      85             : 
      86   250656768 :                       Kuv(i,j) -= soln_factor*phi[j][qp];
      87   274212864 :                       Kuv(j,i) -= soln_factor*phi[j][qp];
      88             :                     }
      89             :                 }
      90             :             }
      91             :         }
      92             :     }
      93             : 
      94      980736 :   return compute_jacobian;
      95             : }
      96             : 
      97             : } // namespace libMesh