1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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/quadrature.h"
30 #include "libmesh/sparse_matrix.h"
31 #include "libmesh/time_solver.h"
32 #include "libmesh/unsteady_solver.h" // For eulerian_residual
33
34
35 namespace libMesh
36 {
37
eulerian_residual(bool request_jacobian,DiffContext &)38 bool FEMPhysics::eulerian_residual (bool request_jacobian,
39 DiffContext & /*c*/)
40 {
41 // Only calculate a mesh movement residual if it's necessary
42 if (!_mesh_sys)
43 return request_jacobian;
44
45 libmesh_not_implemented();
46
47 #if 0
48 FEMContext & context = cast_ref<FEMContext &>(c);
49
50 // This function only supports fully coupled mesh motion for now
51 libmesh_assert_equal_to (_mesh_sys, this);
52
53 unsigned int n_qpoints = (context.get_element_qrule())->n_points();
54
55 const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
56 0 : context.dof_indices_var[_mesh_x_var].size();
57 const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
58 0 : context.dof_indices_var[_mesh_y_var].size();
59 const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
60 0 : context.dof_indices_var[_mesh_z_var].size();
61
62 const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
63 (n_y_dofs ? _mesh_y_var :
64 (n_z_dofs ? _mesh_z_var :
65 libMesh::invalid_uint));
66
67 // If we're our own _mesh_sys, we'd better be in charge of
68 // at least one coordinate, and we'd better have the same
69 // FE type for all coordinates we are in charge of
70 libmesh_assert_not_equal_to (mesh_xyz_var, libMesh::invalid_uint);
71 libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
72 context.element_fe_var[mesh_xyz_var]);
73 libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
74 context.element_fe_var[mesh_xyz_var]);
75 libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
76 context.element_fe_var[mesh_xyz_var]);
77
78 const std::vector<std::vector<Real>> & psi =
79 context.element_fe_var[mesh_xyz_var]->get_phi();
80
81 for (auto var : make_range(context.n_vars()))
82 {
83 // Mesh motion only affects time-evolving variables
84 if (this->is_time_evolving(var))
85 continue;
86
87 // The mesh coordinate variables themselves are Lagrangian,
88 // not Eulerian, and no convective term is desired.
89 if (/*_mesh_sys == this && */
90 (var == _mesh_x_var ||
91 var == _mesh_y_var ||
92 var == _mesh_z_var))
93 continue;
94
95 // Some of this code currently relies on the assumption that
96 // we can pull mesh coordinate data from our own system
97 if (_mesh_sys != this)
98 libmesh_not_implemented();
99
100 // This residual should only be called by unsteady solvers:
101 // if the mesh is steady, there's no mesh convection term!
102 UnsteadySolver * unsteady;
103 if (this->time_solver->is_steady())
104 return request_jacobian;
105 else
106 unsteady = cast_ptr<UnsteadySolver*>(this->time_solver.get());
107
108 const std::vector<Real> & JxW =
109 context.element_fe_var[var]->get_JxW();
110
111 const std::vector<std::vector<Real>> & phi =
112 context.element_fe_var[var]->get_phi();
113
114 const std::vector<std::vector<RealGradient>> & dphi =
115 context.element_fe_var[var]->get_dphi();
116
117 const unsigned int n_u_dofs = context.dof_indices_var[var].size();
118
119 DenseSubVector<Number> & Fu = *context.elem_subresiduals[var];
120 DenseSubMatrix<Number> & Kuu = *context.elem_subjacobians[var][var];
121
122 DenseSubMatrix<Number> * Kux = n_x_dofs ?
123 context.elem_subjacobians[var][_mesh_x_var] : nullptr;
124 DenseSubMatrix<Number> * Kuy = n_y_dofs ?
125 context.elem_subjacobians[var][_mesh_y_var] : nullptr;
126 DenseSubMatrix<Number> * Kuz = n_z_dofs ?
127 context.elem_subjacobians[var][_mesh_z_var] : nullptr;
128
129 std::vector<Real> delta_x(n_x_dofs, 0.);
130 std::vector<Real> delta_y(n_y_dofs, 0.);
131 std::vector<Real> delta_z(n_z_dofs, 0.);
132
133 for (unsigned int i = 0; i != n_x_dofs; ++i)
134 {
135 unsigned int j = context.dof_indices_var[_mesh_x_var][i];
136 delta_x[i] = libmesh_real(this->current_solution(j)) -
137 libmesh_real(unsteady->old_nonlinear_solution(j));
138 }
139
140 for (unsigned int i = 0; i != n_y_dofs; ++i)
141 {
142 unsigned int j = context.dof_indices_var[_mesh_y_var][i];
143 delta_y[i] = libmesh_real(this->current_solution(j)) -
144 libmesh_real(unsteady->old_nonlinear_solution(j));
145 }
146
147 for (unsigned int i = 0; i != n_z_dofs; ++i)
148 {
149 unsigned int j = context.dof_indices_var[_mesh_z_var][i];
150 delta_z[i] = libmesh_real(this->current_solution(j)) -
151 libmesh_real(unsteady->old_nonlinear_solution(j));
152 }
153
154 for (unsigned int qp = 0; qp != n_qpoints; ++qp)
155 {
156 Gradient grad_u = context.interior_gradient(var, qp);
157 RealGradient convection(0.);
158
159 for (unsigned int i = 0; i != n_x_dofs; ++i)
160 convection(0) += delta_x[i] * psi[i][qp];
161 for (unsigned int i = 0; i != n_y_dofs; ++i)
162 convection(1) += delta_y[i] * psi[i][qp];
163 for (unsigned int i = 0; i != n_z_dofs; ++i)
164 convection(2) += delta_z[i] * psi[i][qp];
165
166 for (unsigned int i = 0; i != n_u_dofs; ++i)
167 {
168 Number JxWxPhiI = JxW[qp] * phi[i][qp];
169 Fu(i) += (convection * grad_u) * JxWxPhiI;
170 if (request_jacobian)
171 {
172 Number JxWxPhiI = JxW[qp] * phi[i][qp];
173 for (unsigned int j = 0; j != n_u_dofs; ++j)
174 Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);
175
176 Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;
177
178 Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
179 for (unsigned int j = 0; j != n_x_dofs; ++j)
180 (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];
181
182 Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
183 for (unsigned int j = 0; j != n_y_dofs; ++j)
184 (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];
185
186 Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
187 for (unsigned int j = 0; j != n_z_dofs; ++j)
188 (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
189 }
190 }
191 }
192 }
193 #endif // 0
194
195 return request_jacobian;
196 }
197
198
199
mass_residual(bool request_jacobian,DiffContext & c)200 bool FEMPhysics::mass_residual (bool request_jacobian,
201 DiffContext & c)
202 {
203 FEMContext & context = cast_ref<FEMContext &>(c);
204
205 unsigned int n_qpoints = context.get_element_qrule().n_points();
206
207 for (auto var : make_range(context.n_vars()))
208 {
209 if (!this->is_time_evolving(var))
210 continue;
211
212 FEBase * elem_fe = nullptr;
213 context.get_element_fe( var, elem_fe );
214
215 const std::vector<Real> & JxW = elem_fe->get_JxW();
216
217 const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
218
219 const unsigned int n_dofs = cast_int<unsigned int>
220 (context.get_dof_indices(var).size());
221
222 DenseSubVector<Number> & Fu = context.get_elem_residual(var);
223 DenseSubMatrix<Number> & Kuu = context.get_elem_jacobian( var, var );
224
225 for (unsigned int qp = 0; qp != n_qpoints; ++qp)
226 {
227 Number uprime;
228 context.interior_rate(var, qp, uprime);
229 const Number JxWxU = JxW[qp] * uprime;
230 for (unsigned int i = 0; i != n_dofs; ++i)
231 {
232 Fu(i) -= JxWxU * phi[i][qp];
233 if (request_jacobian && context.elem_solution_rate_derivative)
234 {
235 const Number JxWxPhiIxDeriv = JxW[qp] * phi[i][qp] *
236 context.elem_solution_rate_derivative;
237 Kuu(i,i) -= JxWxPhiIxDeriv * phi[i][qp];
238 for (unsigned int j = i+1; j < n_dofs; ++j)
239 {
240 const Number Kij = JxWxPhiIxDeriv * phi[j][qp];
241 Kuu(i,j) -= Kij;
242 Kuu(j,i) -= Kij;
243 }
244 }
245 }
246 }
247 }
248
249 return request_jacobian;
250 }
251
252 } // namespace libMesh
253