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