libMesh
fem_physics.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 
19 
20 #include "libmesh/dof_map.h"
21 #include "libmesh/elem.h"
22 #include "libmesh/equation_systems.h"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/fem_context.h"
25 #include "libmesh/fem_system.h"
26 #include "libmesh/libmesh_logging.h"
27 #include "libmesh/mesh_base.h"
28 #include "libmesh/numeric_vector.h"
29 #include "libmesh/parallel.h"
30 #include "libmesh/quadrature.h"
31 #include "libmesh/sparse_matrix.h"
32 #include "libmesh/time_solver.h"
33 #include "libmesh/unsteady_solver.h" // For eulerian_residual
34 
35 
36 namespace libMesh
37 {
38 
39 bool FEMPhysics::eulerian_residual (bool request_jacobian,
40  DiffContext & /*c*/)
41 {
42  // Only calculate a mesh movement residual if it's necessary
43  if (!_mesh_sys)
44  return request_jacobian;
45 
46  libmesh_not_implemented();
47 
48 #if 0
49  FEMContext & context = cast_ref<FEMContext &>(c);
50 
51  // This function only supports fully coupled mesh motion for now
52  libmesh_assert_equal_to (_mesh_sys, this);
53 
54  unsigned int n_qpoints = (context.get_element_qrule())->n_points();
55 
56  const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
57  0 : context.dof_indices_var[_mesh_x_var].size();
58  const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
59  0 : context.dof_indices_var[_mesh_y_var].size();
60  const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
61  0 : context.dof_indices_var[_mesh_z_var].size();
62 
63  const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
64  (n_y_dofs ? _mesh_y_var :
65  (n_z_dofs ? _mesh_z_var :
67 
68  // If we're our own _mesh_sys, we'd better be in charge of
69  // at least one coordinate, and we'd better have the same
70  // FE type for all coordinates we are in charge of
71  libmesh_assert_not_equal_to (mesh_xyz_var, libMesh::invalid_uint);
72  libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
73  context.element_fe_var[mesh_xyz_var]);
74  libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
75  context.element_fe_var[mesh_xyz_var]);
76  libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
77  context.element_fe_var[mesh_xyz_var]);
78 
79  const std::vector<std::vector<Real>> & psi =
80  context.element_fe_var[mesh_xyz_var]->get_phi();
81 
82  for (unsigned int var = 0; var != context.n_vars(); ++var)
83  {
84  // Mesh motion only affects time-evolving variables
85  if (this->is_time_evolving(var))
86  continue;
87 
88  // The mesh coordinate variables themselves are Lagrangian,
89  // not Eulerian, and no convective term is desired.
90  if (/*_mesh_sys == this && */
91  (var == _mesh_x_var ||
92  var == _mesh_y_var ||
93  var == _mesh_z_var))
94  continue;
95 
96  // Some of this code currently relies on the assumption that
97  // we can pull mesh coordinate data from our own system
98  if (_mesh_sys != this)
99  libmesh_not_implemented();
100 
101  // This residual should only be called by unsteady solvers:
102  // if the mesh is steady, there's no mesh convection term!
103  UnsteadySolver * unsteady;
104  if (this->time_solver->is_steady())
105  return request_jacobian;
106  else
107  unsteady = cast_ptr<UnsteadySolver*>(this->time_solver.get());
108 
109  const std::vector<Real> & JxW =
110  context.element_fe_var[var]->get_JxW();
111 
112  const std::vector<std::vector<Real>> & phi =
113  context.element_fe_var[var]->get_phi();
114 
115  const std::vector<std::vector<RealGradient>> & dphi =
116  context.element_fe_var[var]->get_dphi();
117 
118  const unsigned int n_u_dofs = context.dof_indices_var[var].size();
119 
120  DenseSubVector<Number> & Fu = *context.elem_subresiduals[var];
121  DenseSubMatrix<Number> & Kuu = *context.elem_subjacobians[var][var];
122 
123  DenseSubMatrix<Number> * Kux = n_x_dofs ?
124  context.elem_subjacobians[var][_mesh_x_var] : libmesh_nullptr;
125  DenseSubMatrix<Number> * Kuy = n_y_dofs ?
126  context.elem_subjacobians[var][_mesh_y_var] : libmesh_nullptr;
127  DenseSubMatrix<Number> * Kuz = n_z_dofs ?
128  context.elem_subjacobians[var][_mesh_z_var] : libmesh_nullptr;
129 
130  std::vector<Real> delta_x(n_x_dofs, 0.);
131  std::vector<Real> delta_y(n_y_dofs, 0.);
132  std::vector<Real> delta_z(n_z_dofs, 0.);
133 
134  for (unsigned int i = 0; i != n_x_dofs; ++i)
135  {
136  unsigned int j = context.dof_indices_var[_mesh_x_var][i];
137  delta_x[i] = libmesh_real(this->current_solution(j)) -
138  libmesh_real(unsteady->old_nonlinear_solution(j));
139  }
140 
141  for (unsigned int i = 0; i != n_y_dofs; ++i)
142  {
143  unsigned int j = context.dof_indices_var[_mesh_y_var][i];
144  delta_y[i] = libmesh_real(this->current_solution(j)) -
145  libmesh_real(unsteady->old_nonlinear_solution(j));
146  }
147 
148  for (unsigned int i = 0; i != n_z_dofs; ++i)
149  {
150  unsigned int j = context.dof_indices_var[_mesh_z_var][i];
151  delta_z[i] = libmesh_real(this->current_solution(j)) -
152  libmesh_real(unsteady->old_nonlinear_solution(j));
153  }
154 
155  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
156  {
157  Gradient grad_u = context.interior_gradient(var, qp);
158  RealGradient convection(0.);
159 
160  for (unsigned int i = 0; i != n_x_dofs; ++i)
161  convection(0) += delta_x[i] * psi[i][qp];
162  for (unsigned int i = 0; i != n_y_dofs; ++i)
163  convection(1) += delta_y[i] * psi[i][qp];
164  for (unsigned int i = 0; i != n_z_dofs; ++i)
165  convection(2) += delta_z[i] * psi[i][qp];
166 
167  for (unsigned int i = 0; i != n_u_dofs; ++i)
168  {
169  Number JxWxPhiI = JxW[qp] * phi[i][qp];
170  Fu(i) += (convection * grad_u) * JxWxPhiI;
171  if (request_jacobian)
172  {
173  Number JxWxPhiI = JxW[qp] * phi[i][qp];
174  for (unsigned int j = 0; j != n_u_dofs; ++j)
175  Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);
176 
177  Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;
178 
179  Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
180  for (unsigned int j = 0; j != n_x_dofs; ++j)
181  (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];
182 
183  Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
184  for (unsigned int j = 0; j != n_y_dofs; ++j)
185  (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];
186 
187  Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
188  for (unsigned int j = 0; j != n_z_dofs; ++j)
189  (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
190  }
191  }
192  }
193  }
194 #endif // 0
195 
196  return request_jacobian;
197 }
198 
199 
200 
201 bool FEMPhysics::mass_residual (bool request_jacobian,
202  DiffContext & c)
203 {
204  FEMContext & context = cast_ref<FEMContext &>(c);
205 
206  unsigned int n_qpoints = context.get_element_qrule().n_points();
207 
208  for (unsigned int var = 0; var != context.n_vars(); ++var)
209  {
210  if (!this->is_time_evolving(var))
211  continue;
212 
213  FEBase * elem_fe = libmesh_nullptr;
214  context.get_element_fe( var, elem_fe );
215 
216  const std::vector<Real> & JxW = elem_fe->get_JxW();
217 
218  const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
219 
220  const unsigned int n_dofs = cast_int<unsigned int>
221  (context.get_dof_indices(var).size());
222 
223  DenseSubVector<Number> & Fu = context.get_elem_residual(var);
224  DenseSubMatrix<Number> & Kuu = context.get_elem_jacobian( var, var );
225 
226  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
227  {
228  Number uprime;
229  context.interior_rate(var, qp, uprime);
230  const Number JxWxU = JxW[qp] * uprime;
231  for (unsigned int i = 0; i != n_dofs; ++i)
232  {
233  Fu(i) -= JxWxU * phi[i][qp];
234  if (request_jacobian && context.elem_solution_rate_derivative)
235  {
236  const Number JxWxPhiIxDeriv = JxW[qp] * phi[i][qp] *
238  Kuu(i,i) -= JxWxPhiIxDeriv * phi[i][qp];
239  for (unsigned int j = i+1; j < n_dofs; ++j)
240  {
241  const Number Kij = JxWxPhiIxDeriv * phi[j][qp];
242  Kuu(i,j) -= Kij;
243  Kuu(j,i) -= Kij;
244  }
245  }
246  }
247  }
248  }
249 
250  return request_jacobian;
251 }
252 
253 } // namespace libMesh
T libmesh_real(T a)
unsigned int _mesh_x_var
Variables from which to acquire moving mesh information.
Definition: diff_physics.h:546
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:184
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
bool is_time_evolving(unsigned int var) const
Definition: diff_physics.h:276
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
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.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
Number old_nonlinear_solution(const dof_id_type global_dof_number) const
Real elem_solution_rate_derivative
The derivative of elem_solution_rate with respect to the current nonlinear solution, for use by systems with non default mass_residual terms.
Definition: diff_context.h:490
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
virtual bool mass_residual(bool request_jacobian, DiffContext &) libmesh_override
Subtracts a mass vector contribution on elem from elem_residual.
Definition: fem_physics.C:201
Defines a dense submatrix for use in Finite Element-type computations.
System * _mesh_sys
System from which to acquire moving mesh information.
Definition: diff_physics.h:541
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
virtual bool eulerian_residual(bool request_jacobian, DiffContext &context) libmesh_override
Adds a pseudo-convection contribution on elem to elem_residual, if the nodes of elem are being transl...
Definition: fem_physics.C:39
This is a generic class that defines a solver to handle time integration of DifferentiableSystems.
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
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