Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : #ifndef LIBMESH_STATIC_CONDENSATION_H
19 : #define LIBMESH_STATIC_CONDENSATION_H
20 :
21 : #include "libmesh/libmesh_config.h"
22 :
23 : // Forward declarations
24 : namespace libMesh {
25 : class Elem;
26 : }
27 :
28 : // shell matrices currently only work with petsc,
29 : // and we currently rely on Eigen for local LU factorizations
30 : #if defined(LIBMESH_HAVE_EIGEN) && defined(LIBMESH_HAVE_PETSC)
31 :
32 : #include "libmesh/petsc_matrix_shell_matrix.h"
33 : #include "libmesh/id_types.h"
34 : #include "libmesh/libmesh_common.h"
35 : #include "libmesh/dense_matrix.h"
36 : #include "libmesh/variable.h"
37 : #include "libmesh/sparsity_pattern.h"
38 :
39 : #include <unordered_map>
40 : #include <memory>
41 : #include <vector>
42 :
43 : // Warnings about stack protection
44 : #include "libmesh/ignore_warnings.h"
45 : #include <Eigen/Dense>
46 : #include "libmesh/restore_warnings.h"
47 :
48 : namespace libMesh
49 : {
50 : class MeshBase;
51 : class System;
52 : class DofMap;
53 : class Elem;
54 : template <typename>
55 : class LinearSolver;
56 : template <typename>
57 : class NumericVector;
58 : template <typename>
59 : class Preconditioner;
60 : class StaticCondensationPreconditioner;
61 : class StaticCondensationDofMap;
62 :
63 : typedef Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> EigenMatrix;
64 : typedef Eigen::Matrix<Number, Eigen::Dynamic, 1> EigenVector;
65 :
66 90 : class StaticCondensation : public PetscMatrixShellMatrix<Number>
67 : {
68 : public:
69 : StaticCondensation(const MeshBase & mesh,
70 : System & system,
71 : const DofMap & full_dof_map,
72 : StaticCondensationDofMap & reduced_dof_map);
73 : virtual ~StaticCondensation();
74 :
75 : //
76 : // SparseMatrix overrides
77 : //
78 :
79 : virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override;
80 :
81 : virtual SolverPackage solver_package() override;
82 :
83 : virtual void init(const numeric_index_type m,
84 : const numeric_index_type n,
85 : const numeric_index_type m_l,
86 : const numeric_index_type n_l,
87 : const numeric_index_type nnz = 30,
88 : const numeric_index_type noz = 10,
89 : const numeric_index_type blocksize = 1) override;
90 :
91 : virtual void init(const ParallelType type) override;
92 :
93 : virtual bool initialized() const override;
94 :
95 : virtual void clear() noexcept override;
96 :
97 : virtual void zero() override;
98 :
99 : virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override;
100 :
101 : virtual std::unique_ptr<SparseMatrix<Number>> clone() const override;
102 :
103 : virtual void close() override;
104 :
105 : virtual numeric_index_type m() const override;
106 :
107 0 : virtual numeric_index_type n() const override { return this->m(); }
108 :
109 : virtual numeric_index_type row_start() const override;
110 :
111 : virtual numeric_index_type row_stop() const override;
112 :
113 0 : virtual numeric_index_type col_start() const override { return this->row_start(); }
114 :
115 0 : virtual numeric_index_type col_stop() const override { return this->row_stop(); }
116 :
117 : virtual void
118 : set(const numeric_index_type i, const numeric_index_type j, const Number value) override;
119 :
120 : virtual void
121 : add(const numeric_index_type i, const numeric_index_type j, const Number value) override;
122 :
123 : virtual void add_matrix(const DenseMatrix<Number> & dm,
124 : const std::vector<numeric_index_type> & rows,
125 : const std::vector<numeric_index_type> & cols) override;
126 :
127 : virtual void add_matrix(const DenseMatrix<Number> & dm,
128 : const std::vector<numeric_index_type> & dof_indices) override;
129 :
130 : virtual void add(const Number a, const SparseMatrix<Number> & X) override;
131 :
132 : virtual Number operator()(const numeric_index_type i, const numeric_index_type j) const override;
133 :
134 : virtual Real l1_norm() const override;
135 :
136 : virtual Real linfty_norm() const override;
137 :
138 : virtual bool closed() const override;
139 :
140 : virtual void print_personal(std::ostream & os = libMesh::out) const override;
141 :
142 : virtual void get_diagonal(NumericVector<Number> & dest) const override;
143 :
144 : virtual void get_transpose(SparseMatrix<Number> & dest) const override;
145 :
146 : virtual void get_row(numeric_index_type i,
147 : std::vector<numeric_index_type> & indices,
148 : std::vector<Number> & values) const override;
149 :
150 : //
151 : // Unique methods
152 : //
153 :
154 : /**
155 : * Size the element matrices
156 : */
157 : void init();
158 :
159 : /**
160 : * A no-op to be consistent with shimming from the \p StaticCondenstionPreconditioner. "setup"
161 : * should take place at the end of matrix assembly, e.g. during a call to \p close()
162 : */
163 : void setup();
164 :
165 : /**
166 : * Perform our three stages, \p forward_elimination(), a \p solve() on the condensed system, and
167 : * finally a \p backwards_substitution()
168 : */
169 : void apply(const NumericVector<Number> & full_rhs, NumericVector<Number> & full_sol);
170 :
171 : const SparseMatrix<Number> & get_condensed_mat() const;
172 :
173 : /**
174 : * Set the current element. This enables fast lookups of local indices from global indices
175 : */
176 : void set_current_elem(const Elem & elem);
177 :
178 : /**
179 : * Get the preconditioning wrapper
180 : */
181 18 : StaticCondensationPreconditioner & get_preconditioner() { return *_scp; }
182 :
183 : /**
184 : * @returns The reduced system linear solver
185 : */
186 : LinearSolver<Number> & reduced_system_solver();
187 :
188 490 : virtual bool require_sparsity_pattern() const override { return false; }
189 :
190 : /**
191 : * Sets whether this matrix represents uncondensed dofs only. In that case when building the Schur
192 : * complement we won't attempt to invert zero element matrices corresponding to the condensed dofs
193 : */
194 : void uncondensed_dofs_only() { _uncondensed_dofs_only = true; }
195 :
196 : /**
197 : * Add \p vars to the list of variables not to condense. This can be useful when some variable's
198 : * equation is discretized with a DG method or if including the variable in the condensed block
199 : * diagonal would result in it being singular
200 : */
201 : void dont_condense_vars(const std::unordered_set<unsigned int> & vars);
202 :
203 : private:
204 : /**
205 : * Retrieves the degree of freedom values from \p global_vector corresponding to \p
206 : * elem_dof_indices, filling both \p elem_dof_values_vec and \p elem_dof_values
207 : */
208 : static void set_local_vectors(const NumericVector<Number> & global_vector,
209 : const std::vector<dof_id_type> & elem_dof_indices,
210 : std::vector<Number> & elem_dof_values_vec,
211 : EigenVector & elem_dof_values);
212 :
213 : /**
214 : * Takes an incoming "full" RHS from the solver, where full means the union of uncondensed and
215 : * condennsed dofs, and condenses it down into the condensed \p _reduced_rhs data member using
216 : * element Schur complements
217 : */
218 : void forward_elimination(const NumericVector<Number> & full_rhs);
219 :
220 : /**
221 : * After performing the solve with the \p _reduced_rhs and Schur complement matrix (\p
222 : * _reduced_sys_mat) to determine the \p _reduced_sol, we use the uncondensed \p full_rhs along
223 : * with the \p _reduced_sol to back substitute into the \p full_sol, which is the final output
224 : * data of the static condensation preconditioner application
225 : */
226 : void backwards_substitution(const NumericVector<Number> & full_rhs,
227 : NumericVector<Number> & full_sol);
228 :
229 : /**
230 : * Data stored on a per-element basis used to compute element Schur complements and their
231 : * applications to vectors
232 : */
233 : struct MatrixData
234 : {
235 : /// condensed-condensed matrix entries
236 : EigenMatrix Acc;
237 : /// condensed-uncondensed matrix entries
238 : EigenMatrix Acu;
239 : /// uncondensed-condensed matrix entries
240 : EigenMatrix Auc;
241 : /// uncondensed-uncondensed matrix entries
242 : EigenMatrix Auu;
243 :
244 : // Acc LU decompositions
245 : typename std::remove_const<decltype(Acc.partialPivLu())>::type AccFactor;
246 : };
247 :
248 : /// A map from element ID to Schur complement data
249 : std::unordered_map<dof_id_type, MatrixData> _elem_to_matrix_data;
250 :
251 : const MeshBase & _mesh;
252 : System & _system;
253 : const DofMap & _full_dof_map;
254 : StaticCondensationDofMap & _reduced_dof_map;
255 :
256 : /// global sparse matrix for the uncondensed degrees of freedom
257 : std::unique_ptr<SparseMatrix<Number>> _reduced_sys_mat;
258 : /// solution for the uncondensed degrees of freedom
259 : std::unique_ptr<NumericVector<Number>> _reduced_sol;
260 : /// RHS corresponding to the uncondensed degrees of freedom. This is constructed by applying the
261 : /// element Schur complements to an incoming RHS which contains both condensed and uncondensed
262 : /// degrees of freedom
263 : std::unique_ptr<NumericVector<Number>> _reduced_rhs;
264 : /// The solver for the uncondensed degrees of freedom
265 : std::unique_ptr<LinearSolver<Number>> _reduced_solver;
266 : /// This is a ghosted representation of the full (uncondensed + condensed) solution. Note that
267 : // this is, in general, *not equal* to the system solution, e.g. this may correspond to the
268 : // solution for the Newton *update*
269 : std::unique_ptr<NumericVector<Number>> _ghosted_full_sol;
270 :
271 : /// The current element ID. This is one half of a key, along with the global index, that maps to a
272 : /// local index
273 : dof_id_type _current_elem_id;
274 :
275 : /// Helper data member for adding individual matrix elements
276 : DenseMatrix<Number> _size_one_mat;
277 :
278 : /// Preconditioner object which will call back to us for the preconditioning action
279 : std::unique_ptr<StaticCondensationPreconditioner> _scp;
280 :
281 : /// Whether our object has been initialized
282 : bool _sc_is_initialized;
283 :
284 : /// Whether we have cached values via add_XXX()
285 : bool _have_cached_values;
286 :
287 : /// The parallel type to use for the reduced matrix
288 : ParallelType _parallel_type;
289 :
290 : /// whether this matrix represents uncondensed dofs only. In that case when building the Schur
291 : /// complement we won't attempt to invert zero element matrices corresponding to the condensed
292 : /// dofs
293 : bool _uncondensed_dofs_only;
294 : };
295 :
296 : inline const SparseMatrix<Number> & StaticCondensation::get_condensed_mat() const
297 : {
298 : libmesh_assert(_reduced_sys_mat);
299 : return *_reduced_sys_mat;
300 : }
301 :
302 : inline LinearSolver<Number> & StaticCondensation::reduced_system_solver()
303 : {
304 : libmesh_assert_msg(_reduced_solver, "Reduced system solver not built yet");
305 : return *_reduced_solver;
306 : }
307 :
308 : }
309 :
310 : #else
311 :
312 : #include "libmesh/sparse_matrix.h"
313 : #include <unordered_set>
314 :
315 : namespace libMesh
316 : {
317 : class MeshBase;
318 : class System;
319 : class DofMap;
320 : class StaticCondensationPreconditioner;
321 : class StaticCondensationDofMap;
322 :
323 : class StaticCondensation : public SparseMatrix<Number>
324 : {
325 : public:
326 : StaticCondensation(const MeshBase &,
327 : const System &,
328 : const DofMap & full_dof_map,
329 : StaticCondensationDofMap & reduced_dof_map);
330 :
331 : const std::unordered_set<unsigned int> & uncondensed_vars() const { libmesh_not_implemented(); }
332 0 : StaticCondensationPreconditioner & get_preconditioner() { libmesh_not_implemented(); }
333 0 : virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override
334 : {
335 0 : libmesh_not_implemented();
336 : }
337 0 : virtual SolverPackage solver_package() override { libmesh_not_implemented(); }
338 0 : virtual void init(const numeric_index_type,
339 : const numeric_index_type,
340 : const numeric_index_type,
341 : const numeric_index_type,
342 : const numeric_index_type = 30,
343 : const numeric_index_type = 10,
344 : const numeric_index_type = 1) override
345 : {
346 0 : libmesh_not_implemented();
347 : }
348 0 : virtual void init(ParallelType) override { libmesh_not_implemented(); }
349 0 : virtual void clear() override { libmesh_not_implemented(); }
350 0 : virtual void zero() override { libmesh_not_implemented(); }
351 0 : virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override
352 : {
353 0 : libmesh_not_implemented();
354 : }
355 0 : virtual std::unique_ptr<SparseMatrix<Number>> clone() const override
356 : {
357 0 : libmesh_not_implemented();
358 : }
359 0 : virtual void close() override { libmesh_not_implemented(); }
360 0 : virtual numeric_index_type m() const override { libmesh_not_implemented(); }
361 0 : virtual numeric_index_type n() const override { libmesh_not_implemented(); }
362 0 : virtual numeric_index_type row_start() const override { libmesh_not_implemented(); }
363 0 : virtual numeric_index_type row_stop() const override { libmesh_not_implemented(); }
364 0 : virtual numeric_index_type col_start() const override { libmesh_not_implemented(); }
365 0 : virtual numeric_index_type col_stop() const override { libmesh_not_implemented(); }
366 0 : virtual void set(const numeric_index_type, const numeric_index_type, const Number) override
367 : {
368 0 : libmesh_not_implemented();
369 : }
370 0 : virtual void add(const numeric_index_type, const numeric_index_type, const Number) override
371 : {
372 0 : libmesh_not_implemented();
373 : }
374 0 : virtual void add_matrix(const DenseMatrix<Number> &,
375 : const std::vector<numeric_index_type> &,
376 : const std::vector<numeric_index_type> &) override
377 : {
378 0 : libmesh_not_implemented();
379 : }
380 0 : virtual void add_matrix(const DenseMatrix<Number> &,
381 : const std::vector<numeric_index_type> &) override
382 : {
383 0 : libmesh_not_implemented();
384 : }
385 0 : virtual void add(const Number, const SparseMatrix<Number> &) override
386 : {
387 0 : libmesh_not_implemented();
388 : }
389 0 : virtual Number operator()(const numeric_index_type, const numeric_index_type) const override
390 : {
391 0 : libmesh_not_implemented();
392 : }
393 0 : virtual Real l1_norm() const override { libmesh_not_implemented(); }
394 0 : virtual Real linfty_norm() const override { libmesh_not_implemented(); }
395 0 : virtual bool closed() const override { libmesh_not_implemented(); }
396 0 : virtual void print_personal(std::ostream & = libMesh::out) const override
397 : {
398 0 : libmesh_not_implemented();
399 : }
400 0 : virtual void get_diagonal(NumericVector<Number> &) const override { libmesh_not_implemented(); }
401 0 : virtual void get_transpose(SparseMatrix<Number> &) const override { libmesh_not_implemented(); }
402 0 : virtual void get_row(numeric_index_type,
403 : std::vector<numeric_index_type> &,
404 : std::vector<Number> &) const override
405 : {
406 0 : libmesh_not_implemented();
407 : }
408 0 : void init() { libmesh_not_implemented(); }
409 0 : void setup() { libmesh_not_implemented(); }
410 0 : void apply(const NumericVector<Number> &, NumericVector<Number> &) { libmesh_not_implemented(); }
411 :
412 : void set_current_elem(const Elem &) {}
413 : void dont_condense_vars(const std::unordered_set<unsigned int> &) {}
414 : };
415 : }
416 :
417 : #endif // LIBMESH_HAVE_EIGEN && LIBMESH_HAVE_PETSC
418 : #endif // LIBMESH_STATIC_CONDENSATION_H
|
