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 :
19 :
20 : #ifndef LIBMESH_FEM_SYSTEM_H
21 : #define LIBMESH_FEM_SYSTEM_H
22 :
23 : // Local Includes
24 : #include "libmesh/diff_system.h"
25 : #include "libmesh/fem_physics.h"
26 :
27 : // C++ includes
28 : #include <cstddef>
29 :
30 : namespace libMesh
31 : {
32 :
33 : // Forward Declarations
34 : class DiffContext;
35 : class FEMContext;
36 :
37 :
38 : /**
39 : * This class provides a specific system class. It aims
40 : * at nonlinear implicit systems, requiring only a
41 : * cell residual calculation from the user.
42 : *
43 : * \note Additional vectors/matrices can be added via parent class
44 : * interfaces.
45 : *
46 : * This class is part of the new DifferentiableSystem framework,
47 : * which is still experimental. Users of this framework should
48 : * beware of bugs and future API changes.
49 : *
50 : * \author Roy H. Stogner
51 : * \date 2006
52 : */
53 300 : class FEMSystem : public DifferentiableSystem,
54 : public FEMPhysics
55 : {
56 : public:
57 :
58 : /**
59 : * Constructor. Optionally initializes required
60 : * data structures.
61 : */
62 : FEMSystem (EquationSystems & es,
63 : const std::string & name,
64 : const unsigned int number);
65 :
66 : /**
67 : * Special functions.
68 : * - This class has the same restrictions as its base class.
69 : * - The destructor is defaulted out-of-line.
70 : */
71 : FEMSystem (const FEMSystem &) = delete;
72 : FEMSystem & operator= (const FEMSystem &) = delete;
73 : FEMSystem (FEMSystem &&) = delete;
74 : FEMSystem & operator= (FEMSystem &&) = delete;
75 : virtual ~FEMSystem ();
76 :
77 : /**
78 : * The type of system.
79 : */
80 : typedef FEMSystem sys_type;
81 :
82 : /**
83 : * The type of the parent.
84 : */
85 : typedef DifferentiableSystem Parent;
86 :
87 : /**
88 : * Prepares \p matrix or \p rhs for matrix assembly.
89 : * Users may reimplement this to add pre- or post-assembly
90 : * code before or after calling FEMSystem::assembly().
91 : * Note that in some cases only
92 : * \link current_local_solution \endlink is used during assembly,
93 : * and, therefore, if \link solution \endlink has been altered
94 : * without \link update() \endlink being called, then the
95 : * user must call \link update() \endlink before calling
96 : * this function.
97 : */
98 : virtual void assembly (bool get_residual,
99 : bool get_jacobian,
100 : bool apply_heterogeneous_constraints = false,
101 : bool apply_no_constraints = false) override;
102 :
103 : /**
104 : * Invokes the solver associated with the system. For steady state
105 : * solvers, this will find a root x where F(x) = 0. For transient
106 : * solvers, this will integrate dx/dt = F(x).
107 : *
108 : * For moving mesh systems, this also translates the mesh to the
109 : * solution position.
110 : */
111 : virtual void solve () override;
112 :
113 : /**
114 : * Tells the FEMSystem to set the degree of freedom coefficients
115 : * which should correspond to mesh nodal coordinates.
116 : */
117 : void mesh_position_get();
118 :
119 : /**
120 : * Tells the FEMSystem to set the mesh nodal coordinates
121 : * which should correspond to degree of freedom coefficients.
122 : */
123 : void mesh_position_set();
124 :
125 : /**
126 : * Builds a FEMContext object with enough information to do
127 : * evaluations on each element.
128 : *
129 : * For most problems, the default FEMSystem implementation is correct; users
130 : * who subclass FEMContext will need to also reimplement this method to build
131 : * it.
132 : */
133 : virtual std::unique_ptr<DiffContext> build_context() override;
134 :
135 : /*
136 : * Prepares the result of a build_context() call for use.
137 : *
138 : * Most FEMSystem-based problems will need to reimplement this in order to
139 : * call FE::get_*() as their particular physics requires.
140 : */
141 : virtual void init_context(DiffContext &) override;
142 :
143 : /**
144 : * Runs a postprocessing loop over all elements, and if
145 : * \p postprocess_sides is true over all sides.
146 : */
147 : virtual void postprocess () override;
148 :
149 : /**
150 : * Runs a qoi assembly loop over all elements, and if
151 : * \p assemble_qoi_sides is true over all sides.
152 : *
153 : * Users may have to override this function if they have any
154 : * quantities of interest that are not expressible as a sum of
155 : * element qois.
156 : */
157 : virtual void assemble_qoi (const QoISet & indices = QoISet()) override;
158 :
159 : /**
160 : * Runs a qoi derivative assembly loop over all elements, and if
161 : * \p assemble_qoi_sides is true over all sides.
162 : *
163 : * Users may have to override this function for quantities of
164 : * interest that are not expressible as a sum of element qois.
165 : */
166 : virtual void assemble_qoi_derivative (const QoISet & qoi_indices = QoISet(),
167 : bool include_liftfunc = true,
168 : bool apply_constraints = true) override;
169 :
170 : /**
171 : * If fe_reinit_during_postprocess is true (it is true by default), FE
172 : * objects will be reinit()ed with their default quadrature rules. If false,
173 : * FE objects will need to be reinit()ed by the user or will be in an
174 : * undefined state.
175 : */
176 : bool fe_reinit_during_postprocess;
177 :
178 : /**
179 : * If calculating numeric jacobians is required, the FEMSystem
180 : * will perturb each solution vector entry by numerical_jacobian_h
181 : * when calculating finite differences. This defaults to the
182 : * libMesh TOLERANCE but can be set manually.
183 : *
184 : * For ALE terms, the FEMSystem will perturb each mesh point in an
185 : * element by numerical_jacobian_h * Elem::hmin()
186 : */
187 : Real numerical_jacobian_h;
188 :
189 : /**
190 : * If numerical_jacobian_h_for_var(var_num) is changed from its
191 : * default value (numerical_jacobian_h), the FEMSystem will perturb
192 : * solution vector entries for variable var_num by that amount when
193 : * calculating finite differences with respect to that variable.
194 : *
195 : * This is useful in multiphysics problems which have not been
196 : * normalized.
197 : */
198 : Real numerical_jacobian_h_for_var(unsigned int var_num) const;
199 :
200 : void set_numerical_jacobian_h_for_var(unsigned int var_num, Real new_h);
201 :
202 : /**
203 : * If verify_analytic_jacobian is equal to zero (as it is by
204 : * default), no numeric jacobians will be calculated unless
205 : * an overridden element_time_derivative(), element_constraint(),
206 : * side_time_derivative(), or side_constraint() function cannot
207 : * provide an analytic jacobian upon request.
208 : *
209 : * If verify_analytic_jacobian is equal to the positive value tol,
210 : * then any time a full analytic element jacobian can be calculated
211 : * it will be tested against a numerical jacobian on the same element,
212 : * and the program will abort if the relative error (in matrix l1 norms)
213 : * exceeds tol.
214 : */
215 : Real verify_analytic_jacobians;
216 :
217 : /**
218 : * Syntax sugar to make numerical_jacobian() declaration easier.
219 : */
220 : typedef bool (TimeSolver::*TimeSolverResPtr)(bool, DiffContext &);
221 :
222 : /**
223 : * Uses the results of multiple \p res calls
224 : * to numerically differentiate the corresponding jacobian.
225 : */
226 : void numerical_jacobian (TimeSolverResPtr res,
227 : FEMContext & context) const;
228 :
229 : /**
230 : * Uses the results of multiple element_residual() calls
231 : * to numerically differentiate the corresponding jacobian
232 : * on an element.
233 : */
234 : void numerical_elem_jacobian (FEMContext & context) const;
235 :
236 : /**
237 : * Uses the results of multiple side_residual() calls
238 : * to numerically differentiate the corresponding jacobian
239 : * on an element's side.
240 : */
241 : void numerical_side_jacobian (FEMContext & context) const;
242 :
243 : /**
244 : * Uses the results of multiple side_residual() calls
245 : * to numerically differentiate the corresponding jacobian
246 : * on nonlocal DoFs.
247 : */
248 : void numerical_nonlocal_jacobian (FEMContext & context) const;
249 :
250 : protected:
251 : /**
252 : * Initializes the member data fields associated with
253 : * the system, so that, e.g., \p assemble() may be used.
254 : */
255 : virtual void init_data () override;
256 :
257 : private:
258 : std::vector<Real> _numerical_jacobian_h_for_var;
259 : };
260 :
261 : // --------------------------------------------------------------
262 : // FEMSystem inline methods
263 : inline
264 : Real
265 298738 : FEMSystem::numerical_jacobian_h_for_var(unsigned int var_num) const
266 : {
267 362082 : if ((var_num >= _numerical_jacobian_h_for_var.size()) ||
268 0 : _numerical_jacobian_h_for_var[var_num] == Real(0))
269 330410 : return numerical_jacobian_h;
270 :
271 0 : return _numerical_jacobian_h_for_var[var_num];
272 : }
273 :
274 : inline
275 : void FEMSystem::set_numerical_jacobian_h_for_var(unsigned int var_num,
276 : Real new_h)
277 : {
278 : if (_numerical_jacobian_h_for_var.size() <= var_num)
279 : _numerical_jacobian_h_for_var.resize(var_num+1,Real(0));
280 :
281 : libmesh_assert_greater(new_h, Real(0));
282 :
283 : _numerical_jacobian_h_for_var[var_num] = new_h;
284 : }
285 :
286 : } // namespace libMesh
287 :
288 :
289 : #endif // LIBMESH_FEM_SYSTEM_H
|
