libMesh
first_order_unsteady_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2016 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 
26 {
27  context.get_elem_solution_accel() = context.get_elem_solution_rate();
28 
30 }
31 
33 {
34  FEMContext & context = cast_ref<FEMContext &>(c);
35 
36  unsigned int n_qpoints = context.get_element_qrule().n_points();
37 
38  for (unsigned int var = 0; var != context.n_vars(); ++var)
39  {
40  if (!this->_system.is_second_order_var(var))
41  continue;
42 
43  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  libmesh_assert( context.get_system().variable(var).type() ==
47  context.get_system().variable(dot_var).type() );
48 
49  FEBase * elem_fe = libmesh_nullptr;
50  context.get_element_fe( var, elem_fe );
51 
52  const std::vector<Real> & JxW = elem_fe->get_JxW();
53 
54  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
55 
56  const unsigned int n_dofs = cast_int<unsigned int>
57  (context.get_dof_indices(dot_var).size());
58 
59  DenseSubVector<Number> & Fu = context.get_elem_residual(var);
60  DenseSubMatrix<Number> & Kuu = context.get_elem_jacobian( var, var );
61  DenseSubMatrix<Number> & Kuv = context.get_elem_jacobian( var, dot_var );
62 
63  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
64  {
65  Number udot, v;
66  context.interior_rate(var, qp, udot);
67  context.interior_value(dot_var, qp, v);
68 
69  for (unsigned int i = 0; i < n_dofs; i++)
70  {
71  Fu(i) += JxW[qp]*(udot-v)*phi[i][qp];
72 
73  if (compute_jacobian)
74  {
75  Number rate_factor = JxW[qp]*context.get_elem_solution_rate_derivative()*phi[i][qp];
76  Number soln_factor = JxW[qp]*context.get_elem_solution_derivative()*phi[i][qp];
77 
78  Kuu(i,i) += rate_factor*phi[i][qp];
79  Kuv(i,i) -= soln_factor*phi[i][qp];
80 
81  for (unsigned int j = i+1; j < n_dofs; j++)
82  {
83  Kuu(i,j) += rate_factor*phi[j][qp];
84  Kuu(j,i) += rate_factor*phi[j][qp];
85 
86  Kuv(i,j) -= soln_factor*phi[j][qp];
87  Kuv(j,i) -= soln_factor*phi[j][qp];
88  }
89  }
90  }
91  }
92  }
93 
94  return compute_jacobian;
95 }
96 
97 } // namespace libMesh
const FEType & type() const
Definition: variable.h:119
bool compute_second_order_eqns(bool compute_jacobian, DiffContext &c)
If there are second order variables, then we need to compute their residual equations and correspondi...
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
void prepare_accel(DiffContext &context)
If there are second order variables in the system, then we also prepare the accel for those variables...
const DenseVector< Number > & get_elem_solution_accel() const
Accessor for element solution accel of change w.r.t.
Definition: diff_context.h:179
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const class libmesh_nullptr_t libmesh_nullptr
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
The libMesh namespace provides an interface to certain functionality in the library.
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1273
libmesh_assert(j)
Defines a dense subvector for use in finite element computations.
unsigned int get_second_order_dot_var(unsigned int var) const
For a given second order (in time) variable var, this method will return the index to the correspondi...
Definition: diff_system.C:313
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
sys_type & _system
A reference to the system we are solving.
Definition: time_solver.h:257
Real get_elem_solution_rate_derivative() const
The derivative of the current elem_solution_rate w.r.t.
Definition: diff_context.h:428
const Variable & variable(unsigned int var) const
Return a constant reference to Variable var.
Definition: system.h:2114
Defines a dense submatrix for use in Finite Element-type computations.
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
const System & get_system() const
Accessor for associated system.
Definition: diff_context.h:104
const DenseVector< Number > & get_elem_solution_rate() const
Accessor for element solution rate of change w.r.t.
Definition: diff_context.h:144
bool is_second_order_var(unsigned int var) const
Definition: diff_physics.h:532
Real get_elem_solution_derivative() const
The derivative of the current elem_solution w.r.t.
Definition: diff_context.h:419
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
Real elem_solution_accel_derivative
The derivative of elem_solution_accel with respect to the current nonlinear solution, for use by systems with non default mass_residual terms.
Definition: diff_context.h:497
unsigned int n_points() const
Definition: quadrature.h:113
This class forms the foundation from which generic finite elements may be derived.
unsigned int n_vars() const
Number of variables in solution.
Definition: diff_context.h:98