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 100 : 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> & parent_rhs, NumericVector<Number> & parent_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 20 : StaticCondensationPreconditioner & get_preconditioner() { return *_scp; }
182 :
183 : /**
184 : * @returns The reduced system linear solver
185 : */
186 : LinearSolver<Number> & reduced_system_solver();
187 :
188 48 : 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 RHS from the solver on the parent space of the static condensation dof map,
215 : * which holds both uncondensed and condensed dofs, and condenses it down into the condensed \p
216 : * _reduced_rhs data member using element Schur complements
217 : */
218 : void forward_elimination(const NumericVector<Number> & parent_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 parent-space \p parent_rhs
223 : * along with the \p _reduced_sol to back substitute into the parent-space \p parent_sol, which
224 : * is the final output data of the static condensation preconditioner application
225 : */
226 : void backwards_substitution(const NumericVector<Number> & parent_rhs,
227 : NumericVector<Number> & parent_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 : StaticCondensationDofMap & _reduced_dof_map;
254 :
255 : /// global sparse matrix for the uncondensed degrees of freedom
256 : std::unique_ptr<SparseMatrix<Number>> _reduced_sys_mat;
257 : /// solution for the uncondensed degrees of freedom
258 : std::unique_ptr<NumericVector<Number>> _reduced_sol;
259 : /// RHS corresponding to the uncondensed degrees of freedom. This is constructed by applying the
260 : /// element Schur complements to an incoming RHS which contains both condensed and uncondensed
261 : /// degrees of freedom
262 : std::unique_ptr<NumericVector<Number>> _reduced_rhs;
263 : /// The solver for the uncondensed degrees of freedom
264 : std::unique_ptr<LinearSolver<Number>> _reduced_solver;
265 : /// This is a ghosted representation of the solution on the parent space, uncondensed and
266 : // condensed dofs together. Note that this is, in general, *not equal* to the system solution,
267 : // e.g. this may correspond to the solution for the Newton *update*
268 : std::unique_ptr<NumericVector<Number>> _ghosted_parent_sol;
269 :
270 : /// The current element ID. This is one half of a key, along with the global index, that maps to a
271 : /// local index
272 : dof_id_type _current_elem_id;
273 :
274 : /// Helper data member for adding individual matrix elements
275 : DenseMatrix<Number> _size_one_mat;
276 :
277 : /// Preconditioner object which will call back to us for the preconditioning action
278 : std::unique_ptr<StaticCondensationPreconditioner> _scp;
279 :
280 : /// Whether our object has been initialized
281 : bool _sc_is_initialized;
282 :
283 : /// Whether we have cached values via add_XXX()
284 : bool _have_cached_values;
285 :
286 : /// The parallel type to use for the reduced matrix
287 : ParallelType _parallel_type;
288 :
289 : /// whether this matrix represents uncondensed dofs only. In that case when building the Schur
290 : /// complement we won't attempt to invert zero element matrices corresponding to the condensed
291 : /// dofs
292 : bool _uncondensed_dofs_only;
293 : };
294 :
295 : inline const SparseMatrix<Number> & StaticCondensation::get_condensed_mat() const
296 : {
297 : libmesh_assert(_reduced_sys_mat);
298 : return *_reduced_sys_mat;
299 : }
300 :
301 : inline LinearSolver<Number> & StaticCondensation::reduced_system_solver()
302 : {
303 : libmesh_assert_msg(_reduced_solver, "Reduced system solver not built yet");
304 : return *_reduced_solver;
305 : }
306 :
307 : }
308 :
309 : #else
310 :
311 : #include "libmesh/sparse_matrix.h"
312 : #include <unordered_set>
313 :
314 : namespace libMesh
315 : {
316 : class MeshBase;
317 : class System;
318 : class DofMap;
319 : class StaticCondensationPreconditioner;
320 : class StaticCondensationDofMap;
321 :
322 : class StaticCondensation : public SparseMatrix<Number>
323 : {
324 : public:
325 : StaticCondensation(const MeshBase &,
326 : const System &,
327 : const DofMap & full_dof_map,
328 : StaticCondensationDofMap & reduced_dof_map);
329 :
330 : const std::unordered_set<unsigned int> & uncondensed_vars() const { libmesh_not_implemented(); }
331 0 : StaticCondensationPreconditioner & get_preconditioner() { libmesh_not_implemented(); }
332 0 : virtual SparseMatrix<Number> & operator=(const SparseMatrix<Number> &) override
333 : {
334 0 : libmesh_not_implemented();
335 : }
336 0 : virtual SolverPackage solver_package() override { libmesh_not_implemented(); }
337 0 : virtual void init(const numeric_index_type,
338 : const numeric_index_type,
339 : const numeric_index_type,
340 : const numeric_index_type,
341 : const numeric_index_type = 30,
342 : const numeric_index_type = 10,
343 : const numeric_index_type = 1) override
344 : {
345 0 : libmesh_not_implemented();
346 : }
347 0 : virtual void init(ParallelType) override { libmesh_not_implemented(); }
348 0 : virtual void clear() override { libmesh_not_implemented(); }
349 0 : virtual void zero() override { libmesh_not_implemented(); }
350 0 : virtual std::unique_ptr<SparseMatrix<Number>> zero_clone() const override
351 : {
352 0 : libmesh_not_implemented();
353 : }
354 0 : virtual std::unique_ptr<SparseMatrix<Number>> clone() const override
355 : {
356 0 : libmesh_not_implemented();
357 : }
358 0 : virtual void close() override { libmesh_not_implemented(); }
359 0 : virtual numeric_index_type m() const override { libmesh_not_implemented(); }
360 0 : virtual numeric_index_type n() const override { libmesh_not_implemented(); }
361 0 : virtual numeric_index_type row_start() const override { libmesh_not_implemented(); }
362 0 : virtual numeric_index_type row_stop() const override { libmesh_not_implemented(); }
363 0 : virtual numeric_index_type col_start() const override { libmesh_not_implemented(); }
364 0 : virtual numeric_index_type col_stop() const override { libmesh_not_implemented(); }
365 0 : virtual void set(const numeric_index_type, const numeric_index_type, const Number) override
366 : {
367 0 : libmesh_not_implemented();
368 : }
369 0 : virtual void add(const numeric_index_type, const numeric_index_type, const Number) override
370 : {
371 0 : libmesh_not_implemented();
372 : }
373 0 : virtual void add_matrix(const DenseMatrix<Number> &,
374 : const std::vector<numeric_index_type> &,
375 : const std::vector<numeric_index_type> &) override
376 : {
377 0 : libmesh_not_implemented();
378 : }
379 0 : virtual void add_matrix(const DenseMatrix<Number> &,
380 : const std::vector<numeric_index_type> &) override
381 : {
382 0 : libmesh_not_implemented();
383 : }
384 0 : virtual void add(const Number, const SparseMatrix<Number> &) override
385 : {
386 0 : libmesh_not_implemented();
387 : }
388 0 : virtual Number operator()(const numeric_index_type, const numeric_index_type) const override
389 : {
390 0 : libmesh_not_implemented();
391 : }
392 0 : virtual Real l1_norm() const override { libmesh_not_implemented(); }
393 0 : virtual Real linfty_norm() const override { libmesh_not_implemented(); }
394 0 : virtual bool closed() const override { libmesh_not_implemented(); }
395 0 : virtual void print_personal(std::ostream & = libMesh::out) const override
396 : {
397 0 : libmesh_not_implemented();
398 : }
399 0 : virtual void get_diagonal(NumericVector<Number> &) const override { libmesh_not_implemented(); }
400 0 : virtual void get_transpose(SparseMatrix<Number> &) const override { libmesh_not_implemented(); }
401 0 : virtual void get_row(numeric_index_type,
402 : std::vector<numeric_index_type> &,
403 : std::vector<Number> &) const override
404 : {
405 0 : libmesh_not_implemented();
406 : }
407 0 : void init() { libmesh_not_implemented(); }
408 0 : void setup() { libmesh_not_implemented(); }
409 0 : void apply(const NumericVector<Number> &, NumericVector<Number> &) { libmesh_not_implemented(); }
410 :
411 : void set_current_elem(const Elem &) {}
412 : void dont_condense_vars(const std::unordered_set<unsigned int> &) {}
413 : };
414 : }
415 :
416 : #endif // LIBMESH_HAVE_EIGEN && LIBMESH_HAVE_PETSC
417 : #endif // LIBMESH_STATIC_CONDENSATION_H
|
