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 
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 = 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
bool is_second_order_var(unsigned int var) const
Definition: diff_physics.h:533
Real get_elem_solution_derivative() const
Definition: diff_context.h:436
const DenseMatrix< Number > & get_elem_jacobian() const
Definition: diff_context.h:283
bool compute_second_order_eqns(bool compute_jacobian, DiffContext &c)
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:306
const Variable & variable(unsigned int var) const
Definition: system.h:2133
const DenseVector< Number > & get_elem_solution_rate() const
Definition: diff_context.h:145
sys_type & _system
Definition: time_solver.h:258
unsigned int get_second_order_dot_var(unsigned int var) const
Definition: diff_system.C:306
const System & get_system() const
Definition: diff_context.h:105
const std::vector< dof_id_type > & get_dof_indices() const
Definition: diff_context.h:367
unsigned int n_points() const
Definition: quadrature.h:127
const DenseVector< Number > & get_elem_residual() const
Definition: diff_context.h:249
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1253
const DenseVector< Number > & get_elem_solution_accel() const
Definition: diff_context.h:180
Real get_elem_solution_rate_derivative() const
Definition: diff_context.h:445
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Definition: fem_context.h:262
unsigned int n_vars() const
Definition: diff_context.h:99
Real elem_solution_accel_derivative
Definition: diff_context.h:514
const QBase & get_element_qrule() const
Definition: fem_context.h:765
const FEType & type() const
Definition: variable.h:119