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 : // Local includes
19 : #include "libmesh/dof_map.h"
20 :
21 : // libMesh includes
22 : #include "libmesh/boundary_info.h" // needed for dirichlet constraints
23 : #include "libmesh/dense_matrix.h"
24 : #include "libmesh/dense_vector.h"
25 : #include "libmesh/dirichlet_boundaries.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/elem_range.h"
28 : #include "libmesh/fe_base.h"
29 : #include "libmesh/fe_interface.h"
30 : #include "libmesh/fe_type.h"
31 : #include "libmesh/function_base.h"
32 : #include "libmesh/int_range.h"
33 : #include "libmesh/libmesh_logging.h"
34 : #include "libmesh/linear_solver.h" // for spline Dirichlet projection solves
35 : #include "libmesh/mesh_base.h"
36 : #include "libmesh/null_output_iterator.h"
37 : #include "libmesh/mesh_tools.h" // for libmesh_assert_valid_boundary_ids()
38 : #include "libmesh/nonlinear_implicit_system.h"
39 : #include "libmesh/numeric_vector.h" // for enforce_constraints_exactly()
40 : #include "libmesh/parallel_algebra.h"
41 : #include "libmesh/parallel_elem.h"
42 : #include "libmesh/parallel_node.h"
43 : #include "libmesh/periodic_boundaries.h"
44 : #include "libmesh/periodic_boundary.h"
45 : #include "libmesh/periodic_boundary_base.h"
46 : #include "libmesh/point_locator_base.h"
47 : #include "libmesh/quadrature.h" // for dirichlet constraints
48 : #include "libmesh/raw_accessor.h"
49 : #include "libmesh/sparse_matrix.h" // needed to constrain adjoint rhs
50 : #include "libmesh/static_condensation_dof_map.h"
51 : #include "libmesh/system.h" // needed by enforce_constraints_exactly()
52 : #include "libmesh/tensor_tools.h"
53 : #include "libmesh/threads.h"
54 : #include "libmesh/enum_to_string.h"
55 : #include "libmesh/coupling_matrix.h"
56 :
57 : // TIMPI includes
58 : #include "timpi/parallel_implementation.h"
59 : #include "timpi/parallel_sync.h"
60 :
61 : // C++ Includes
62 : #include <set>
63 : #include <algorithm> // for std::count, std::fill
64 : #include <sstream>
65 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
66 : #include <cmath>
67 : #include <memory>
68 : #include <numeric>
69 : #include <unordered_set>
70 :
71 : // Anonymous namespace to hold helper classes
72 : namespace {
73 :
74 : using namespace libMesh;
75 :
76 : class ComputeConstraints
77 : {
78 : public:
79 9406 : ComputeConstraints (DofConstraints & constraints,
80 : DofMap & dof_map,
81 : #ifdef LIBMESH_ENABLE_PERIODIC
82 : PeriodicBoundaries & periodic_boundaries,
83 : #endif
84 : const MeshBase & mesh,
85 329528 : const unsigned int variable_number) :
86 310716 : _constraints(constraints),
87 310716 : _dof_map(dof_map),
88 : #ifdef LIBMESH_ENABLE_PERIODIC
89 310716 : _periodic_boundaries(periodic_boundaries),
90 : #endif
91 310716 : _mesh(mesh),
92 329528 : _variable_number(variable_number)
93 9406 : {}
94 :
95 329731 : void operator()(const ConstElemRange & range) const
96 : {
97 329731 : const Variable & var_description = _dof_map.variable(_variable_number);
98 :
99 : #ifdef LIBMESH_ENABLE_PERIODIC
100 320258 : std::unique_ptr<PointLocatorBase> point_locator;
101 : const bool have_periodic_boundaries =
102 329731 : !_periodic_boundaries.empty();
103 329731 : if (have_periodic_boundaries && !range.empty())
104 1414 : point_locator = _mesh.sub_point_locator();
105 : #endif
106 :
107 8185314 : for (const auto & elem : range)
108 7855583 : if (var_description.active_on_subdomain(elem->subdomain_id()))
109 : {
110 : #ifdef LIBMESH_ENABLE_AMR
111 8537681 : FEInterface::compute_constraints (_constraints,
112 : _dof_map,
113 7841055 : _variable_number,
114 : elem);
115 : #endif
116 : #ifdef LIBMESH_ENABLE_PERIODIC
117 : // FIXME: periodic constraints won't work on a non-serial
118 : // mesh unless it's kept ghost elements from opposing
119 : // boundaries!
120 7841055 : if (have_periodic_boundaries)
121 532492 : FEInterface::compute_periodic_constraints (_constraints,
122 : _dof_map,
123 264164 : _periodic_boundaries,
124 : _mesh,
125 67082 : point_locator.get(),
126 264164 : _variable_number,
127 : elem);
128 : #endif
129 : }
130 329731 : }
131 :
132 : private:
133 : DofConstraints & _constraints;
134 : DofMap & _dof_map;
135 : #ifdef LIBMESH_ENABLE_PERIODIC
136 : PeriodicBoundaries & _periodic_boundaries;
137 : #endif
138 : const MeshBase & _mesh;
139 : const unsigned int _variable_number;
140 : };
141 :
142 :
143 :
144 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
145 : class ComputeNodeConstraints
146 : {
147 : public:
148 8448 : ComputeNodeConstraints (NodeConstraints & node_constraints,
149 : #ifdef LIBMESH_ENABLE_PERIODIC
150 : PeriodicBoundaries & periodic_boundaries,
151 : #endif
152 16896 : const MeshBase & mesh) :
153 : _node_constraints(node_constraints),
154 : #ifdef LIBMESH_ENABLE_PERIODIC
155 : _periodic_boundaries(periodic_boundaries),
156 : #endif
157 16896 : _mesh(mesh)
158 8448 : {}
159 :
160 16986 : void operator()(const ConstElemRange & range) const
161 : {
162 : #ifdef LIBMESH_ENABLE_PERIODIC
163 8493 : std::unique_ptr<PointLocatorBase> point_locator;
164 16986 : bool have_periodic_boundaries = !_periodic_boundaries.empty();
165 16986 : if (have_periodic_boundaries && !range.empty())
166 244 : point_locator = _mesh.sub_point_locator();
167 : #endif
168 :
169 1052142 : for (const auto & elem : range)
170 : {
171 : #ifdef LIBMESH_ENABLE_AMR
172 1035156 : FEBase::compute_node_constraints (_node_constraints, elem);
173 : #endif
174 : #ifdef LIBMESH_ENABLE_PERIODIC
175 : // FIXME: periodic constraints won't work on a non-serial
176 : // mesh unless it's kept ghost elements from opposing
177 : // boundaries!
178 1035156 : if (have_periodic_boundaries)
179 324310 : FEBase::compute_periodic_node_constraints (_node_constraints,
180 129724 : _periodic_boundaries,
181 : _mesh,
182 64862 : point_locator.get(),
183 : elem);
184 : #endif
185 : }
186 16986 : }
187 :
188 : private:
189 : NodeConstraints & _node_constraints;
190 : #ifdef LIBMESH_ENABLE_PERIODIC
191 : PeriodicBoundaries & _periodic_boundaries;
192 : #endif
193 : const MeshBase & _mesh;
194 : };
195 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
196 :
197 :
198 : #ifdef LIBMESH_ENABLE_DIRICHLET
199 :
200 : /**
201 : * This functor class hierarchy adds a constraint row to a DofMap
202 : */
203 : class AddConstraint
204 : {
205 : protected:
206 : DofMap & dof_map;
207 :
208 : public:
209 20314 : AddConstraint(DofMap & dof_map_in) : dof_map(dof_map_in) {}
210 :
211 : virtual void operator()(dof_id_type dof_number,
212 : const DofConstraintRow & constraint_row,
213 : const Number constraint_rhs) const = 0;
214 : };
215 :
216 : class AddPrimalConstraint : public AddConstraint
217 : {
218 : public:
219 18129 : AddPrimalConstraint(DofMap & dof_map_in) : AddConstraint(dof_map_in) {}
220 :
221 649521 : virtual void operator()(dof_id_type dof_number,
222 : const DofConstraintRow & constraint_row,
223 : const Number constraint_rhs) const
224 : {
225 649521 : if (!dof_map.is_constrained_dof(dof_number))
226 320173 : dof_map.add_constraint_row (dof_number, constraint_row,
227 : constraint_rhs, true);
228 649521 : }
229 : };
230 :
231 : class AddAdjointConstraint : public AddConstraint
232 : {
233 : private:
234 : const unsigned int qoi_index;
235 :
236 : public:
237 62 : AddAdjointConstraint(DofMap & dof_map_in, unsigned int qoi_index_in)
238 2247 : : AddConstraint(dof_map_in), qoi_index(qoi_index_in) {}
239 :
240 14304 : virtual void operator()(dof_id_type dof_number,
241 : const DofConstraintRow & constraint_row,
242 : const Number constraint_rhs) const
243 : {
244 14304 : dof_map.add_adjoint_constraint_row
245 14304 : (qoi_index, dof_number, constraint_row, constraint_rhs,
246 : true);
247 14304 : }
248 : };
249 :
250 :
251 :
252 : /**
253 : * This class implements turning an arbitrary
254 : * boundary function into Dirichlet constraints. It
255 : * may be executed in parallel on multiple threads.
256 : */
257 : class ConstrainDirichlet
258 : {
259 : private:
260 : DofMap & dof_map;
261 : const MeshBase & mesh;
262 : const Real time;
263 : const DirichletBoundaries & dirichlets;
264 :
265 : const AddConstraint & add_fn;
266 :
267 664676 : static Number f_component (FunctionBase<Number> * f,
268 : FEMFunctionBase<Number> * f_fem,
269 : const FEMContext * c,
270 : unsigned int i,
271 : const Point & p,
272 : Real time)
273 : {
274 664676 : if (f_fem)
275 : {
276 0 : if (c)
277 0 : return f_fem->component(*c, i, p, time);
278 : else
279 0 : return std::numeric_limits<Real>::quiet_NaN();
280 : }
281 664676 : return f->component(i, p, time);
282 : }
283 :
284 30966 : static Gradient g_component (FunctionBase<Gradient> * g,
285 : FEMFunctionBase<Gradient> * g_fem,
286 : const FEMContext * c,
287 : unsigned int i,
288 : const Point & p,
289 : Real time)
290 : {
291 30966 : if (g_fem)
292 : {
293 0 : if (c)
294 0 : return g_fem->component(*c, i, p, time);
295 : else
296 0 : return std::numeric_limits<Number>::quiet_NaN();
297 : }
298 30966 : return g->component(i, p, time);
299 : }
300 :
301 :
302 :
303 : /**
304 : * Handy struct to pass around BoundaryInfo for a single Elem. Must
305 : * be created with a reference to a BoundaryInfo object and a map
306 : * from boundary_id -> set<DirichletBoundary *> objects involving
307 : * that id.
308 : */
309 : struct SingleElemBoundaryInfo
310 : {
311 19784 : SingleElemBoundaryInfo(const BoundaryInfo & bi,
312 20380 : const std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>> & ordered_map_in) :
313 19188 : boundary_info(bi),
314 19188 : boundary_id_to_ordered_dirichlet_boundaries(ordered_map_in),
315 19188 : elem(nullptr),
316 19188 : n_sides(0),
317 19188 : n_edges(0),
318 20380 : n_nodes(0)
319 19784 : {}
320 :
321 : const BoundaryInfo & boundary_info;
322 : const std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>> & boundary_id_to_ordered_dirichlet_boundaries;
323 : const Elem * elem;
324 :
325 : unsigned short n_sides;
326 : unsigned short n_edges;
327 : unsigned short n_nodes;
328 :
329 : // Mapping from DirichletBoundary objects which are active on this
330 : // element to sides/nodes/edges/shellfaces of this element which
331 : // they are active on.
332 : std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_node_map;
333 : std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_side_map;
334 : std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_edge_map;
335 : std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_shellface_map;
336 :
337 : std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_nodeset_map;
338 :
339 : // The set of (dirichlet_id, DirichletBoundary) pairs which have at least one boundary
340 : // id related to this Elem.
341 : std::set<std::pair<unsigned int, DirichletBoundary *>> ordered_dbs;
342 :
343 : /**
344 : * Given a single Elem, fills the SingleElemBoundaryInfo struct with
345 : * required data.
346 : *
347 : * @return true if this Elem has _any_ boundary ids associated with
348 : * it, false otherwise.
349 : */
350 1171243 : bool reinit(const Elem * elem_in)
351 : {
352 1171243 : elem = elem_in;
353 :
354 1171243 : n_sides = elem->n_sides();
355 1171243 : n_edges = elem->n_edges();
356 1171243 : n_nodes = elem->n_nodes();
357 :
358 : // objects and node/side/edge/shellface ids.
359 100487 : is_boundary_node_map.clear();
360 100487 : is_boundary_side_map.clear();
361 100487 : is_boundary_edge_map.clear();
362 100487 : is_boundary_shellface_map.clear();
363 100487 : is_boundary_nodeset_map.clear();
364 :
365 : // Clear any DirichletBoundaries from the previous Elem
366 100487 : ordered_dbs.clear();
367 :
368 : // Update has_dirichlet_constraint below, and if it remains false then
369 : // we can skip this element since there are not constraints to impose.
370 100487 : bool has_dirichlet_constraint = false;
371 :
372 : // Container to catch boundary ids handed back for sides,
373 : // nodes, and edges in the loops below.
374 100487 : std::vector<boundary_id_type> ids_vec;
375 :
376 6219454 : for (unsigned char s = 0; s != n_sides; ++s)
377 : {
378 : // First see if this side has been requested
379 5048211 : boundary_info.boundary_ids (elem, s, ids_vec);
380 :
381 432115 : bool do_this_side = false;
382 5300847 : for (const auto & bc_id : ids_vec)
383 : {
384 273988 : if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
385 21352 : it != boundary_id_to_ordered_dirichlet_boundaries.end())
386 : {
387 13154 : do_this_side = true;
388 :
389 : // Associate every DirichletBoundary object that has this bc_id with the current Elem
390 154338 : ordered_dbs.insert(it->second.begin(), it->second.end());
391 :
392 : // Turn on the flag for the current side for each DirichletBoundary
393 316969 : for (const auto & db_pair : it->second)
394 : {
395 : // Attempt to emplace an empty vector. If there
396 : // is already an entry, the insertion will fail
397 : // and we'll get an iterator back to the
398 : // existing entry. Either way, we'll then set
399 : // index s of that vector to "true".
400 297500 : auto pr = is_boundary_side_map.emplace(db_pair.second, std::vector<bool>(n_sides, false));
401 27762 : pr.first->second[s] = true;
402 : }
403 : }
404 : }
405 :
406 5048211 : if (!do_this_side)
407 4474912 : continue;
408 :
409 13154 : has_dirichlet_constraint = true;
410 :
411 : // Then determine what nodes are on this side
412 1380023 : for (unsigned int n = 0; n != n_nodes; ++n)
413 1225685 : if (elem->is_node_on_side(n,s))
414 : {
415 : // Attempt to emplace an empty vector. If there is
416 : // already an entry, the insertion will fail and we'll
417 : // get an iterator back to the existing entry. Either
418 : // way, we'll then set index n of that vector to
419 : // "true".
420 1091891 : for (const auto & db_pair : ordered_dbs)
421 : {
422 : // Only add this as a boundary node for this db if
423 : // it is also a boundary side for this db.
424 559687 : if (auto side_it = is_boundary_side_map.find(db_pair.second);
425 559687 : side_it != is_boundary_side_map.end() && side_it->second[s])
426 : {
427 1020116 : auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
428 94434 : pr.first->second[n] = true;
429 : }
430 : }
431 : }
432 :
433 : // Finally determine what edges are on this side
434 1288749 : for (unsigned int e = 0; e != n_edges; ++e)
435 1134411 : if (elem->is_edge_on_side(e,s))
436 : {
437 : // Attempt to emplace an empty vector. If there is
438 : // already an entry, the insertion will fail and we'll
439 : // get an iterator back to the existing entry. Either
440 : // way, we'll then set index e of that vector to
441 : // "true".
442 737841 : for (const auto & db_pair : ordered_dbs)
443 : {
444 : // Only add this as a boundary edge for this db if
445 : // it is also a boundary side for this db.
446 377667 : if (auto side_it = is_boundary_side_map.find(db_pair.second);
447 377667 : side_it != is_boundary_side_map.end() && side_it->second[s])
448 : {
449 690132 : auto pr = is_boundary_edge_map.emplace(db_pair.second, std::vector<bool>(n_edges, false));
450 63530 : pr.first->second[e] = true;
451 : }
452 : }
453 : }
454 : } // for (s = 0..n_sides)
455 :
456 : // We can also impose Dirichlet boundary conditions on nodes, so we should
457 : // also independently check whether the nodes have been requested
458 8963730 : for (unsigned int n=0; n != n_nodes; ++n)
459 : {
460 8454918 : boundary_info.boundary_ids (elem->node_ptr(n), ids_vec);
461 :
462 8356291 : for (const auto & bc_id : ids_vec)
463 : {
464 611011 : if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
465 47207 : it != boundary_id_to_ordered_dirichlet_boundaries.end())
466 : {
467 : // Associate every DirichletBoundary object that has this bc_id with the current Elem
468 215612 : ordered_dbs.insert(it->second.begin(), it->second.end());
469 :
470 : // Turn on the flag for the current node for each DirichletBoundary
471 448942 : for (const auto & db_pair : it->second)
472 : {
473 427306 : auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
474 39354 : pr.first->second[n] = true;
475 :
476 427306 : auto pr2 = is_boundary_nodeset_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
477 39354 : pr2.first->second[n] = true;
478 : }
479 :
480 18153 : has_dirichlet_constraint = true;
481 : }
482 : }
483 : } // for (n = 0..n_nodes)
484 :
485 : // We can also impose Dirichlet boundary conditions on edges, so we should
486 : // also independently check whether the edges have been requested
487 8052034 : for (unsigned short e=0; e != n_edges; ++e)
488 : {
489 6880791 : boundary_info.edge_boundary_ids (elem, e, ids_vec);
490 :
491 586235 : bool do_this_side = false;
492 6881031 : for (const auto & bc_id : ids_vec)
493 : {
494 260 : if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
495 20 : it != boundary_id_to_ordered_dirichlet_boundaries.end())
496 : {
497 20 : do_this_side = true;
498 :
499 : // We need to loop over all DirichletBoundary objects associated with bc_id
500 240 : ordered_dbs.insert(it->second.begin(), it->second.end());
501 :
502 : // Turn on the flag for the current edge for each DirichletBoundary
503 768 : for (const auto & db_pair : it->second)
504 : {
505 968 : auto pr = is_boundary_edge_map.emplace(db_pair.second, std::vector<bool>(n_edges, false));
506 88 : pr.first->second[e] = true;
507 : }
508 : }
509 : }
510 :
511 6880791 : if (!do_this_side)
512 6294336 : continue;
513 :
514 20 : has_dirichlet_constraint = true;
515 :
516 : // Then determine what nodes are on this edge
517 2160 : for (unsigned int n = 0; n != n_nodes; ++n)
518 1920 : if (elem->is_node_on_edge(n,e))
519 : {
520 : // Attempt to emplace an empty vector. If there is
521 : // already an entry, the insertion will fail and we'll
522 : // get an iterator back to the existing entry. Either
523 : // way, we'll then set index n of that vector to
524 : // "true".
525 1536 : for (const auto & db_pair : ordered_dbs)
526 : {
527 : // Only add this as a boundary node for this db if
528 : // it is also a boundary edge for this db.
529 1056 : if (auto edge_it = is_boundary_edge_map.find(db_pair.second);
530 1056 : edge_it != is_boundary_edge_map.end() && edge_it->second[e])
531 : {
532 1936 : auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
533 176 : pr.first->second[n] = true;
534 : }
535 : }
536 : }
537 : }
538 :
539 : // We can also impose Dirichlet boundary conditions on shellfaces, so we should
540 : // also independently check whether the shellfaces have been requested
541 3513729 : for (unsigned short shellface=0; shellface != 2; ++shellface)
542 : {
543 2342486 : boundary_info.shellface_boundary_ids (elem, shellface, ids_vec);
544 200974 : bool do_this_shellface = false;
545 :
546 2382446 : for (const auto & bc_id : ids_vec)
547 : {
548 43290 : if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
549 3330 : it != boundary_id_to_ordered_dirichlet_boundaries.end())
550 : {
551 1 : has_dirichlet_constraint = true;
552 1 : do_this_shellface = true;
553 :
554 : // We need to loop over all DirichletBoundary objects associated with bc_id
555 12 : ordered_dbs.insert(it->second.begin(), it->second.end());
556 :
557 : // Turn on the flag for the current shellface for each DirichletBoundary
558 24 : for (const auto & db_pair : it->second)
559 : {
560 22 : auto pr = is_boundary_shellface_map.emplace(db_pair.second, std::vector<bool>(/*n_shellfaces=*/2, false));
561 2 : pr.first->second[shellface] = true;
562 : }
563 : }
564 : }
565 :
566 2342486 : if (do_this_shellface)
567 : {
568 : // Shellface BCs induce BCs on all the nodes of a shell Elem
569 60 : for (unsigned int n = 0; n != n_nodes; ++n)
570 96 : for (const auto & db_pair : ordered_dbs)
571 : {
572 : // Only add this as a boundary node for this db if
573 : // it is also a boundary shellface for this db.
574 48 : if (auto side_it = is_boundary_shellface_map.find(db_pair.second);
575 48 : side_it != is_boundary_shellface_map.end() && side_it->second[shellface])
576 : {
577 88 : auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
578 8 : pr.first->second[n] = true;
579 : }
580 : }
581 : }
582 : } // for (shellface = 0..2)
583 :
584 1271730 : return has_dirichlet_constraint;
585 : } // SingleElemBoundaryInfo::reinit()
586 :
587 : }; // struct SingleElemBoundaryInfo
588 :
589 :
590 :
591 : template<typename OutputType>
592 166234 : void apply_lagrange_dirichlet_impl(const SingleElemBoundaryInfo & sebi,
593 : const Variable & variable,
594 : const DirichletBoundary & dirichlet,
595 : FEMContext & fem_context) const
596 : {
597 : // Get pointer to the Elem we are currently working on
598 166234 : const Elem * elem = sebi.elem;
599 :
600 : // Per-subdomain variables don't need to be projected on
601 : // elements where they're not active
602 166234 : if (!variable.active_on_subdomain(elem->subdomain_id()))
603 0 : return;
604 :
605 14264 : FunctionBase<Number> * f = dirichlet.f.get();
606 14264 : FEMFunctionBase<Number> * f_fem = dirichlet.f_fem.get();
607 :
608 166234 : const System * f_system = dirichlet.f_system;
609 :
610 : // We need data to project
611 14264 : libmesh_assert(f || f_fem);
612 14264 : libmesh_assert(!(f && f_fem));
613 :
614 : // Iff our data depends on a system, we should have it.
615 14264 : libmesh_assert(!(f && f_system));
616 14264 : libmesh_assert(!(f_fem && !f_system));
617 :
618 : // The new element coefficients. For Lagrange FEs, these are the
619 : // nodal values.
620 166234 : DenseVector<Number> Ue;
621 :
622 : // Get a reference to the fe_type associated with this variable
623 14264 : const FEType & fe_type = variable.type();
624 :
625 : // Dimension of the vector-valued FE (1 for scalar-valued FEs)
626 166234 : unsigned int n_vec_dim = FEInterface::n_vec_dim(mesh, fe_type);
627 :
628 : const unsigned int var_component =
629 28528 : variable.first_scalar_number();
630 :
631 : // Get this Variable's number, as determined by the System.
632 28528 : const unsigned int var = variable.number();
633 :
634 : // If our supplied functions require a FEMContext, and if we have
635 : // an initialized solution to use with that FEMContext, then
636 : // create one. We're not going to use elem_jacobian or
637 : // subjacobians here so don't allocate them.
638 151970 : std::unique_ptr<FEMContext> context;
639 166234 : if (f_fem)
640 : {
641 0 : libmesh_assert(f_system);
642 0 : if (f_system->current_local_solution->initialized())
643 : {
644 0 : context = std::make_unique<FEMContext>
645 0 : (*f_system, nullptr,
646 0 : /* allocate local_matrices = */ false);
647 0 : f_fem->init_context(*context);
648 : }
649 : }
650 :
651 166234 : if (f_system && context.get())
652 0 : context->pre_fe_reinit(*f_system, elem);
653 :
654 : // Also pre-init the fem_context() we were passed on the current Elem.
655 166234 : fem_context.pre_fe_reinit(fem_context.get_system(), elem);
656 :
657 : // Get a reference to the DOF indices for the current element
658 : const std::vector<dof_id_type> & dof_indices =
659 14264 : fem_context.get_dof_indices(var);
660 :
661 : // The number of DOFs on the element
662 : const unsigned int n_dofs =
663 28528 : cast_int<unsigned int>(dof_indices.size());
664 :
665 : // Fixed vs. free DoFs on edge/face projections
666 180498 : std::vector<char> dof_is_fixed(n_dofs, false); // bools
667 :
668 : // Zero the interpolated values
669 151970 : Ue.resize (n_dofs); Ue.zero();
670 :
671 : // For Lagrange elements, side, edge, and shellface BCs all
672 : // "induce" boundary conditions on the nodes of those entities.
673 : // In SingleElemBoundaryInfo::reinit(), we therefore set entries
674 : // in the "is_boundary_node_map" container based on side and
675 : // shellface BCs, Then, when we actually apply constraints, we
676 : // only have to check whether any Nodes are in this container, and
677 : // compute values as necessary.
678 14264 : unsigned int current_dof = 0;
679 1483675 : for (unsigned int n=0; n!= sebi.n_nodes; ++n)
680 : {
681 : // For Lagrange this can return 0 (in case of a lower-order FE
682 : // on a higher-order Elem) or 1. This function accounts for the
683 : // Elem::p_level() internally.
684 112367 : const unsigned int nc =
685 1317441 : FEInterface::n_dofs_at_node (fe_type, elem, n);
686 :
687 : // If there are no DOFs at this node, then it doesn't matter
688 : // if it's technically a boundary node or not, there's nothing
689 : // to constrain.
690 1317441 : if (!nc)
691 101308 : continue;
692 :
693 : // Check whether the current node is a boundary node
694 1271037 : auto is_boundary_node_it = sebi.is_boundary_node_map.find(&dirichlet);
695 108481 : const bool is_boundary_node =
696 1379518 : (is_boundary_node_it != sebi.is_boundary_node_map.end() &&
697 108481 : is_boundary_node_it->second[n]);
698 :
699 : // Check whether the current node is in a boundary nodeset
700 1271037 : auto is_boundary_nodeset_it = sebi.is_boundary_nodeset_map.find(&dirichlet);
701 108481 : const bool is_boundary_nodeset =
702 1337589 : (is_boundary_nodeset_it != sebi.is_boundary_nodeset_map.end() &&
703 66552 : is_boundary_nodeset_it->second[n]);
704 :
705 : // If node is neither a boundary node or from a boundary nodeset, go to the next one.
706 1271037 : if ( !(is_boundary_node || is_boundary_nodeset) )
707 : {
708 687967 : current_dof += nc;
709 687967 : continue;
710 : }
711 :
712 : // Compute function values, storing them in Ue
713 49691 : libmesh_assert_equal_to (nc, n_vec_dim);
714 1166140 : for (unsigned int c = 0; c < n_vec_dim; c++)
715 : {
716 632761 : Ue(current_dof+c) =
717 583070 : f_component(f, f_fem, context.get(), var_component+c,
718 583070 : elem->point(n), time);
719 632761 : dof_is_fixed[current_dof+c] = true;
720 : }
721 583070 : current_dof += n_vec_dim;
722 : } // end for (n=0..n_nodes)
723 :
724 : // Lock the DofConstraints since it is shared among threads.
725 : {
726 28528 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
727 :
728 1437271 : for (unsigned int i = 0; i < n_dofs; i++)
729 : {
730 216962 : DofConstraintRow empty_row;
731 1379518 : if (dof_is_fixed[i] && !libmesh_isnan(Ue(i)))
732 632761 : add_fn (dof_indices[i], empty_row, Ue(i));
733 : }
734 : }
735 :
736 137706 : } // apply_lagrange_dirichlet_impl
737 :
738 :
739 :
740 : template<typename OutputType>
741 13501 : void apply_dirichlet_impl(const SingleElemBoundaryInfo & sebi,
742 : const Variable & variable,
743 : const DirichletBoundary & dirichlet,
744 : FEMContext & fem_context) const
745 : {
746 : // Get pointer to the Elem we are currently working on
747 13501 : const Elem * elem = sebi.elem;
748 :
749 : // Per-subdomain variables don't need to be projected on
750 : // elements where they're not active
751 13501 : if (!variable.active_on_subdomain(elem->subdomain_id()))
752 0 : return;
753 :
754 : typedef OutputType OutputShape;
755 : typedef typename TensorTools::IncrementRank<OutputShape>::type OutputGradient;
756 : //typedef typename TensorTools::IncrementRank<OutputGradient>::type OutputTensor;
757 : typedef typename TensorTools::MakeNumber<OutputShape>::type OutputNumber;
758 : typedef typename TensorTools::IncrementRank<OutputNumber>::type OutputNumberGradient;
759 : //typedef typename TensorTools::IncrementRank<OutputNumberGradient>::type OutputNumberTensor;
760 :
761 1127 : FunctionBase<Number> * f = dirichlet.f.get();
762 1127 : FunctionBase<Gradient> * g = dirichlet.g.get();
763 :
764 1127 : FEMFunctionBase<Number> * f_fem = dirichlet.f_fem.get();
765 1127 : FEMFunctionBase<Gradient> * g_fem = dirichlet.g_fem.get();
766 :
767 13501 : const System * f_system = dirichlet.f_system;
768 :
769 : // We need data to project
770 1127 : libmesh_assert(f || f_fem);
771 1127 : libmesh_assert(!(f && f_fem));
772 :
773 : // Iff our data depends on a system, we should have it.
774 1127 : libmesh_assert(!(f && f_system));
775 1127 : libmesh_assert(!(f_fem && !f_system));
776 :
777 : // The element matrix and RHS for projections.
778 : // Note that Ke is always real-valued, whereas
779 : // Fe may be complex valued if complex number
780 : // support is enabled
781 15755 : DenseMatrix<Real> Ke;
782 13501 : DenseVector<Number> Fe;
783 : // The new element coefficients
784 13501 : DenseVector<Number> Ue;
785 :
786 : // The dimensionality of the current mesh
787 13501 : const unsigned int dim = mesh.mesh_dimension();
788 :
789 : // Get a reference to the fe_type associated with this variable
790 1127 : const FEType & fe_type = variable.type();
791 :
792 : // Dimension of the vector-valued FE (1 for scalar-valued FEs)
793 13501 : unsigned int n_vec_dim = FEInterface::n_vec_dim(mesh, fe_type);
794 :
795 : const unsigned int var_component =
796 2254 : variable.first_scalar_number();
797 :
798 : // Get this Variable's number, as determined by the System.
799 2254 : const unsigned int var = variable.number();
800 :
801 : // The type of projections done depend on the FE's continuity.
802 13501 : FEContinuity cont = FEInterface::get_continuity(fe_type);
803 :
804 : // Make sure we have the right data available for C1 projections
805 1127 : if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
806 : {
807 : // We'll need gradient data for a C1 projection
808 459 : libmesh_assert(g || g_fem);
809 :
810 : // We currently demand that either neither nor both function
811 : // object depend on current FEM data.
812 459 : libmesh_assert(!(g && g_fem));
813 459 : libmesh_assert(!(f && g_fem));
814 459 : libmesh_assert(!(f_fem && g));
815 : }
816 :
817 : // If our supplied functions require a FEMContext, and if we have
818 : // an initialized solution to use with that FEMContext, then
819 : // create one. We're not going to use elem_jacobian or
820 : // subjacobians here so don't allocate them.
821 12374 : std::unique_ptr<FEMContext> context;
822 13501 : if (f_fem)
823 : {
824 0 : libmesh_assert(f_system);
825 0 : if (f_system->current_local_solution->initialized())
826 : {
827 0 : context = std::make_unique<FEMContext>
828 0 : (*f_system, nullptr,
829 0 : /* allocate local_matrices = */ false);
830 0 : f_fem->init_context(*context);
831 0 : if (g_fem)
832 0 : g_fem->init_context(*context);
833 : }
834 : }
835 :
836 : // There's a chicken-and-egg problem with FEMFunction-based
837 : // Dirichlet constraints: we can't evaluate the FEMFunction
838 : // until we have an initialized local solution vector, we
839 : // can't initialize local solution vectors until we have a
840 : // send list, and we can't generate a send list until we know
841 : // all our constraints
842 : //
843 : // We don't generate constraints on uninitialized systems;
844 : // currently user code will have to reinit() before any
845 : // FEMFunction-based constraints will be correct. This should
846 : // be fine, since user code would want to reinit() after
847 : // setting initial conditions anyway.
848 13501 : if (f_system && context.get())
849 0 : context->pre_fe_reinit(*f_system, elem);
850 :
851 : // Also pre-init the fem_context() we were passed on the current Elem.
852 13501 : fem_context.pre_fe_reinit(fem_context.get_system(), elem);
853 :
854 : // Get a reference to the DOF indices for the current element
855 : const std::vector<dof_id_type> & dof_indices =
856 1127 : fem_context.get_dof_indices(var);
857 :
858 : // The number of DOFs on the element
859 : const unsigned int n_dofs =
860 2254 : cast_int<unsigned int>(dof_indices.size());
861 :
862 : // Fixed vs. free DoFs on edge/face projections
863 14628 : std::vector<char> dof_is_fixed(n_dofs, false); // bools
864 14628 : std::vector<int> free_dof(n_dofs, 0);
865 :
866 : // Zero the interpolated values
867 12374 : Ue.resize (n_dofs); Ue.zero();
868 :
869 : // In general, we need a series of
870 : // projections to ensure a unique and continuous
871 : // solution. We start by interpolating boundary nodes, then
872 : // hold those fixed and project boundary edges, then hold
873 : // those fixed and project boundary faces,
874 :
875 : // Interpolate node values first.
876 : //
877 : // If we have non-vertex nodes that have a boundary nodeset, we
878 : // need to interpolate them directly, but that had better be
879 : // happening in the apply_lagrange_dirichlet_impl code path,
880 : // because you *can't* interpolate to evaluate non-nodal FE
881 : // coefficients.
882 1127 : unsigned int current_dof = 0;
883 117811 : for (unsigned int n=0; n!= sebi.n_nodes; ++n)
884 : {
885 : // FIXME: this should go through the DofMap,
886 : // not duplicate dof_indices code badly!
887 :
888 : // Get the number of DOFs at this node, accounting for
889 : // Elem::p_level() internally.
890 8704 : const unsigned int nc =
891 104310 : FEInterface::n_dofs_at_node (fe_type, elem, n);
892 :
893 : // Get a reference to the "is_boundary_node" flags for the
894 : // current DirichletBoundary object. In case the map does not
895 : // contain an entry for this DirichletBoundary object, it
896 : // means there are no boundary nodes active.
897 104310 : auto is_boundary_node_it = sebi.is_boundary_node_map.find(&dirichlet);
898 :
899 : // The current n is not a boundary node if either there is no
900 : // boundary_node_map for this DirichletBoundary object, or if
901 : // there is but the entry in the corresponding vector is
902 : // false.
903 8704 : const bool not_boundary_node =
904 113014 : (is_boundary_node_it == sebi.is_boundary_node_map.end() ||
905 17408 : !is_boundary_node_it->second[n]);
906 :
907 104310 : if (!elem->is_vertex(n) || not_boundary_node)
908 : {
909 72844 : current_dof += nc;
910 6078 : continue;
911 : }
912 31466 : if (cont == DISCONTINUOUS)
913 : {
914 0 : libmesh_assert_equal_to (nc, 0);
915 : }
916 : // Assume that C_ZERO elements have a single nodal
917 : // value shape function
918 31466 : else if ((cont == C_ZERO) || (fe_type.family == SUBDIVISION))
919 : {
920 1726 : libmesh_assert_equal_to (nc, n_vec_dim);
921 53136 : for (unsigned int c = 0; c < n_vec_dim; c++)
922 : {
923 35126 : Ue(current_dof+c) =
924 32424 : f_component(f, f_fem, context.get(), var_component+c,
925 32424 : elem->point(n), time);
926 35126 : dof_is_fixed[current_dof+c] = true;
927 : }
928 20712 : current_dof += n_vec_dim;
929 17260 : }
930 : // The hermite element vertex shape functions are weird
931 10754 : else if (fe_type.family == HERMITE)
932 : {
933 4272 : Ue(current_dof) =
934 4272 : f_component(f, f_fem, context.get(), var_component,
935 4272 : elem->point(n), time);
936 4272 : dof_is_fixed[current_dof] = true;
937 4272 : current_dof++;
938 : Gradient grad =
939 4272 : g_component(g, g_fem, context.get(), var_component,
940 4272 : elem->point(n), time);
941 : // x derivative
942 4272 : Ue(current_dof) = grad(0);
943 4272 : dof_is_fixed[current_dof] = true;
944 4272 : current_dof++;
945 4272 : if (dim > 1)
946 : {
947 : // We'll finite difference mixed derivatives
948 3744 : Point nxminus = elem->point(n),
949 3744 : nxplus = elem->point(n);
950 3744 : nxminus(0) -= TOLERANCE;
951 3744 : nxplus(0) += TOLERANCE;
952 : Gradient gxminus =
953 3744 : g_component(g, g_fem, context.get(), var_component,
954 3744 : nxminus, time);
955 : Gradient gxplus =
956 3744 : g_component(g, g_fem, context.get(), var_component,
957 3744 : nxplus, time);
958 : // y derivative
959 3744 : Ue(current_dof) = grad(1);
960 3744 : dof_is_fixed[current_dof] = true;
961 3744 : current_dof++;
962 : // xy derivative
963 4056 : Ue(current_dof) = (gxplus(1) - gxminus(1))
964 3432 : / 2. / TOLERANCE;
965 3744 : dof_is_fixed[current_dof] = true;
966 3744 : current_dof++;
967 :
968 3744 : if (dim > 2)
969 : {
970 : // z derivative
971 0 : Ue(current_dof) = grad(2);
972 0 : dof_is_fixed[current_dof] = true;
973 0 : current_dof++;
974 : // xz derivative
975 0 : Ue(current_dof) = (gxplus(2) - gxminus(2))
976 0 : / 2. / TOLERANCE;
977 0 : dof_is_fixed[current_dof] = true;
978 0 : current_dof++;
979 : // We need new points for yz
980 0 : Point nyminus = elem->point(n),
981 0 : nyplus = elem->point(n);
982 0 : nyminus(1) -= TOLERANCE;
983 0 : nyplus(1) += TOLERANCE;
984 : Gradient gyminus =
985 0 : g_component(g, g_fem, context.get(), var_component,
986 0 : nyminus, time);
987 : Gradient gyplus =
988 0 : g_component(g, g_fem, context.get(), var_component,
989 0 : nyplus, time);
990 : // xz derivative
991 0 : Ue(current_dof) = (gyplus(2) - gyminus(2))
992 0 : / 2. / TOLERANCE;
993 0 : dof_is_fixed[current_dof] = true;
994 0 : current_dof++;
995 : // Getting a 2nd order xyz is more tedious
996 0 : Point nxmym = elem->point(n),
997 0 : nxmyp = elem->point(n),
998 0 : nxpym = elem->point(n),
999 0 : nxpyp = elem->point(n);
1000 0 : nxmym(0) -= TOLERANCE;
1001 0 : nxmym(1) -= TOLERANCE;
1002 0 : nxmyp(0) -= TOLERANCE;
1003 0 : nxmyp(1) += TOLERANCE;
1004 0 : nxpym(0) += TOLERANCE;
1005 0 : nxpym(1) -= TOLERANCE;
1006 0 : nxpyp(0) += TOLERANCE;
1007 0 : nxpyp(1) += TOLERANCE;
1008 : Gradient gxmym =
1009 0 : g_component(g, g_fem, context.get(), var_component,
1010 0 : nxmym, time);
1011 : Gradient gxmyp =
1012 0 : g_component(g, g_fem, context.get(), var_component,
1013 0 : nxmyp, time);
1014 : Gradient gxpym =
1015 0 : g_component(g, g_fem, context.get(), var_component,
1016 0 : nxpym, time);
1017 : Gradient gxpyp =
1018 0 : g_component(g, g_fem, context.get(), var_component,
1019 0 : nxpyp, time);
1020 0 : Number gxzplus = (gxpyp(2) - gxmyp(2))
1021 0 : / 2. / TOLERANCE;
1022 0 : Number gxzminus = (gxpym(2) - gxmym(2))
1023 0 : / 2. / TOLERANCE;
1024 : // xyz derivative
1025 0 : Ue(current_dof) = (gxzplus - gxzminus)
1026 0 : / 2. / TOLERANCE;
1027 0 : dof_is_fixed[current_dof] = true;
1028 0 : current_dof++;
1029 : }
1030 : }
1031 : }
1032 : // Assume that other C_ONE elements have a single nodal
1033 : // value shape function and nodal gradient component
1034 : // shape functions
1035 6482 : else if (cont == C_ONE)
1036 : {
1037 544 : libmesh_assert_equal_to (nc, 1 + dim);
1038 7026 : Ue(current_dof) =
1039 6482 : f_component(f, f_fem, context.get(), var_component,
1040 6482 : elem->point(n), time);
1041 6482 : dof_is_fixed[current_dof] = true;
1042 6482 : current_dof++;
1043 : Gradient grad =
1044 6482 : g_component(g, g_fem, context.get(), var_component,
1045 6482 : elem->point(n), time);
1046 19446 : for (unsigned int i=0; i!= dim; ++i)
1047 : {
1048 12964 : Ue(current_dof) = grad(i);
1049 12964 : dof_is_fixed[current_dof] = true;
1050 12964 : current_dof++;
1051 : }
1052 : }
1053 : else
1054 0 : libmesh_error_msg("Unknown continuity cont = " << cont);
1055 : } // end for (n=0..n_nodes)
1056 :
1057 : // In 3D, project any edge values next
1058 13501 : if (dim > 2 && cont != DISCONTINUOUS)
1059 : {
1060 : // Get a pointer to the 1 dimensional (edge) FE for the current
1061 : // var which is stored in the fem_context. This will only be
1062 : // different from side_fe in 3D.
1063 98 : FEGenericBase<OutputType> * edge_fe = nullptr;
1064 98 : fem_context.get_edge_fe(var, edge_fe);
1065 :
1066 : // Set tolerance on underlying FEMap object. This will allow us to
1067 : // avoid spurious negative Jacobian errors while imposing BCs by
1068 : // simply ignoring them. This should only be required in certain
1069 : // special cases, see the DirichletBoundaries comments on this
1070 : // parameter for more information.
1071 1176 : edge_fe->get_fe_map().set_jacobian_tolerance(dirichlet.jacobian_tolerance);
1072 :
1073 : // Pre-request FE data
1074 98 : const std::vector<std::vector<OutputShape>> & phi = edge_fe->get_phi();
1075 1176 : const std::vector<Point> & xyz_values = edge_fe->get_xyz();
1076 294 : const std::vector<Real> & JxW = edge_fe->get_JxW();
1077 :
1078 : // Only pre-request gradients for C1 projections
1079 98 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1080 1176 : if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
1081 : {
1082 0 : const std::vector<std::vector<OutputGradient>> & ref_dphi = edge_fe->get_dphi();
1083 0 : dphi = &ref_dphi;
1084 : }
1085 :
1086 : // Vector to hold edge local DOF indices
1087 196 : std::vector<unsigned int> edge_dofs;
1088 :
1089 : // Get a reference to the "is_boundary_edge" flags for the
1090 : // current DirichletBoundary object. In case the map does not
1091 : // contain an entry for this DirichletBoundary object, it
1092 : // means there are no boundary edges active.
1093 1176 : if (const auto is_boundary_edge_it = sebi.is_boundary_edge_map.find(&dirichlet);
1094 98 : is_boundary_edge_it != sebi.is_boundary_edge_map.end())
1095 : {
1096 15288 : for (unsigned int e=0; e != sebi.n_edges; ++e)
1097 : {
1098 15288 : if (!is_boundary_edge_it->second[e])
1099 14112 : continue;
1100 :
1101 6480 : FEInterface::dofs_on_edge(elem, dim, fe_type, e,
1102 : edge_dofs);
1103 :
1104 : const unsigned int n_edge_dofs =
1105 1080 : cast_int<unsigned int>(edge_dofs.size());
1106 :
1107 : // Some edge dofs are on nodes and already
1108 : // fixed, others are free to calculate
1109 540 : unsigned int free_dofs = 0;
1110 45360 : for (unsigned int i=0; i != n_edge_dofs; ++i)
1111 45360 : if (!dof_is_fixed[edge_dofs[i]])
1112 0 : free_dof[free_dofs++] = i;
1113 :
1114 : // There may be nothing to project
1115 6480 : if (!free_dofs)
1116 5940 : continue;
1117 :
1118 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1119 0 : Fe.resize (free_dofs); Fe.zero();
1120 : // The new edge coefficients
1121 0 : DenseVector<Number> Uedge(free_dofs);
1122 :
1123 : // Initialize FE data on the edge
1124 0 : edge_fe->edge_reinit(elem, e);
1125 0 : const unsigned int n_qp = fem_context.get_edge_qrule().n_points();
1126 :
1127 : // Loop over the quadrature points
1128 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1129 : {
1130 : // solution at the quadrature point
1131 0 : OutputNumber fineval(0);
1132 0 : libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
1133 :
1134 0 : for (unsigned int c = 0; c < n_vec_dim; c++)
1135 0 : f_accessor(c) =
1136 0 : f_component(f, f_fem, context.get(), var_component+c,
1137 0 : xyz_values[qp], time);
1138 :
1139 : // solution grad at the quadrature point
1140 0 : OutputNumberGradient finegrad;
1141 0 : libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
1142 :
1143 : unsigned int g_rank;
1144 0 : switch( FEInterface::field_type( fe_type ) )
1145 : {
1146 0 : case TYPE_SCALAR:
1147 : {
1148 0 : g_rank = 1;
1149 0 : break;
1150 : }
1151 0 : case TYPE_VECTOR:
1152 : {
1153 0 : g_rank = 2;
1154 0 : break;
1155 : }
1156 0 : default:
1157 0 : libmesh_error_msg("Unknown field type!");
1158 : }
1159 :
1160 0 : if (cont == C_ONE)
1161 0 : for (unsigned int c = 0; c < n_vec_dim; c++)
1162 0 : for (unsigned int d = 0; d < g_rank; d++)
1163 0 : g_accessor(c + d*dim ) =
1164 0 : g_component(g, g_fem, context.get(), var_component,
1165 0 : xyz_values[qp], time)(c);
1166 :
1167 : // Form edge projection matrix
1168 0 : for (unsigned int sidei=0, freei=0; sidei != n_edge_dofs; ++sidei)
1169 : {
1170 0 : unsigned int i = edge_dofs[sidei];
1171 : // fixed DoFs aren't test functions
1172 0 : if (dof_is_fixed[i])
1173 0 : continue;
1174 0 : for (unsigned int sidej=0, freej=0; sidej != n_edge_dofs; ++sidej)
1175 : {
1176 0 : unsigned int j = edge_dofs[sidej];
1177 0 : if (dof_is_fixed[j])
1178 0 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
1179 0 : JxW[qp] * Ue(j);
1180 : else
1181 0 : Ke(freei,freej) += phi[i][qp] *
1182 0 : phi[j][qp] * JxW[qp];
1183 0 : if (cont == C_ONE)
1184 : {
1185 0 : if (dof_is_fixed[j])
1186 0 : Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp]) ) *
1187 0 : JxW[qp] * Ue(j);
1188 : else
1189 0 : Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1190 0 : * JxW[qp];
1191 : }
1192 0 : if (!dof_is_fixed[j])
1193 0 : freej++;
1194 : }
1195 0 : Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1196 0 : if (cont == C_ONE)
1197 0 : Fe(freei) += (finegrad.contract( (*dphi)[i][qp]) ) *
1198 : JxW[qp];
1199 0 : freei++;
1200 : }
1201 : }
1202 :
1203 0 : Ke.cholesky_solve(Fe, Uedge);
1204 :
1205 : // Transfer new edge solutions to element
1206 0 : for (unsigned int i=0; i != free_dofs; ++i)
1207 : {
1208 0 : Number & ui = Ue(edge_dofs[free_dof[i]]);
1209 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1210 : std::abs(ui - Uedge(i)) < TOLERANCE);
1211 0 : ui = Uedge(i);
1212 0 : dof_is_fixed[edge_dofs[free_dof[i]]] = true;
1213 : }
1214 : } // end for (e = 0..n_edges)
1215 : } // end if (is_boundary_edge_it != sebi.is_boundary_edge_map.end())
1216 : } // end if (dim > 2 && cont != DISCONTINUOUS)
1217 :
1218 : // Project any side values (edges in 2D, faces in 3D)
1219 13501 : if (dim > 1 && cont != DISCONTINUOUS)
1220 : {
1221 11582 : FEGenericBase<OutputType> * side_fe = nullptr;
1222 11582 : fem_context.get_side_fe(var, side_fe);
1223 :
1224 : // Set tolerance on underlying FEMap object. This will allow us to
1225 : // avoid spurious negative Jacobian errors while imposing BCs by
1226 : // simply ignoring them. This should only be required in certain
1227 : // special cases, see the DirichletBoundaries comments on this
1228 : // parameter for more information.
1229 12637 : side_fe->get_fe_map().set_jacobian_tolerance(dirichlet.jacobian_tolerance);
1230 :
1231 : // Pre-request FE data
1232 1055 : const std::vector<std::vector<OutputShape>> & phi = side_fe->get_phi();
1233 12637 : const std::vector<Point> & xyz_values = side_fe->get_xyz();
1234 3165 : const std::vector<Real> & JxW = side_fe->get_JxW();
1235 :
1236 : // Only pre-request gradients for C1 projections
1237 1055 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1238 12637 : if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
1239 : {
1240 415 : const std::vector<std::vector<OutputGradient>> & ref_dphi = side_fe->get_dphi();
1241 415 : dphi = &ref_dphi;
1242 : }
1243 :
1244 : // Vector to hold side local DOF indices
1245 2110 : std::vector<unsigned int> side_dofs;
1246 :
1247 : // Get a reference to the "is_boundary_side" flags for the
1248 : // current DirichletBoundary object. In case the map does not
1249 : // contain an entry for this DirichletBoundary object, it
1250 : // means there are no boundary sides active.
1251 12637 : if (const auto is_boundary_side_it = sebi.is_boundary_side_map.find(&dirichlet);
1252 1055 : is_boundary_side_it != sebi.is_boundary_side_map.end())
1253 : {
1254 60484 : for (unsigned int s=0; s != sebi.n_sides; ++s)
1255 : {
1256 52230 : if (!is_boundary_side_it->second[s])
1257 37010 : continue;
1258 :
1259 15061 : FEInterface::dofs_on_side(elem, dim, fe_type, s,
1260 : side_dofs);
1261 :
1262 : const unsigned int n_side_dofs =
1263 2514 : cast_int<unsigned int>(side_dofs.size());
1264 :
1265 : // Some side dofs are on nodes/edges and already
1266 : // fixed, others are free to calculate
1267 1257 : unsigned int free_dofs = 0;
1268 101144 : for (unsigned int i=0; i != n_side_dofs; ++i)
1269 100457 : if (!dof_is_fixed[side_dofs[i]])
1270 13926 : free_dof[free_dofs++] = i;
1271 :
1272 : // There may be nothing to project
1273 15061 : if (!free_dofs)
1274 3542 : continue;
1275 :
1276 10262 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1277 10262 : Fe.resize (free_dofs); Fe.zero();
1278 : // The new side coefficients
1279 11197 : DenseVector<Number> Uside(free_dofs);
1280 :
1281 : // Initialize FE data on the side
1282 11197 : side_fe->reinit(elem, s);
1283 935 : const unsigned int n_qp = fem_context.get_side_qrule().n_points();
1284 :
1285 : // Loop over the quadrature points
1286 49625 : for (unsigned int qp=0; qp<n_qp; qp++)
1287 : {
1288 : // solution at the quadrature point
1289 38428 : OutputNumber fineval(0);
1290 3210 : libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
1291 :
1292 76856 : for (unsigned int c = 0; c < n_vec_dim; c++)
1293 38428 : f_accessor(c) =
1294 38428 : f_component(f, f_fem, context.get(), var_component+c,
1295 38428 : xyz_values[qp], time);
1296 :
1297 : // solution grad at the quadrature point
1298 3210 : OutputNumberGradient finegrad;
1299 3210 : libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
1300 :
1301 : unsigned int g_rank;
1302 38428 : switch( FEInterface::field_type( fe_type ) )
1303 : {
1304 3210 : case TYPE_SCALAR:
1305 : {
1306 3210 : g_rank = 1;
1307 3210 : break;
1308 : }
1309 0 : case TYPE_VECTOR:
1310 : {
1311 0 : g_rank = 2;
1312 0 : break;
1313 : }
1314 0 : default:
1315 0 : libmesh_error_msg("Unknown field type!");
1316 : }
1317 :
1318 38428 : if (cont == C_ONE)
1319 25448 : for (unsigned int c = 0; c < n_vec_dim; c++)
1320 25448 : for (unsigned int d = 0; d < g_rank; d++)
1321 13792 : g_accessor(c + d*dim ) =
1322 13792 : g_component(g, g_fem, context.get(), var_component,
1323 13792 : xyz_values[qp], time)(c);
1324 :
1325 : // Form side projection matrix
1326 212432 : for (unsigned int sidei=0, freei=0; sidei != n_side_dofs; ++sidei)
1327 : {
1328 174004 : unsigned int i = side_dofs[sidei];
1329 : // fixed DoFs aren't test functions
1330 188558 : if (dof_is_fixed[i])
1331 117060 : continue;
1332 258416 : for (unsigned int sidej=0, freej=0; sidej != n_side_dofs; ++sidej)
1333 : {
1334 212164 : unsigned int j = side_dofs[sidej];
1335 229898 : if (dof_is_fixed[j])
1336 227372 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
1337 131404 : JxW[qp] * Ue(j);
1338 : else
1339 91716 : Ke(freei,freej) += phi[i][qp] *
1340 80240 : phi[j][qp] * JxW[qp];
1341 212164 : if (cont == C_ONE)
1342 : {
1343 96544 : if (dof_is_fixed[j])
1344 114792 : Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1345 69936 : JxW[qp] * Ue(j);
1346 : else
1347 18064 : Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1348 14860 : * JxW[qp];
1349 : }
1350 229898 : if (!dof_is_fixed[j])
1351 68764 : freej++;
1352 : }
1353 57838 : Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1354 46252 : if (cont == C_ONE)
1355 14860 : Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1356 : JxW[qp];
1357 46252 : freei++;
1358 : }
1359 : }
1360 :
1361 11197 : Ke.cholesky_solve(Fe, Uside);
1362 :
1363 : // Transfer new side solutions to element
1364 24050 : for (unsigned int i=0; i != free_dofs; ++i)
1365 : {
1366 14999 : Number & ui = Ue(side_dofs[free_dof[i]]);
1367 :
1368 1073 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1369 : std::abs(ui - Uside(i)) < TOLERANCE);
1370 12853 : ui = Uside(i);
1371 :
1372 13926 : dof_is_fixed[side_dofs[free_dof[i]]] = true;
1373 : }
1374 : } // end for (s = 0..n_sides)
1375 : } // end if (is_boundary_side_it != sebi.is_boundary_side_map.end())
1376 : } // end if (dim > 1 && cont != DISCONTINUOUS)
1377 :
1378 : // Project any shellface values
1379 13501 : if (dim == 2 && cont != DISCONTINUOUS)
1380 : {
1381 957 : FEGenericBase<OutputType> * fe = nullptr;
1382 957 : fem_context.get_element_fe(var, fe, dim);
1383 :
1384 : // Set tolerance on underlying FEMap object. This will allow us to
1385 : // avoid spurious negative Jacobian errors while imposing BCs by
1386 : // simply ignoring them. This should only be required in certain
1387 : // special cases, see the DirichletBoundaries comments on this
1388 : // parameter for more information.
1389 11461 : fe->get_fe_map().set_jacobian_tolerance(dirichlet.jacobian_tolerance);
1390 :
1391 : // Pre-request FE data
1392 957 : const std::vector<std::vector<OutputShape>> & phi = fe->get_phi();
1393 11461 : const std::vector<Point> & xyz_values = fe->get_xyz();
1394 2871 : const std::vector<Real> & JxW = fe->get_JxW();
1395 :
1396 : // Only pre-request gradients for C1 projections
1397 957 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1398 11461 : if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
1399 : {
1400 415 : const std::vector<std::vector<OutputGradient>> & ref_dphi = fe->get_dphi();
1401 415 : dphi = &ref_dphi;
1402 : }
1403 :
1404 : // Get a reference to the "is_boundary_shellface" flags for the
1405 : // current DirichletBoundary object. In case the map does not
1406 : // contain an entry for this DirichletBoundary object, it
1407 : // means there are no boundary shellfaces active.
1408 11461 : if (const auto is_boundary_shellface_it = sebi.is_boundary_shellface_map.find(&dirichlet);
1409 957 : is_boundary_shellface_it != sebi.is_boundary_shellface_map.end())
1410 : {
1411 0 : for (unsigned int shellface=0; shellface != 2; ++shellface)
1412 : {
1413 0 : if (!is_boundary_shellface_it->second[shellface])
1414 0 : continue;
1415 :
1416 : // A shellface has the same dof indices as the element itself
1417 0 : std::vector<unsigned int> shellface_dofs(n_dofs);
1418 0 : std::iota(shellface_dofs.begin(), shellface_dofs.end(), 0);
1419 :
1420 : // Some shellface dofs are on nodes/edges and already
1421 : // fixed, others are free to calculate
1422 0 : unsigned int free_dofs = 0;
1423 0 : for (unsigned int i=0; i != n_dofs; ++i)
1424 0 : if (!dof_is_fixed[shellface_dofs[i]])
1425 0 : free_dof[free_dofs++] = i;
1426 :
1427 : // There may be nothing to project
1428 0 : if (!free_dofs)
1429 0 : continue;
1430 :
1431 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1432 0 : Fe.resize (free_dofs); Fe.zero();
1433 : // The new shellface coefficients
1434 0 : DenseVector<Number> Ushellface(free_dofs);
1435 :
1436 : // Initialize FE data on the element
1437 0 : fe->reinit (elem);
1438 0 : const unsigned int n_qp = fem_context.get_element_qrule().n_points();
1439 :
1440 : // Loop over the quadrature points
1441 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1442 : {
1443 : // solution at the quadrature point
1444 0 : OutputNumber fineval(0);
1445 0 : libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
1446 :
1447 0 : for (unsigned int c = 0; c < n_vec_dim; c++)
1448 0 : f_accessor(c) =
1449 0 : f_component(f, f_fem, context.get(), var_component+c,
1450 0 : xyz_values[qp], time);
1451 :
1452 : // solution grad at the quadrature point
1453 0 : OutputNumberGradient finegrad;
1454 0 : libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
1455 :
1456 : unsigned int g_rank;
1457 0 : switch( FEInterface::field_type( fe_type ) )
1458 : {
1459 0 : case TYPE_SCALAR:
1460 : {
1461 0 : g_rank = 1;
1462 0 : break;
1463 : }
1464 0 : case TYPE_VECTOR:
1465 : {
1466 0 : g_rank = 2;
1467 0 : break;
1468 : }
1469 0 : default:
1470 0 : libmesh_error_msg("Unknown field type!");
1471 : }
1472 :
1473 0 : if (cont == C_ONE)
1474 0 : for (unsigned int c = 0; c < n_vec_dim; c++)
1475 0 : for (unsigned int d = 0; d < g_rank; d++)
1476 0 : g_accessor(c + d*dim ) =
1477 0 : g_component(g, g_fem, context.get(), var_component,
1478 0 : xyz_values[qp], time)(c);
1479 :
1480 : // Form shellface projection matrix
1481 0 : for (unsigned int shellfacei=0, freei=0;
1482 0 : shellfacei != n_dofs; ++shellfacei)
1483 : {
1484 0 : unsigned int i = shellface_dofs[shellfacei];
1485 : // fixed DoFs aren't test functions
1486 0 : if (dof_is_fixed[i])
1487 0 : continue;
1488 0 : for (unsigned int shellfacej=0, freej=0;
1489 0 : shellfacej != n_dofs; ++shellfacej)
1490 : {
1491 0 : unsigned int j = shellface_dofs[shellfacej];
1492 0 : if (dof_is_fixed[j])
1493 0 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
1494 0 : JxW[qp] * Ue(j);
1495 : else
1496 0 : Ke(freei,freej) += phi[i][qp] *
1497 0 : phi[j][qp] * JxW[qp];
1498 0 : if (cont == C_ONE)
1499 : {
1500 0 : if (dof_is_fixed[j])
1501 0 : Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1502 0 : JxW[qp] * Ue(j);
1503 : else
1504 0 : Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1505 0 : * JxW[qp];
1506 : }
1507 0 : if (!dof_is_fixed[j])
1508 0 : freej++;
1509 : }
1510 0 : Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1511 0 : if (cont == C_ONE)
1512 0 : Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1513 : JxW[qp];
1514 0 : freei++;
1515 : }
1516 : }
1517 :
1518 0 : Ke.cholesky_solve(Fe, Ushellface);
1519 :
1520 : // Transfer new shellface solutions to element
1521 0 : for (unsigned int i=0; i != free_dofs; ++i)
1522 : {
1523 0 : Number & ui = Ue(shellface_dofs[free_dof[i]]);
1524 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1525 : std::abs(ui - Ushellface(i)) < TOLERANCE);
1526 0 : ui = Ushellface(i);
1527 0 : dof_is_fixed[shellface_dofs[free_dof[i]]] = true;
1528 : }
1529 : } // end for (shellface = 0..2)
1530 : } // end if (is_boundary_shellface_it != sebi.is_boundary_shellface_map.end())
1531 : } // end if (dim == 2 && cont != DISCONTINUOUS)
1532 :
1533 : // Lock the DofConstraints since it is shared among threads.
1534 : {
1535 2254 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1536 :
1537 177221 : for (unsigned int i = 0; i < n_dofs; i++)
1538 : {
1539 27332 : DofConstraintRow empty_row;
1540 177386 : if (dof_is_fixed[i] && !libmesh_isnan(Ue(i)))
1541 87498 : add_fn (dof_indices[i], empty_row, Ue(i));
1542 : }
1543 : }
1544 22494 : } // apply_dirichlet_impl
1545 :
1546 : public:
1547 574 : ConstrainDirichlet (DofMap & dof_map_in,
1548 : const MeshBase & mesh_in,
1549 : const Real time_in,
1550 : const DirichletBoundaries & dirichlets_in,
1551 20314 : const AddConstraint & add_in) :
1552 19166 : dof_map(dof_map_in),
1553 19166 : mesh(mesh_in),
1554 19166 : time(time_in),
1555 19166 : dirichlets(dirichlets_in),
1556 20314 : add_fn(add_in) { }
1557 :
1558 : // This class can be default copy/move constructed.
1559 : ConstrainDirichlet (ConstrainDirichlet &&) = default;
1560 : ConstrainDirichlet (const ConstrainDirichlet &) = default;
1561 :
1562 : // This class cannot be default copy/move assigned because it
1563 : // contains reference members.
1564 : ConstrainDirichlet & operator= (const ConstrainDirichlet &) = delete;
1565 : ConstrainDirichlet & operator= (ConstrainDirichlet &&) = delete;
1566 :
1567 20380 : void operator()(const ConstElemRange & range) const
1568 : {
1569 : /**
1570 : * This method examines an arbitrary boundary solution to calculate
1571 : * corresponding Dirichlet constraints on the current mesh. The
1572 : * input function \p f gives the arbitrary solution.
1573 : */
1574 :
1575 : // Figure out which System the DirichletBoundary objects are
1576 : // defined for. We break out of the loop as soon as we encounter a
1577 : // valid System pointer, the assumption is thus that all Variables
1578 : // are defined on the same System.
1579 596 : System * system = nullptr;
1580 :
1581 : // Map from boundary_id -> set<pair<id,DirichletBoundary*>> objects which
1582 : // are active on that boundary_id. Later we will use this to determine
1583 : // which DirichletBoundary objects to loop over for each Elem.
1584 : std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>>
1585 1192 : boundary_id_to_ordered_dirichlet_boundaries;
1586 :
1587 48374 : for (auto dirichlet_id : index_range(dirichlets))
1588 : {
1589 : // Pointer to the current DirichletBoundary object
1590 27994 : const auto & dirichlet = dirichlets[dirichlet_id];
1591 :
1592 : // Construct mapping from boundary_id -> (dirichlet_id, DirichletBoundary)
1593 85575 : for (const auto & b_id : dirichlet->b)
1594 57581 : boundary_id_to_ordered_dirichlet_boundaries[b_id].emplace(dirichlet_id, dirichlet.get());
1595 :
1596 65047 : for (const auto & var : dirichlet->variables)
1597 : {
1598 37053 : const Variable & variable = dof_map.variable(var);
1599 2156 : auto current_system = variable.system();
1600 :
1601 37053 : if (!system)
1602 596 : system = current_system;
1603 : else
1604 16673 : libmesh_error_msg_if(current_system != system,
1605 : "All variables should be defined on the same System");
1606 : }
1607 : }
1608 :
1609 : // If we found no System, it could be because none of the
1610 : // Variables have one defined, or because there are
1611 : // DirichletBoundary objects with no Variables defined on
1612 : // them. These situations both indicate a likely error in the
1613 : // setup of a problem, so let's throw an error in this case.
1614 20380 : libmesh_error_msg_if(!system, "Valid System not found for any Variables.");
1615 :
1616 : // Construct a FEMContext object for the System on which the
1617 : // Variables in our DirichletBoundary objects are defined. This
1618 : // will be used in the apply_dirichlet_impl() function.
1619 : // We're not going to use elem_jacobian or subjacobians here so
1620 : // don't allocate them.
1621 : auto fem_context = std::make_unique<FEMContext>
1622 20976 : (*system, nullptr, /* allocate local_matrices = */ false);
1623 :
1624 : // At the time we are using this FEMContext, the current_local_solution
1625 : // vector is not initialized, but also we don't need it, so set
1626 : // the algebraic_type flag to DOFS_ONLY.
1627 596 : fem_context->set_algebraic_type(FEMContext::DOFS_ONLY);
1628 :
1629 : // Boundary info for the current mesh
1630 20380 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1631 :
1632 : // This object keeps track of the BoundaryInfo for a single Elem
1633 20976 : SingleElemBoundaryInfo sebi(boundary_info, boundary_id_to_ordered_dirichlet_boundaries);
1634 :
1635 : // Iterate over all the elements in the range
1636 1568315 : for (const auto & elem : range)
1637 : {
1638 : // We only calculate Dirichlet constraints on active
1639 : // elements
1640 1547935 : if (!elem->active())
1641 344533 : continue;
1642 :
1643 : // Reinitialize BoundaryInfo data structures for the current elem
1644 1171243 : bool has_dirichlet_constraint = sebi.reinit(elem);
1645 :
1646 : // If this Elem has no boundary ids, go to the next one.
1647 1171243 : if (!has_dirichlet_constraint)
1648 943713 : continue;
1649 :
1650 285950 : for (const auto & db_pair : sebi.ordered_dbs)
1651 : {
1652 : // Get pointer to the DirichletBoundary object
1653 12555 : const auto & dirichlet = db_pair.second;
1654 :
1655 : // Loop over all the variables which this DirichletBoundary object is responsible for
1656 326801 : for (const auto & var : dirichlet->variables)
1657 : {
1658 179735 : const Variable & variable = dof_map.variable(var);
1659 :
1660 : // Make sure that the Variable and the DofMap agree on
1661 : // what number this variable is.
1662 15391 : libmesh_assert_equal_to(variable.number(), var);
1663 :
1664 15391 : const FEType & fe_type = variable.type();
1665 :
1666 179735 : if (fe_type.family == SCALAR)
1667 0 : continue;
1668 :
1669 179735 : switch( FEInterface::field_type( fe_type ) )
1670 : {
1671 178559 : case TYPE_SCALAR:
1672 : {
1673 : // For Lagrange FEs we don't need to do a full
1674 : // blown projection, we can just interpolate
1675 : // values directly.
1676 178559 : if (fe_type.family == LAGRANGE)
1677 166234 : this->apply_lagrange_dirichlet_impl<Real>( sebi, variable, *dirichlet, *fem_context );
1678 : else
1679 12325 : this->apply_dirichlet_impl<Real>( sebi, variable, *dirichlet, *fem_context );
1680 15293 : break;
1681 : }
1682 98 : case TYPE_VECTOR:
1683 : {
1684 1176 : this->apply_dirichlet_impl<RealGradient>( sebi, variable, *dirichlet, *fem_context );
1685 98 : break;
1686 : }
1687 0 : default:
1688 0 : libmesh_error_msg("Unknown field type!");
1689 : }
1690 : } // for (var : variables)
1691 : } // for (db_pair : ordered_dbs)
1692 : } // for (elem : range)
1693 39568 : } // operator()
1694 :
1695 : }; // class ConstrainDirichlet
1696 :
1697 :
1698 : #endif // LIBMESH_ENABLE_DIRICHLET
1699 :
1700 :
1701 : } // anonymous namespace
1702 :
1703 :
1704 :
1705 : namespace libMesh
1706 : {
1707 :
1708 : // ------------------------------------------------------------
1709 : // DofMap member functions
1710 :
1711 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1712 :
1713 :
1714 1453455 : dof_id_type DofMap::n_constrained_dofs() const
1715 : {
1716 35458 : parallel_object_only();
1717 :
1718 1453455 : dof_id_type nc_dofs = this->n_local_constrained_dofs();
1719 1453455 : this->comm().sum(nc_dofs);
1720 1453455 : return nc_dofs;
1721 : }
1722 :
1723 :
1724 1473973 : dof_id_type DofMap::n_local_constrained_dofs() const
1725 : {
1726 : const DofConstraints::const_iterator lower =
1727 36038 : _dof_constraints.lower_bound(this->first_dof()),
1728 : upper =
1729 36038 : _dof_constraints.lower_bound(this->end_dof());
1730 :
1731 1510011 : return cast_int<dof_id_type>(std::distance(lower, upper));
1732 : }
1733 :
1734 :
1735 :
1736 335191 : void DofMap::create_dof_constraints(const MeshBase & mesh, Real time)
1737 : {
1738 9608 : parallel_object_only();
1739 :
1740 18056 : LOG_SCOPE("create_dof_constraints()", "DofMap");
1741 :
1742 9608 : libmesh_assert (mesh.is_prepared());
1743 :
1744 : // The user might have set boundary conditions after the mesh was
1745 : // prepared; we should double-check that those boundary conditions
1746 : // are still consistent.
1747 : #ifdef DEBUG
1748 9608 : MeshTools::libmesh_assert_valid_boundary_ids(mesh);
1749 : #endif
1750 :
1751 : // We might get constraint equations from AMR hanging nodes in
1752 : // 2D/3D, or from spline constraint rows or boundary conditions in
1753 : // any dimension
1754 : const bool possible_local_constraints = false
1755 335191 : || !mesh.n_elem()
1756 335191 : || mesh.n_constraint_rows()
1757 : #ifdef LIBMESH_ENABLE_AMR
1758 333971 : || mesh.mesh_dimension() > 1
1759 : #endif
1760 : #ifdef LIBMESH_ENABLE_PERIODIC
1761 46064 : || !_periodic_boundaries->empty()
1762 : #endif
1763 : #ifdef LIBMESH_ENABLE_DIRICHLET
1764 389293 : || !_dirichlet_boundaries->empty()
1765 : #endif
1766 : ;
1767 :
1768 : // Even if we don't have constraints, another processor might.
1769 335191 : bool possible_global_constraints = possible_local_constraints;
1770 : #if defined(LIBMESH_ENABLE_PERIODIC) || defined(LIBMESH_ENABLE_DIRICHLET) || defined(LIBMESH_ENABLE_AMR)
1771 9608 : libmesh_assert(this->comm().verify(mesh.is_serial()));
1772 :
1773 335191 : this->comm().max(possible_global_constraints);
1774 : #endif
1775 :
1776 : // Recalculate dof constraints from scratch. (Or just clear them,
1777 : // if the user has just deleted their last dirichlet/periodic/user
1778 : // constraint)
1779 : // Note: any _stashed_dof_constraints are not cleared as it
1780 : // may be the user's intention to restore them later.
1781 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1782 9608 : _dof_constraints.clear();
1783 9608 : _primal_constraint_values.clear();
1784 9608 : _adjoint_constraint_values.clear();
1785 : #endif
1786 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1787 9608 : _node_constraints.clear();
1788 : #endif
1789 :
1790 335191 : if (!possible_global_constraints)
1791 41169 : return;
1792 :
1793 : // Here we build the hanging node constraints. This is done
1794 : // by enforcing the condition u_a = u_b along hanging sides.
1795 : // u_a = u_b is collocated at the nodes of side a, which gives
1796 : // one row of the constraint matrix.
1797 :
1798 : // Processors only compute their local constraints
1799 588044 : ConstElemRange range (mesh.local_elements_begin(),
1800 882066 : mesh.local_elements_end());
1801 :
1802 : // Global computation fails if we're using a FEMFunctionBase BC on a
1803 : // ReplicatedMesh in parallel
1804 : // ConstElemRange range (mesh.elements_begin(),
1805 : // mesh.elements_end());
1806 :
1807 : // compute_periodic_constraints requires a point_locator() from our
1808 : // Mesh, but point_locator() construction is parallel and threaded.
1809 : // Rather than nest threads within threads we'll make sure it's
1810 : // preconstructed.
1811 : #ifdef LIBMESH_ENABLE_PERIODIC
1812 570701 : bool need_point_locator = !_periodic_boundaries->empty() && !range.empty();
1813 :
1814 294022 : this->comm().max(need_point_locator);
1815 :
1816 294022 : if (need_point_locator)
1817 805 : mesh.sub_point_locator();
1818 : #endif
1819 :
1820 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1821 25344 : Threads::parallel_for (range,
1822 16896 : ComputeNodeConstraints (_node_constraints,
1823 : #ifdef LIBMESH_ENABLE_PERIODIC
1824 8448 : *_periodic_boundaries,
1825 : #endif
1826 : mesh));
1827 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1828 :
1829 :
1830 : // Look at all the variables in the system. Reset the element
1831 : // range at each iteration -- there is no need to reconstruct it.
1832 294022 : const auto n_vars = this->n_variables();
1833 623550 : for (unsigned int variable_number=0; variable_number<n_vars;
1834 320122 : ++variable_number, range.reset())
1835 338934 : Threads::parallel_for (range,
1836 338934 : ComputeConstraints (_dof_constraints,
1837 : *this,
1838 : #ifdef LIBMESH_ENABLE_PERIODIC
1839 9406 : *_periodic_boundaries,
1840 : #endif
1841 : mesh,
1842 : variable_number));
1843 :
1844 : #ifdef LIBMESH_ENABLE_DIRICHLET
1845 :
1846 294022 : if (!_dirichlet_boundaries->empty())
1847 : {
1848 : // Sanity check that the boundary ids associated with the
1849 : // DirichletBoundary objects are actually present in the
1850 : // mesh. We do this check by default, but in cases where you
1851 : // intentionally add "inconsistent but valid" DirichletBoundary
1852 : // objects in parallel, this check can deadlock since it does a
1853 : // collective communication internally. In that case it is
1854 : // possible to disable this check by setting the flag to false.
1855 18129 : if (_verify_dirichlet_bc_consistency)
1856 42740 : for (const auto & dirichlet : *_dirichlet_boundaries)
1857 24611 : this->check_dirichlet_bcid_consistency(mesh, *dirichlet);
1858 :
1859 : // Threaded loop over local over elems applying all Dirichlet BCs
1860 : Threads::parallel_for
1861 18641 : (range,
1862 18129 : ConstrainDirichlet(*this, mesh, time, *_dirichlet_boundaries,
1863 18129 : AddPrimalConstraint(*this)));
1864 :
1865 : // Threaded loop over local over elems per QOI applying all adjoint
1866 : // Dirichlet BCs. Note that the ConstElemRange is reset before each
1867 : // execution of Threads::parallel_for().
1868 :
1869 20314 : for (auto qoi_index : index_range(_adjoint_dirichlet_boundaries))
1870 : {
1871 : const DirichletBoundaries & adb_q =
1872 124 : *(_adjoint_dirichlet_boundaries[qoi_index]);
1873 :
1874 2185 : if (!adb_q.empty())
1875 : Threads::parallel_for
1876 2247 : (range.reset(),
1877 2185 : ConstrainDirichlet(*this, mesh, time, adb_q,
1878 2247 : AddAdjointConstraint(*this, qoi_index)));
1879 : }
1880 : }
1881 :
1882 : #endif // LIBMESH_ENABLE_DIRICHLET
1883 :
1884 : // Handle spline node constraints last, so we can try to move
1885 : // existing constraints onto the spline basis if necessary.
1886 294022 : this->process_mesh_constraint_rows(mesh);
1887 : }
1888 :
1889 :
1890 :
1891 294022 : void DofMap::process_mesh_constraint_rows(const MeshBase & mesh)
1892 : {
1893 : // If we already have simple Dirichlet constraints (with right hand
1894 : // sides but with no coupling between DoFs) on spline-constrained FE
1895 : // nodes, then we'll need a solve to compute the corresponding
1896 : // constraints on the relevant spline nodes. (If we already have
1897 : // constraints with coupling between DoFs on spline-constrained FE
1898 : // nodes, then we'll need to go sit down and cry until we figure out
1899 : // how to handle that.)
1900 :
1901 8448 : const auto & constraint_rows = mesh.get_constraint_rows();
1902 :
1903 : // This routine is too expensive to use unless we really might
1904 : // need it
1905 : #ifdef DEBUG
1906 8448 : bool constraint_rows_empty = constraint_rows.empty();
1907 8448 : this->comm().min(constraint_rows_empty);
1908 8448 : libmesh_assert_equal_to(constraint_rows_empty, !mesh.n_constraint_rows());
1909 : #endif
1910 :
1911 294022 : if (!mesh.n_constraint_rows())
1912 292802 : return;
1913 :
1914 : // We can't handle periodic boundary conditions on spline meshes
1915 : // yet.
1916 : #ifdef LIBMESH_ENABLE_PERIODIC
1917 1220 : libmesh_error_msg_if (!_periodic_boundaries->empty(),
1918 : "Periodic boundary conditions are not yet implemented for spline meshes");
1919 : #endif
1920 :
1921 : // We can handle existing Dirichlet constraints, but we'll need
1922 : // to do solves to project them down onto the spline basis.
1923 1186 : std::unique_ptr<SparsityPattern::Build> sp;
1924 1186 : std::unique_ptr<SparseMatrix<Number>> mat;
1925 :
1926 : const unsigned int n_adjoint_rhs =
1927 1220 : _adjoint_constraint_values.size();
1928 :
1929 : // [0] for primal rhs, [q+1] for adjoint qoi q
1930 : std::vector<std::unique_ptr<NumericVector<Number>>>
1931 1288 : solve_rhs(n_adjoint_rhs+1);
1932 :
1933 : // Keep track of which spline DoFs will be Dirichlet.
1934 : // We use a set here to make it easier to find what processors
1935 : // to send the DoFs to later.
1936 68 : std::set<dof_id_type> my_dirichlet_spline_dofs;
1937 :
1938 : // And keep track of which non-spline Dofs were Dirichlet
1939 68 : std::unordered_set<dof_id_type> was_previously_constrained;
1940 :
1941 68 : const unsigned int sys_num = this->sys_number();
1942 192741 : for (auto & node_row : constraint_rows)
1943 : {
1944 191521 : const Node * node = node_row.first;
1945 16528 : libmesh_assert(node == mesh.node_ptr(node->id()));
1946 :
1947 : // Each processor only computes its own (and in distributed
1948 : // cases, is only guaranteed to have the dependency data to
1949 : // compute its own) constraints here.
1950 208049 : if (node->processor_id() != mesh.processor_id())
1951 83793 : continue;
1952 :
1953 329608 : for (auto var_num : IntRange<unsigned int>(0, this->n_variables()))
1954 : {
1955 19154 : const FEFamily & fe_family = this->variable_type(var_num).family;
1956 :
1957 : // constraint_rows only applies to nodal variables
1958 230144 : if (fe_family != LAGRANGE &&
1959 18788 : fe_family != RATIONAL_BERNSTEIN)
1960 0 : continue;
1961 :
1962 38308 : DofConstraintRow dc_row;
1963 :
1964 : const dof_id_type constrained_id =
1965 230144 : node->dof_number(sys_num, var_num, 0);
1966 729015 : for (const auto & [pr, val] : node_row.second)
1967 : {
1968 498871 : const Elem * spline_elem = pr.first;
1969 41623 : libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1970 :
1971 : const Node & spline_node =
1972 498871 : spline_elem->node_ref(pr.second);
1973 :
1974 : const dof_id_type spline_dof_id =
1975 498871 : spline_node.dof_number(sys_num, var_num, 0);
1976 498871 : dc_row[spline_dof_id] = val;
1977 : }
1978 :
1979 : // See if we already have a constraint here.
1980 230144 : if (this->is_constrained_dof(constrained_id))
1981 : {
1982 396 : was_previously_constrained.insert(constrained_id);
1983 :
1984 : // Keep track of which spline DoFs will be
1985 : // inheriting this non-spline DoF's constraints
1986 11664 : for (auto & row_entry : dc_row)
1987 6912 : my_dirichlet_spline_dofs.insert(row_entry.first);
1988 :
1989 : // If it wasn't a simple Dirichlet-type constraint
1990 : // then I don't know what to do with it. We'll make
1991 : // this an assertion only because this should only
1992 : // crop up with periodic boundary conditions, which
1993 : // we've already made sure we don't have.
1994 396 : libmesh_assert(_dof_constraints[constrained_id].empty());
1995 : }
1996 :
1997 : // Add the constraint, replacing any previous, so we can
1998 : // use the new constraint in setting up the solve below
1999 38308 : this->add_constraint_row(constrained_id, dc_row, false);
2000 : }
2001 : }
2002 :
2003 : // my_dirichlet_spline_dofs may now include DoFs whose owners
2004 : // don't know they need to become spline DoFs! We need to push
2005 : // this data to them.
2006 1220 : if (this->comm().size() > 1)
2007 : {
2008 : std::unordered_map
2009 : <processor_id_type, std::vector<dof_id_type>>
2010 68 : their_dirichlet_spline_dofs;
2011 :
2012 : // If we ever change the underlying container here then we'd
2013 : // better do some kind of sort before using it; we'll rely
2014 : // on sorting to make the processor id lookup efficient.
2015 34 : libmesh_assert(std::is_sorted(my_dirichlet_spline_dofs.begin(),
2016 : my_dirichlet_spline_dofs.end()));
2017 1189 : processor_id_type destination_pid = 0;
2018 3475 : for (auto d : my_dirichlet_spline_dofs)
2019 : {
2020 216 : libmesh_assert_less(d, this->end_dof(this->comm().size()-1));
2021 3378 : while (d >= this->end_dof(destination_pid))
2022 872 : destination_pid++;
2023 :
2024 2502 : if (destination_pid != this->processor_id())
2025 612 : their_dirichlet_spline_dofs[destination_pid].push_back(d);
2026 : }
2027 :
2028 : auto receive_dof_functor =
2029 88 : [& my_dirichlet_spline_dofs]
2030 : (processor_id_type,
2031 92 : const std::vector<dof_id_type> & dofs)
2032 : {
2033 96 : my_dirichlet_spline_dofs.insert(dofs.begin(), dofs.end());
2034 1193 : };
2035 :
2036 : Parallel::push_parallel_vector_data
2037 1189 : (this->comm(), their_dirichlet_spline_dofs, receive_dof_functor);
2038 : }
2039 :
2040 :
2041 : // If anyone had any prior constraints in effect, then we need
2042 : // to convert them to constraints on the spline nodes.
2043 : //
2044 : // NOT simply testing prior_constraints here; maybe it turned
2045 : // out that all our constraints were on non-spline-constrained
2046 : // parts of a hybrid mesh?
2047 : bool important_prior_constraints =
2048 1220 : !was_previously_constrained.empty();
2049 1220 : this->comm().max(important_prior_constraints);
2050 :
2051 1220 : if (important_prior_constraints)
2052 : {
2053 : // Now that we have the spline constraints added, we can
2054 : // finally construct a sparsity pattern that correctly
2055 : // accounts for those constraints!
2056 552 : mat = SparseMatrix<Number>::build(this->comm());
2057 568 : for (auto q : IntRange<unsigned int>(0, n_adjoint_rhs+1))
2058 : {
2059 284 : solve_rhs[q] = NumericVector<Number>::build(this->comm());
2060 300 : solve_rhs[q]->init(this->n_dofs(), this->n_local_dofs(),
2061 16 : false, PARALLEL);
2062 : }
2063 :
2064 : // We need to compute our own sparsity pattern, to take into
2065 : // account the particularly non-sparse rows that can be
2066 : // created by the spline constraints we just added.
2067 284 : mat->attach_dof_map(*this);
2068 552 : sp = this->build_sparsity(mesh);
2069 284 : mat->attach_sparsity_pattern(*sp);
2070 284 : mat->init();
2071 :
2072 92892 : for (auto & node_row : constraint_rows)
2073 : {
2074 92608 : const Node * node = node_row.first;
2075 8712 : libmesh_assert(node == mesh.node_ptr(node->id()));
2076 :
2077 370432 : for (auto var_num : IntRange<unsigned int>(0, this->n_variables()))
2078 : {
2079 26136 : const FEFamily & fe_family = this->variable_type(var_num).family;
2080 :
2081 : // constraint_rows only applies to nodal variables
2082 277824 : if (fe_family != LAGRANGE &&
2083 26136 : fe_family != RATIONAL_BERNSTEIN)
2084 0 : continue;
2085 :
2086 : const dof_id_type constrained_id =
2087 277824 : node->dof_number(sys_num, var_num, 0);
2088 :
2089 52272 : if (was_previously_constrained.count(constrained_id))
2090 : {
2091 9504 : for (auto q : IntRange<int>(0, n_adjoint_rhs+1))
2092 : {
2093 5544 : DenseMatrix<Number> K(1,1);
2094 4752 : DenseVector<Number> F(1);
2095 5148 : std::vector<dof_id_type> dof_indices(1, constrained_id);
2096 :
2097 4356 : K(0,0) = 1;
2098 :
2099 4752 : DofConstraintValueMap & vals = q ?
2100 4356 : _adjoint_constraint_values[q-1] :
2101 396 : _primal_constraint_values;
2102 :
2103 : DofConstraintValueMap::const_iterator rhsit =
2104 396 : vals.find(constrained_id);
2105 4752 : F(0) = (rhsit == vals.end()) ? 0 : rhsit->second;
2106 :
2107 : // We no longer need any rhs values here directly.
2108 4752 : if (rhsit != vals.end())
2109 4356 : vals.erase(rhsit);
2110 :
2111 : this->heterogeneously_constrain_element_matrix_and_vector
2112 4752 : (K, F, dof_indices, false, q ? (q-1) : -1);
2113 4752 : if (!q)
2114 4752 : mat->add_matrix(K, dof_indices);
2115 4752 : solve_rhs[q]->add_vector(F, dof_indices);
2116 3960 : }
2117 : }
2118 : }
2119 : }
2120 :
2121 : // Any DoFs that aren't part of any constraint, directly or
2122 : // indirectly, need a diagonal term to make the matrix
2123 : // here invertible.
2124 17224 : for (dof_id_type d : IntRange<dof_id_type>(this->first_dof(),
2125 220712 : this->end_dof()))
2126 50472 : if (!was_previously_constrained.count(d) &&
2127 16560 : !my_dirichlet_spline_dofs.count(d))
2128 196128 : mat->add(d,d,1);
2129 :
2130 : // At this point, we're finally ready to solve for Dirichlet
2131 : // constraint values on spline nodes.
2132 : std::unique_ptr<LinearSolver<Number>> linear_solver =
2133 292 : LinearSolver<Number>::build(this->comm());
2134 :
2135 : std::unique_ptr<NumericVector<Number>> projected_vals =
2136 292 : NumericVector<Number>::build(this->comm());
2137 :
2138 292 : projected_vals->init(this->n_dofs(), this->n_local_dofs(),
2139 16 : false, PARALLEL);
2140 :
2141 16 : DofConstraintRow empty_row;
2142 3488 : for (auto sd : my_dirichlet_spline_dofs)
2143 2832 : if (this->local_index(sd))
2144 216 : this->add_constraint_row(sd, empty_row);
2145 :
2146 568 : for (auto q : IntRange<unsigned int>(0, n_adjoint_rhs+1))
2147 : {
2148 : // FIXME: we don't have an EquationSystems here, but I'd
2149 : // rather not hardcode these...
2150 8 : const double tol = double(TOLERANCE * TOLERANCE);
2151 8 : const unsigned int max_its = 5000;
2152 :
2153 284 : linear_solver->solve(*mat, *projected_vals,
2154 308 : *(solve_rhs[q]), tol, max_its);
2155 :
2156 284 : DofConstraintValueMap & vals = q ?
2157 276 : _adjoint_constraint_values[q-1] :
2158 8 : _primal_constraint_values;
2159 :
2160 3488 : for (auto sd : my_dirichlet_spline_dofs)
2161 2832 : if (this->local_index(sd))
2162 : {
2163 2592 : Number constraint_rhs = (*projected_vals)(sd);
2164 :
2165 : std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
2166 216 : vals.emplace(sd, constraint_rhs);
2167 2592 : if (!rhs_it.second)
2168 2592 : rhs_it.first->second = constraint_rhs;
2169 : }
2170 : }
2171 268 : }
2172 1152 : }
2173 :
2174 :
2175 :
2176 593929 : void DofMap::add_constraint_row (const dof_id_type dof_number,
2177 : const DofConstraintRow & constraint_row,
2178 : const Number constraint_rhs,
2179 : const bool forbid_constraint_overwrite)
2180 : {
2181 : // Optionally allow the user to overwrite constraints. Defaults to false.
2182 593929 : libmesh_error_msg_if(forbid_constraint_overwrite && this->is_constrained_dof(dof_number),
2183 : "ERROR: DOF " << dof_number << " was already constrained!");
2184 :
2185 48175 : libmesh_assert_less(dof_number, this->n_dofs());
2186 :
2187 : // There is an implied "1" on the diagonal of the constraint row, and the user
2188 : // should not try to manually set _any_ value on the diagonal.
2189 48175 : libmesh_assert_msg(!constraint_row.count(dof_number),
2190 : "Error: constraint_row for dof " << dof_number <<
2191 : " should not contain an entry for dof " << dof_number);
2192 :
2193 : #ifndef NDEBUG
2194 93300 : for (const auto & pr : constraint_row)
2195 45125 : libmesh_assert_less(pr.first, this->n_dofs());
2196 : #endif
2197 :
2198 : // Store the constraint_row in the map
2199 593929 : _dof_constraints.insert_or_assign(dof_number, constraint_row);
2200 :
2201 : std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
2202 48175 : _primal_constraint_values.emplace(dof_number, constraint_rhs);
2203 593929 : if (!rhs_it.second)
2204 4752 : rhs_it.first->second = constraint_rhs;
2205 593929 : }
2206 :
2207 :
2208 14304 : void DofMap::add_adjoint_constraint_row (const unsigned int qoi_index,
2209 : const dof_id_type dof_number,
2210 : const DofConstraintRow & /*constraint_row*/,
2211 : const Number constraint_rhs,
2212 : const bool forbid_constraint_overwrite)
2213 : {
2214 : // Optionally allow the user to overwrite constraints. Defaults to false.
2215 14304 : if (forbid_constraint_overwrite)
2216 : {
2217 14304 : libmesh_error_msg_if(!this->is_constrained_dof(dof_number),
2218 : "ERROR: DOF " << dof_number << " has no corresponding primal constraint!");
2219 : #ifndef NDEBUG
2220 : // No way to do this without a non-normalized tolerance?
2221 :
2222 : // // If the user passed in more than just the rhs, let's check the
2223 : // // coefficients for consistency
2224 : // if (!constraint_row.empty())
2225 : // {
2226 : // DofConstraintRow row = _dof_constraints[dof_number];
2227 : // for (const auto & [dof, val] : row)
2228 : // libmesh_assert(constraint_row.find(dof)->second == val);
2229 : // }
2230 : //
2231 : // if (_adjoint_constraint_values[qoi_index].find(dof_number) !=
2232 : // _adjoint_constraint_values[qoi_index].end())
2233 : // libmesh_assert_equal_to(_adjoint_constraint_values[qoi_index][dof_number],
2234 : // constraint_rhs);
2235 :
2236 : #endif
2237 : }
2238 :
2239 : // Creates the map of rhs values if it doesn't already exist; then
2240 : // adds the current value to that map
2241 :
2242 : // Store the rhs value in the map
2243 14304 : _adjoint_constraint_values[qoi_index].insert_or_assign(dof_number, constraint_rhs);
2244 14304 : }
2245 :
2246 :
2247 :
2248 :
2249 0 : void DofMap::print_dof_constraints(std::ostream & os,
2250 : bool print_nonlocal) const
2251 : {
2252 0 : parallel_object_only();
2253 :
2254 : std::string local_constraints =
2255 0 : this->get_local_constraints(print_nonlocal);
2256 :
2257 0 : if (this->processor_id())
2258 : {
2259 0 : this->comm().send(0, local_constraints);
2260 : }
2261 : else
2262 : {
2263 0 : os << "Processor 0:\n";
2264 0 : os << local_constraints;
2265 :
2266 0 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
2267 : {
2268 0 : this->comm().receive(p, local_constraints);
2269 0 : os << "Processor " << p << ":\n";
2270 0 : os << local_constraints;
2271 : }
2272 : }
2273 0 : }
2274 :
2275 0 : std::string DofMap::get_local_constraints(bool print_nonlocal) const
2276 : {
2277 0 : std::ostringstream os;
2278 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2279 0 : if (print_nonlocal)
2280 0 : os << "All ";
2281 : else
2282 0 : os << "Local ";
2283 :
2284 0 : os << "Node Constraints:"
2285 0 : << std::endl;
2286 :
2287 0 : for (const auto & [node, pr] : _node_constraints)
2288 : {
2289 : // Skip non-local nodes if requested
2290 0 : if (!print_nonlocal &&
2291 0 : node->processor_id() != this->processor_id())
2292 0 : continue;
2293 :
2294 0 : const NodeConstraintRow & row = pr.first;
2295 0 : const Point & offset = pr.second;
2296 :
2297 0 : os << "Constraints for Node id " << node->id()
2298 0 : << ": \t";
2299 :
2300 0 : for (const auto & [cnode, val] : row)
2301 0 : os << " (" << cnode->id() << "," << val << ")\t";
2302 :
2303 0 : os << "rhs: " << offset;
2304 :
2305 0 : os << std::endl;
2306 : }
2307 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2308 :
2309 0 : if (print_nonlocal)
2310 0 : os << "All ";
2311 : else
2312 0 : os << "Local ";
2313 :
2314 0 : os << "DoF Constraints:"
2315 0 : << std::endl;
2316 :
2317 0 : for (const auto & [i, row] : _dof_constraints)
2318 : {
2319 : // Skip non-local dofs if requested
2320 0 : if (!print_nonlocal && !this->local_index(i))
2321 0 : continue;
2322 :
2323 : DofConstraintValueMap::const_iterator rhsit =
2324 0 : _primal_constraint_values.find(i);
2325 0 : const Number rhs = (rhsit == _primal_constraint_values.end()) ?
2326 0 : 0 : rhsit->second;
2327 :
2328 0 : os << "Constraints for DoF " << i
2329 0 : << ": \t";
2330 :
2331 0 : for (const auto & item : row)
2332 0 : os << " (" << item.first << "," << item.second << ")\t";
2333 :
2334 0 : os << "rhs: " << rhs;
2335 0 : os << std::endl;
2336 : }
2337 :
2338 0 : for (unsigned int qoi_index = 0,
2339 0 : n_qois = cast_int<unsigned int>(_adjoint_dirichlet_boundaries.size());
2340 0 : qoi_index != n_qois; ++qoi_index)
2341 : {
2342 0 : os << "Adjoint " << qoi_index << " DoF rhs values:"
2343 0 : << std::endl;
2344 :
2345 0 : if (auto adjoint_map_it = _adjoint_constraint_values.find(qoi_index);
2346 0 : adjoint_map_it != _adjoint_constraint_values.end())
2347 0 : for (const auto & [i, rhs] : adjoint_map_it->second)
2348 : {
2349 : // Skip non-local dofs if requested
2350 0 : if (!print_nonlocal && !this->local_index(i))
2351 0 : continue;
2352 :
2353 0 : os << "RHS for DoF " << i
2354 0 : << ": " << rhs;
2355 :
2356 0 : os << std::endl;
2357 : }
2358 : }
2359 :
2360 0 : return os.str();
2361 0 : }
2362 :
2363 :
2364 :
2365 27029125 : void DofMap::constrain_element_matrix (DenseMatrix<Number> & matrix,
2366 : std::vector<dof_id_type> & elem_dofs,
2367 : bool asymmetric_constraint_rows) const
2368 : {
2369 1838984 : libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2370 1838984 : libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2371 :
2372 : // check for easy return
2373 27029125 : if (this->_dof_constraints.empty())
2374 11593719 : return;
2375 :
2376 : // The constrained matrix is built up as C^T K C.
2377 19080142 : DenseMatrix<Number> C;
2378 :
2379 :
2380 15435406 : this->build_constraint_matrix (C, elem_dofs);
2381 :
2382 5467104 : LOG_SCOPE("constrain_elem_matrix()", "DofMap");
2383 :
2384 : // It is possible that the matrix is not constrained at all.
2385 15454938 : if ((C.m() == matrix.m()) &&
2386 221615 : (C.n() == elem_dofs.size())) // It the matrix is constrained
2387 : {
2388 : // Compute the matrix-matrix-matrix product C^T K C
2389 221615 : matrix.left_multiply_transpose (C);
2390 221615 : matrix.right_multiply (C);
2391 :
2392 :
2393 19532 : libmesh_assert_equal_to (matrix.m(), matrix.n());
2394 19532 : libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2395 19532 : libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2396 :
2397 :
2398 1805915 : for (unsigned int i=0,
2399 39064 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2400 1844979 : i != n_elem_dofs; i++)
2401 : // If the DOF is constrained
2402 1767466 : if (this->is_constrained_dof(elem_dofs[i]))
2403 : {
2404 6981808 : for (auto j : make_range(matrix.n()))
2405 6365431 : matrix(i,j) = 0.;
2406 :
2407 583153 : matrix(i,i) = 1.;
2408 :
2409 586265 : if (asymmetric_constraint_rows)
2410 : {
2411 : DofConstraints::const_iterator
2412 3684 : pos = _dof_constraints.find(elem_dofs[i]);
2413 :
2414 3684 : libmesh_assert (pos != _dof_constraints.end());
2415 :
2416 3684 : const DofConstraintRow & constraint_row = pos->second;
2417 :
2418 : // This is an overzealous assertion in the presence of
2419 : // heterogeneous constraints: we now can constrain "u_i = c"
2420 : // with no other u_j terms involved.
2421 : //
2422 : // libmesh_assert (!constraint_row.empty());
2423 :
2424 104892 : for (const auto & item : constraint_row)
2425 350886 : for (unsigned int j=0; j != n_elem_dofs; j++)
2426 319140 : if (elem_dofs[j] == item.first)
2427 53670 : matrix(i,j) = -item.second;
2428 : }
2429 : }
2430 : } // end if is constrained...
2431 11790670 : }
2432 :
2433 :
2434 :
2435 14699830 : void DofMap::constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
2436 : DenseVector<Number> & rhs,
2437 : std::vector<dof_id_type> & elem_dofs,
2438 : bool asymmetric_constraint_rows) const
2439 : {
2440 1301536 : libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2441 1301536 : libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2442 1301536 : libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2443 :
2444 : // check for easy return
2445 14699830 : if (this->_dof_constraints.empty())
2446 7942431 : return;
2447 :
2448 : // The constrained matrix is built up as C^T K C.
2449 : // The constrained RHS is built up as C^T F
2450 8176901 : DenseMatrix<Number> C;
2451 :
2452 6757399 : this->build_constraint_matrix (C, elem_dofs);
2453 :
2454 2129253 : LOG_SCOPE("cnstrn_elem_mat_vec()", "DofMap");
2455 :
2456 : // It is possible that the matrix is not constrained at all.
2457 6878833 : if ((C.m() == matrix.m()) &&
2458 1239199 : (C.n() == elem_dofs.size())) // It the matrix is constrained
2459 : {
2460 : // Compute the matrix-matrix-matrix product C^T K C
2461 1239199 : matrix.left_multiply_transpose (C);
2462 1239199 : matrix.right_multiply (C);
2463 :
2464 :
2465 121434 : libmesh_assert_equal_to (matrix.m(), matrix.n());
2466 121434 : libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2467 121434 : libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2468 :
2469 :
2470 17863628 : for (unsigned int i=0,
2471 242868 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2472 18106496 : i != n_elem_dofs; i++)
2473 18447018 : if (this->is_constrained_dof(elem_dofs[i]))
2474 : {
2475 219793792 : for (auto j : make_range(matrix.n()))
2476 202244889 : matrix(i,j) = 0.;
2477 :
2478 : // If the DOF is constrained
2479 5610516 : matrix(i,i) = 1.;
2480 :
2481 : // This will put a nonsymmetric entry in the constraint
2482 : // row to ensure that the linear system produces the
2483 : // correct value for the constrained DOF.
2484 5930073 : if (asymmetric_constraint_rows)
2485 : {
2486 : DofConstraints::const_iterator
2487 179995 : pos = _dof_constraints.find(elem_dofs[i]);
2488 :
2489 179995 : libmesh_assert (pos != _dof_constraints.end());
2490 :
2491 179995 : const DofConstraintRow & constraint_row = pos->second;
2492 :
2493 : // p refinement creates empty constraint rows
2494 : // libmesh_assert (!constraint_row.empty());
2495 :
2496 4948016 : for (const auto & item : constraint_row)
2497 50902690 : for (unsigned int j=0; j != n_elem_dofs; j++)
2498 52327335 : if (elem_dofs[j] == item.first)
2499 3402021 : matrix(i,j) = -item.second;
2500 : }
2501 : }
2502 :
2503 :
2504 : // Compute the matrix-vector product C^T F
2505 242868 : DenseVector<Number> old_rhs(rhs);
2506 :
2507 : // compute matrix/vector product
2508 1239199 : C.vector_mult_transpose(rhs, old_rhs);
2509 : } // end if is constrained...
2510 5337897 : }
2511 :
2512 :
2513 :
2514 684225 : void DofMap::heterogeneously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
2515 : DenseVector<Number> & rhs,
2516 : std::vector<dof_id_type> & elem_dofs,
2517 : bool asymmetric_constraint_rows,
2518 : int qoi_index) const
2519 : {
2520 67251 : libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2521 67251 : libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2522 67251 : libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2523 :
2524 : // check for easy return
2525 684225 : if (this->_dof_constraints.empty())
2526 30122 : return;
2527 :
2528 : // The constrained matrix is built up as C^T K C.
2529 : // The constrained RHS is built up as C^T (F - K H)
2530 787351 : DenseMatrix<Number> C;
2531 654103 : DenseVector<Number> H;
2532 :
2533 654103 : this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
2534 :
2535 199872 : LOG_SCOPE("hetero_cnstrn_elem_mat_vec()", "DofMap");
2536 :
2537 : // It is possible that the matrix is not constrained at all.
2538 670393 : if ((C.m() == matrix.m()) &&
2539 172593 : (C.n() == elem_dofs.size())) // It the matrix is constrained
2540 : {
2541 : // We may have rhs values to use later
2542 16290 : const DofConstraintValueMap * rhs_values = nullptr;
2543 172593 : if (qoi_index < 0)
2544 172593 : rhs_values = &_primal_constraint_values;
2545 0 : else if (auto it = _adjoint_constraint_values.find(qoi_index);
2546 0 : it != _adjoint_constraint_values.end())
2547 0 : rhs_values = &it->second;
2548 :
2549 : // Compute matrix/vector product K H
2550 172593 : DenseVector<Number> KH;
2551 172593 : matrix.vector_mult(KH, H);
2552 :
2553 : // Compute the matrix-vector product C^T (F - KH)
2554 32580 : DenseVector<Number> F_minus_KH(rhs);
2555 156303 : F_minus_KH -= KH;
2556 172593 : C.vector_mult_transpose(rhs, F_minus_KH);
2557 :
2558 : // Compute the matrix-matrix-matrix product C^T K C
2559 172593 : matrix.left_multiply_transpose (C);
2560 172593 : matrix.right_multiply (C);
2561 :
2562 16290 : libmesh_assert_equal_to (matrix.m(), matrix.n());
2563 16290 : libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2564 16290 : libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2565 :
2566 3154276 : for (unsigned int i=0,
2567 32580 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2568 3186856 : i != n_elem_dofs; i++)
2569 : {
2570 3308973 : const dof_id_type dof_id = elem_dofs[i];
2571 :
2572 2719553 : if (this->is_constrained_dof(dof_id))
2573 : {
2574 19116891 : for (auto j : make_range(matrix.n()))
2575 16936421 : matrix(i,j) = 0.;
2576 :
2577 : // If the DOF is constrained
2578 838223 : matrix(i,i) = 1.;
2579 :
2580 : // This will put a nonsymmetric entry in the constraint
2581 : // row to ensure that the linear system produces the
2582 : // correct value for the constrained DOF.
2583 902874 : if (asymmetric_constraint_rows)
2584 : {
2585 : DofConstraints::const_iterator
2586 85547 : pos = _dof_constraints.find(dof_id);
2587 :
2588 85547 : libmesh_assert (pos != _dof_constraints.end());
2589 :
2590 85547 : const DofConstraintRow & constraint_row = pos->second;
2591 :
2592 969202 : for (const auto & item : constraint_row)
2593 1675760 : for (unsigned int j=0; j != n_elem_dofs; j++)
2594 1719300 : if (elem_dofs[j] == item.first)
2595 95250 : matrix(i,j) = -item.second;
2596 :
2597 881231 : if (rhs_values)
2598 : {
2599 : const DofConstraintValueMap::const_iterator valpos =
2600 85547 : rhs_values->find(dof_id);
2601 :
2602 904168 : rhs(i) = (valpos == rhs_values->end()) ?
2603 22937 : 0 : valpos->second;
2604 : }
2605 : }
2606 : else
2607 19883 : rhs(i) = 0.;
2608 : }
2609 : }
2610 :
2611 : } // end if is constrained...
2612 520855 : }
2613 :
2614 :
2615 1760 : void DofMap::heterogeneously_constrain_element_jacobian_and_residual
2616 : (DenseMatrix<Number> & matrix,
2617 : DenseVector<Number> & rhs,
2618 : std::vector<dof_id_type> & elem_dofs,
2619 : NumericVector<Number> & solution_local) const
2620 : {
2621 160 : libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2622 160 : libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2623 160 : libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2624 :
2625 160 : libmesh_assert (solution_local.type() == SERIAL ||
2626 : solution_local.type() == GHOSTED);
2627 :
2628 : // check for easy return
2629 1760 : if (this->_dof_constraints.empty())
2630 0 : return;
2631 :
2632 : // The constrained matrix is built up as C^T K C.
2633 : // The constrained RHS is built up as C^T F
2634 : // Asymmetric residual terms are added if we do not have x = Cx+h
2635 1920 : DenseMatrix<Number> C;
2636 1600 : DenseVector<Number> H;
2637 :
2638 1760 : this->build_constraint_matrix_and_vector (C, H, elem_dofs);
2639 :
2640 320 : LOG_SCOPE("hetero_cnstrn_elem_jac_res()", "DofMap");
2641 :
2642 : // It is possible that the matrix is not constrained at all.
2643 1920 : if ((C.m() != matrix.m()) ||
2644 1760 : (C.n() != elem_dofs.size()))
2645 0 : return;
2646 :
2647 : // Compute the matrix-vector product C^T F
2648 320 : DenseVector<Number> old_rhs(rhs);
2649 1760 : C.vector_mult_transpose(rhs, old_rhs);
2650 :
2651 : // Compute the matrix-matrix-matrix product C^T K C
2652 1760 : matrix.left_multiply_transpose (C);
2653 1760 : matrix.right_multiply (C);
2654 :
2655 160 : libmesh_assert_equal_to (matrix.m(), matrix.n());
2656 160 : libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2657 160 : libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2658 :
2659 8480 : for (unsigned int i=0,
2660 320 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2661 8800 : i != n_elem_dofs; i++)
2662 : {
2663 7680 : const dof_id_type dof_id = elem_dofs[i];
2664 :
2665 7040 : if (auto pos = _dof_constraints.find(dof_id);
2666 640 : pos != _dof_constraints.end())
2667 : {
2668 26400 : for (auto j : make_range(matrix.n()))
2669 21120 : matrix(i,j) = 0.;
2670 :
2671 : // If the DOF is constrained
2672 5280 : matrix(i,i) = 1.;
2673 :
2674 : // This will put a nonsymmetric entry in the constraint
2675 : // row to ensure that the linear system produces the
2676 : // correct value for the constrained DOF.
2677 480 : const DofConstraintRow & constraint_row = pos->second;
2678 :
2679 5280 : for (const auto & item : constraint_row)
2680 0 : for (unsigned int j=0; j != n_elem_dofs; j++)
2681 0 : if (elem_dofs[j] == item.first)
2682 0 : matrix(i,j) = -item.second;
2683 :
2684 : const DofConstraintValueMap::const_iterator valpos =
2685 480 : _primal_constraint_values.find(dof_id);
2686 :
2687 480 : Number & rhs_val = rhs(i);
2688 5760 : rhs_val = (valpos == _primal_constraint_values.end()) ?
2689 5280 : 0 : -valpos->second;
2690 5280 : for (const auto & [constraining_dof, coef] : constraint_row)
2691 0 : rhs_val -= coef * solution_local(constraining_dof);
2692 5280 : rhs_val += solution_local(dof_id);
2693 : }
2694 : }
2695 1440 : }
2696 :
2697 :
2698 1760 : void DofMap::heterogeneously_constrain_element_residual
2699 : (DenseVector<Number> & rhs,
2700 : std::vector<dof_id_type> & elem_dofs,
2701 : NumericVector<Number> & solution_local) const
2702 : {
2703 160 : libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2704 :
2705 160 : libmesh_assert (solution_local.type() == SERIAL ||
2706 : solution_local.type() == GHOSTED);
2707 :
2708 : // check for easy return
2709 1760 : if (this->_dof_constraints.empty())
2710 0 : return;
2711 :
2712 : // The constrained RHS is built up as C^T F
2713 : // Asymmetric residual terms are added if we do not have x = Cx+h
2714 1920 : DenseMatrix<Number> C;
2715 1600 : DenseVector<Number> H;
2716 :
2717 1760 : this->build_constraint_matrix_and_vector (C, H, elem_dofs);
2718 :
2719 320 : LOG_SCOPE("hetero_cnstrn_elem_res()", "DofMap");
2720 :
2721 : // It is possible that the element is not constrained at all.
2722 2080 : if ((C.m() != rhs.size()) ||
2723 1760 : (C.n() != elem_dofs.size()))
2724 0 : return;
2725 :
2726 : // Compute the matrix-vector product C^T F
2727 320 : DenseVector<Number> old_rhs(rhs);
2728 1760 : C.vector_mult_transpose(rhs, old_rhs);
2729 :
2730 8480 : for (unsigned int i=0,
2731 320 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2732 8800 : i != n_elem_dofs; i++)
2733 : {
2734 7680 : const dof_id_type dof_id = elem_dofs[i];
2735 :
2736 7040 : if (auto pos = _dof_constraints.find(dof_id);
2737 640 : pos != _dof_constraints.end())
2738 : {
2739 : // This will put a nonsymmetric entry in the constraint
2740 : // row to ensure that the linear system produces the
2741 : // correct value for the constrained DOF.
2742 480 : const DofConstraintRow & constraint_row = pos->second;
2743 :
2744 : const DofConstraintValueMap::const_iterator valpos =
2745 480 : _primal_constraint_values.find(dof_id);
2746 :
2747 480 : Number & rhs_val = rhs(i);
2748 10080 : rhs_val = (valpos == _primal_constraint_values.end()) ?
2749 5280 : 0 : -valpos->second;
2750 5280 : for (const auto & [constraining_dof, coef] : constraint_row)
2751 0 : rhs_val -= coef * solution_local(constraining_dof);
2752 5280 : rhs_val += solution_local(dof_id);
2753 : }
2754 : }
2755 1440 : }
2756 :
2757 :
2758 0 : void DofMap::constrain_element_residual
2759 : (DenseVector<Number> & rhs,
2760 : std::vector<dof_id_type> & elem_dofs,
2761 : NumericVector<Number> & solution_local) const
2762 : {
2763 0 : libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2764 :
2765 0 : libmesh_assert (solution_local.type() == SERIAL ||
2766 : solution_local.type() == GHOSTED);
2767 :
2768 : // check for easy return
2769 0 : if (this->_dof_constraints.empty())
2770 0 : return;
2771 :
2772 : // The constrained RHS is built up as C^T F
2773 0 : DenseMatrix<Number> C;
2774 :
2775 0 : this->build_constraint_matrix (C, elem_dofs);
2776 :
2777 0 : LOG_SCOPE("cnstrn_elem_residual()", "DofMap");
2778 :
2779 : // It is possible that the matrix is not constrained at all.
2780 0 : if (C.n() != elem_dofs.size())
2781 0 : return;
2782 :
2783 : // Compute the matrix-vector product C^T F
2784 0 : DenseVector<Number> old_rhs(rhs);
2785 0 : C.vector_mult_transpose(rhs, old_rhs);
2786 :
2787 0 : for (unsigned int i=0,
2788 0 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2789 0 : i != n_elem_dofs; i++)
2790 : {
2791 0 : const dof_id_type dof_id = elem_dofs[i];
2792 :
2793 0 : if (auto pos = _dof_constraints.find(dof_id);
2794 0 : pos != _dof_constraints.end())
2795 : {
2796 : // This will put a nonsymmetric entry in the constraint
2797 : // row to ensure that the linear system produces the
2798 : // correct value for the constrained DOF.
2799 0 : const DofConstraintRow & constraint_row = pos->second;
2800 :
2801 0 : Number & rhs_val = rhs(i);
2802 0 : rhs_val = 0;
2803 0 : for (const auto & [constraining_dof, coef] : constraint_row)
2804 0 : rhs_val -= coef * solution_local(constraining_dof);
2805 0 : rhs_val += solution_local(dof_id);
2806 : }
2807 : }
2808 0 : }
2809 :
2810 :
2811 1680 : void DofMap::heterogeneously_constrain_element_vector (const DenseMatrix<Number> & matrix,
2812 : DenseVector<Number> & rhs,
2813 : std::vector<dof_id_type> & elem_dofs,
2814 : bool asymmetric_constraint_rows,
2815 : int qoi_index) const
2816 : {
2817 140 : libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2818 140 : libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2819 140 : libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2820 :
2821 : // check for easy return
2822 1680 : if (this->_dof_constraints.empty())
2823 0 : return;
2824 :
2825 : // The constrained matrix is built up as C^T K C.
2826 : // The constrained RHS is built up as C^T (F - K H)
2827 1960 : DenseMatrix<Number> C;
2828 1680 : DenseVector<Number> H;
2829 :
2830 1680 : this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
2831 :
2832 420 : LOG_SCOPE("hetero_cnstrn_elem_vec()", "DofMap");
2833 :
2834 : // It is possible that the matrix is not constrained at all.
2835 1820 : if ((C.m() == matrix.m()) &&
2836 1680 : (C.n() == elem_dofs.size())) // It the matrix is constrained
2837 : {
2838 : // We may have rhs values to use later
2839 140 : const DofConstraintValueMap * rhs_values = nullptr;
2840 1680 : if (qoi_index < 0)
2841 0 : rhs_values = &_primal_constraint_values;
2842 1820 : else if (auto it = _adjoint_constraint_values.find(qoi_index);
2843 140 : it != _adjoint_constraint_values.end())
2844 1680 : rhs_values = &it->second;
2845 :
2846 : // Compute matrix/vector product K H
2847 1680 : DenseVector<Number> KH;
2848 1680 : matrix.vector_mult(KH, H);
2849 :
2850 : // Compute the matrix-vector product C^T (F - KH)
2851 280 : DenseVector<Number> F_minus_KH(rhs);
2852 1540 : F_minus_KH -= KH;
2853 1680 : C.vector_mult_transpose(rhs, F_minus_KH);
2854 :
2855 16520 : for (unsigned int i=0,
2856 280 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2857 16800 : i != n_elem_dofs; i++)
2858 : {
2859 16380 : const dof_id_type dof_id = elem_dofs[i];
2860 :
2861 9420 : if (this->is_constrained_dof(dof_id))
2862 : {
2863 : // This will put a nonsymmetric entry in the constraint
2864 : // row to ensure that the linear system produces the
2865 : // correct value for the constrained DOF.
2866 5328 : if (asymmetric_constraint_rows && rhs_values)
2867 : {
2868 : const DofConstraintValueMap::const_iterator valpos =
2869 0 : rhs_values->find(dof_id);
2870 :
2871 0 : rhs(i) = (valpos == rhs_values->end()) ?
2872 0 : 0 : valpos->second;
2873 : }
2874 : else
2875 4884 : rhs(i) = 0.;
2876 : }
2877 : }
2878 :
2879 : } // end if is constrained...
2880 1400 : }
2881 :
2882 :
2883 :
2884 :
2885 0 : void DofMap::constrain_element_matrix (DenseMatrix<Number> & matrix,
2886 : std::vector<dof_id_type> & row_dofs,
2887 : std::vector<dof_id_type> & col_dofs,
2888 : bool asymmetric_constraint_rows) const
2889 : {
2890 0 : libmesh_assert_equal_to (row_dofs.size(), matrix.m());
2891 0 : libmesh_assert_equal_to (col_dofs.size(), matrix.n());
2892 :
2893 : // check for easy return
2894 0 : if (this->_dof_constraints.empty())
2895 0 : return;
2896 :
2897 : // The constrained matrix is built up as R^T K C.
2898 0 : DenseMatrix<Number> R;
2899 0 : DenseMatrix<Number> C;
2900 :
2901 : // Safeguard against the user passing us the same
2902 : // object for row_dofs and col_dofs. If that is done
2903 : // the calls to build_matrix would fail
2904 0 : std::vector<dof_id_type> orig_row_dofs(row_dofs);
2905 0 : std::vector<dof_id_type> orig_col_dofs(col_dofs);
2906 :
2907 0 : this->build_constraint_matrix (R, orig_row_dofs);
2908 0 : this->build_constraint_matrix (C, orig_col_dofs);
2909 :
2910 0 : LOG_SCOPE("constrain_elem_matrix()", "DofMap");
2911 :
2912 0 : row_dofs = orig_row_dofs;
2913 0 : col_dofs = orig_col_dofs;
2914 :
2915 0 : bool constraint_found = false;
2916 :
2917 : // K_constrained = R^T K C
2918 :
2919 0 : if ((R.m() == matrix.m()) &&
2920 0 : (R.n() == row_dofs.size()))
2921 : {
2922 0 : matrix.left_multiply_transpose (R);
2923 0 : constraint_found = true;
2924 : }
2925 :
2926 0 : if ((C.m() == matrix.n()) &&
2927 0 : (C.n() == col_dofs.size()))
2928 : {
2929 0 : matrix.right_multiply (C);
2930 0 : constraint_found = true;
2931 : }
2932 :
2933 : // It is possible that the matrix is not constrained at all.
2934 0 : if (constraint_found)
2935 : {
2936 0 : libmesh_assert_equal_to (matrix.m(), row_dofs.size());
2937 0 : libmesh_assert_equal_to (matrix.n(), col_dofs.size());
2938 :
2939 :
2940 0 : for (unsigned int i=0,
2941 0 : n_row_dofs = cast_int<unsigned int>(row_dofs.size());
2942 0 : i != n_row_dofs; i++)
2943 0 : if (this->is_constrained_dof(row_dofs[i]))
2944 : {
2945 0 : for (auto j : make_range(matrix.n()))
2946 : {
2947 0 : if (row_dofs[i] != col_dofs[j])
2948 0 : matrix(i,j) = 0.;
2949 : else // If the DOF is constrained
2950 0 : matrix(i,j) = 1.;
2951 : }
2952 :
2953 0 : if (asymmetric_constraint_rows)
2954 : {
2955 : DofConstraints::const_iterator
2956 0 : pos = _dof_constraints.find(row_dofs[i]);
2957 :
2958 0 : libmesh_assert (pos != _dof_constraints.end());
2959 :
2960 0 : const DofConstraintRow & constraint_row = pos->second;
2961 :
2962 0 : libmesh_assert (!constraint_row.empty());
2963 :
2964 0 : for (const auto & item : constraint_row)
2965 0 : for (unsigned int j=0,
2966 0 : n_col_dofs = cast_int<unsigned int>(col_dofs.size());
2967 0 : j != n_col_dofs; j++)
2968 0 : if (col_dofs[j] == item.first)
2969 0 : matrix(i,j) = -item.second;
2970 : }
2971 : }
2972 : } // end if is constrained...
2973 0 : }
2974 :
2975 :
2976 :
2977 60111280 : void DofMap::constrain_element_vector (DenseVector<Number> & rhs,
2978 : std::vector<dof_id_type> & row_dofs,
2979 : bool) const
2980 : {
2981 4367348 : libmesh_assert_equal_to (rhs.size(), row_dofs.size());
2982 :
2983 : // check for easy return
2984 60111280 : if (this->_dof_constraints.empty())
2985 21979090 : return;
2986 :
2987 : // The constrained RHS is built up as R^T F.
2988 46549654 : DenseMatrix<Number> R;
2989 :
2990 38132190 : this->build_constraint_matrix (R, row_dofs);
2991 :
2992 12730644 : LOG_SCOPE("constrain_elem_vector()", "DofMap");
2993 :
2994 : // It is possible that the vector is not constrained at all.
2995 42377523 : if ((R.m() == rhs.size()) &&
2996 1621220 : (R.n() == row_dofs.size())) // if the RHS is constrained
2997 : {
2998 : // Compute the matrix-vector product
2999 282098 : DenseVector<Number> old_rhs(rhs);
3000 1621220 : R.vector_mult_transpose(rhs, old_rhs);
3001 :
3002 141049 : libmesh_assert_equal_to (row_dofs.size(), rhs.size());
3003 :
3004 22202451 : for (unsigned int i=0,
3005 281894 : n_row_dofs = cast_int<unsigned int>(row_dofs.size());
3006 22484345 : i != n_row_dofs; i++)
3007 22686588 : if (this->is_constrained_dof(row_dofs[i]))
3008 : {
3009 : // If the DOF is constrained
3010 687331 : libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
3011 :
3012 7321145 : rhs(i) = 0;
3013 : }
3014 : } // end if the RHS is constrained.
3015 29714726 : }
3016 :
3017 :
3018 :
3019 20369 : void DofMap::constrain_element_dyad_matrix (DenseVector<Number> & v,
3020 : DenseVector<Number> & w,
3021 : std::vector<dof_id_type> & row_dofs,
3022 : bool) const
3023 : {
3024 1795 : libmesh_assert_equal_to (v.size(), row_dofs.size());
3025 1795 : libmesh_assert_equal_to (w.size(), row_dofs.size());
3026 :
3027 : // check for easy return
3028 20369 : if (this->_dof_constraints.empty())
3029 0 : return;
3030 :
3031 : // The constrained RHS is built up as R^T F.
3032 23959 : DenseMatrix<Number> R;
3033 :
3034 20369 : this->build_constraint_matrix (R, row_dofs);
3035 :
3036 5385 : LOG_SCOPE("cnstrn_elem_dyad_mat()", "DofMap");
3037 :
3038 : // It is possible that the vector is not constrained at all.
3039 22959 : if ((R.m() == v.size()) &&
3040 8489 : (R.n() == row_dofs.size())) // if the RHS is constrained
3041 : {
3042 : // Compute the matrix-vector products
3043 1590 : DenseVector<Number> old_v(v);
3044 1590 : DenseVector<Number> old_w(w);
3045 :
3046 : // compute matrix/vector product
3047 8489 : R.vector_mult_transpose(v, old_v);
3048 8489 : R.vector_mult_transpose(w, old_w);
3049 :
3050 795 : libmesh_assert_equal_to (row_dofs.size(), v.size());
3051 795 : libmesh_assert_equal_to (row_dofs.size(), w.size());
3052 :
3053 : /* Constrain only v, not w. */
3054 :
3055 50829 : for (unsigned int i=0,
3056 1590 : n_row_dofs = cast_int<unsigned int>(row_dofs.size());
3057 52419 : i != n_row_dofs; i++)
3058 48026 : if (this->is_constrained_dof(row_dofs[i]))
3059 : {
3060 : // If the DOF is constrained
3061 916 : libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
3062 :
3063 9974 : v(i) = 0;
3064 : }
3065 : } // end if the RHS is constrained.
3066 16779 : }
3067 :
3068 :
3069 :
3070 2379157 : void DofMap::constrain_nothing (std::vector<dof_id_type> & dofs) const
3071 : {
3072 : // check for easy return
3073 2379157 : if (this->_dof_constraints.empty())
3074 199777 : return;
3075 :
3076 : // All the work is done by \p build_constraint_matrix. We just need
3077 : // a dummy matrix.
3078 2593540 : DenseMatrix<Number> R;
3079 2179380 : this->build_constraint_matrix (R, dofs);
3080 1765220 : }
3081 :
3082 :
3083 :
3084 890246 : void DofMap::enforce_constraints_exactly (const System & system,
3085 : NumericVector<Number> * v,
3086 : bool homogeneous) const
3087 : {
3088 23126 : parallel_object_only();
3089 :
3090 890246 : if (!this->n_constrained_dofs())
3091 312861 : return;
3092 :
3093 42828 : LOG_SCOPE("enforce_constraints_exactly()","DofMap");
3094 :
3095 577385 : if (!v)
3096 5806 : v = system.solution.get();
3097 :
3098 577385 : if (!v->closed())
3099 0 : v->close();
3100 :
3101 14412 : NumericVector<Number> * v_local = nullptr; // will be initialized below
3102 14412 : NumericVector<Number> * v_global = nullptr; // will be initialized below
3103 563381 : std::unique_ptr<NumericVector<Number>> v_built;
3104 577385 : if (v->type() == SERIAL)
3105 : {
3106 2506 : v_built = NumericVector<Number>::build(this->comm());
3107 1253 : v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
3108 1253 : v_built->close();
3109 :
3110 1798756 : for (dof_id_type i=v_built->first_local_index();
3111 1798756 : i<v_built->last_local_index(); i++)
3112 1797503 : v_built->set(i, (*v)(i));
3113 1253 : v_built->close();
3114 0 : v_global = v_built.get();
3115 :
3116 0 : v_local = v;
3117 0 : libmesh_assert (v_local->closed());
3118 : }
3119 576132 : else if (v->type() == PARALLEL)
3120 : {
3121 507778 : v_built = NumericVector<Number>::build(this->comm());
3122 269113 : v_built->init (v->size(), v->local_size(),
3123 : this->get_send_list(), true,
3124 15224 : GHOSTED);
3125 261501 : v->localize(*v_built, this->get_send_list());
3126 261501 : v_built->close();
3127 7612 : v_local = v_built.get();
3128 :
3129 7612 : v_global = v;
3130 : }
3131 314631 : else if (v->type() == GHOSTED)
3132 : {
3133 6800 : v_local = v;
3134 6800 : v_global = v;
3135 : }
3136 : else // unknown v->type()
3137 0 : libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(v->type()));
3138 :
3139 : // We should never hit these asserts because we should error-out in
3140 : // else clause above. Just to be sure we don't try to use v_local
3141 : // and v_global uninitialized...
3142 14412 : libmesh_assert(v_local);
3143 14412 : libmesh_assert(v_global);
3144 14412 : libmesh_assert_equal_to (this, &(system.get_dof_map()));
3145 :
3146 10689956 : for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3147 : {
3148 10112571 : if (!this->local_index(constrained_dof))
3149 1460088 : continue;
3150 :
3151 760468 : Number exact_value = 0;
3152 8450775 : if (!homogeneous)
3153 : {
3154 6624651 : if (auto rhsit = _primal_constraint_values.find(constrained_dof);
3155 590274 : rhsit != _primal_constraint_values.end())
3156 200612 : exact_value = rhsit->second;
3157 : }
3158 16464805 : for (const auto & [dof, val] : constraint_row)
3159 8014030 : exact_value += val * (*v_local)(dof);
3160 :
3161 8450775 : v_global->set(constrained_dof, exact_value);
3162 : }
3163 :
3164 : // If the old vector was serial, we probably need to send our values
3165 : // to other processors
3166 577385 : if (v->type() == SERIAL)
3167 : {
3168 : #ifndef NDEBUG
3169 0 : v_global->close();
3170 : #endif
3171 1253 : v_global->localize (*v);
3172 : }
3173 577385 : v->close();
3174 548969 : }
3175 :
3176 281159 : void DofMap::enforce_constraints_on_residual (const NonlinearImplicitSystem & system,
3177 : NumericVector<Number> * rhs,
3178 : NumericVector<Number> const * solution,
3179 : bool homogeneous) const
3180 : {
3181 5660 : parallel_object_only();
3182 :
3183 281159 : if (!this->n_constrained_dofs())
3184 3265 : return;
3185 :
3186 277798 : if (!rhs)
3187 0 : rhs = system.rhs;
3188 277798 : if (!solution)
3189 0 : solution = system.solution.get();
3190 :
3191 5564 : NumericVector<Number> const * solution_local = nullptr; // will be initialized below
3192 272642 : std::unique_ptr<NumericVector<Number>> solution_built;
3193 277798 : if (solution->type() == SERIAL || solution->type() == GHOSTED)
3194 5564 : solution_local = solution;
3195 0 : else if (solution->type() == PARALLEL)
3196 : {
3197 0 : solution_built = NumericVector<Number>::build(this->comm());
3198 0 : solution_built->init (solution->size(), solution->local_size(),
3199 0 : this->get_send_list(), true, GHOSTED);
3200 0 : solution->localize(*solution_built, this->get_send_list());
3201 0 : solution_built->close();
3202 0 : solution_local = solution_built.get();
3203 : }
3204 : else // unknown solution->type()
3205 0 : libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(solution->type()));
3206 :
3207 : // We should never hit these asserts because we should error-out in
3208 : // else clause above. Just to be sure we don't try to use solution_local
3209 5564 : libmesh_assert(solution_local);
3210 5564 : libmesh_assert_equal_to (this, &(system.get_dof_map()));
3211 :
3212 474070 : for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3213 : {
3214 196272 : if (!this->local_index(constrained_dof))
3215 72529 : continue;
3216 :
3217 8833 : Number exact_value = 0;
3218 201557 : for (const auto & [dof, val] : constraint_row)
3219 82936 : exact_value -= val * (*solution_local)(dof);
3220 118621 : exact_value += (*solution_local)(constrained_dof);
3221 118621 : if (!homogeneous)
3222 : {
3223 0 : if (auto rhsit = _primal_constraint_values.find(constrained_dof);
3224 0 : rhsit != _primal_constraint_values.end())
3225 0 : exact_value += rhsit->second;
3226 : }
3227 :
3228 118621 : rhs->set(constrained_dof, exact_value);
3229 : }
3230 267078 : }
3231 :
3232 158472 : void DofMap::enforce_constraints_on_jacobian (const NonlinearImplicitSystem & system,
3233 : SparseMatrix<Number> * jac) const
3234 : {
3235 3240 : parallel_object_only();
3236 :
3237 158472 : if (!this->n_constrained_dofs())
3238 76 : return;
3239 :
3240 155811 : if (!jac)
3241 0 : jac = system.matrix;
3242 :
3243 3164 : libmesh_assert_equal_to (this, &(system.get_dof_map()));
3244 :
3245 284535 : for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3246 : {
3247 128724 : if (!this->local_index(constrained_dof))
3248 42357 : continue;
3249 :
3250 130893 : for (const auto & j : constraint_row)
3251 47658 : jac->set(constrained_dof, j.first, -j.second);
3252 83235 : jac->set(constrained_dof, constrained_dof, 1);
3253 : }
3254 : }
3255 :
3256 :
3257 60624 : void DofMap::enforce_adjoint_constraints_exactly (NumericVector<Number> & v,
3258 : unsigned int q) const
3259 : {
3260 1896 : parallel_object_only();
3261 :
3262 60624 : if (!this->n_constrained_dofs())
3263 2111 : return;
3264 :
3265 5454 : LOG_SCOPE("enforce_adjoint_constraints_exactly()", "DofMap");
3266 :
3267 1818 : NumericVector<Number> * v_local = nullptr; // will be initialized below
3268 1818 : NumericVector<Number> * v_global = nullptr; // will be initialized below
3269 56695 : std::unique_ptr<NumericVector<Number>> v_built;
3270 58513 : if (v.type() == SERIAL)
3271 : {
3272 198 : v_built = NumericVector<Number>::build(this->comm());
3273 99 : v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
3274 99 : v_built->close();
3275 :
3276 87318 : for (dof_id_type i=v_built->first_local_index();
3277 87318 : i<v_built->last_local_index(); i++)
3278 87219 : v_built->set(i, v(i));
3279 99 : v_built->close();
3280 0 : v_global = v_built.get();
3281 :
3282 0 : v_local = &v;
3283 0 : libmesh_assert (v_local->closed());
3284 : }
3285 58414 : else if (v.type() == PARALLEL)
3286 : {
3287 10920 : v_built = NumericVector<Number>::build(this->comm());
3288 6088 : v_built->init (v.size(), v.local_size(),
3289 628 : this->get_send_list(), true, GHOSTED);
3290 5774 : v.localize(*v_built, this->get_send_list());
3291 5774 : v_built->close();
3292 314 : v_local = v_built.get();
3293 :
3294 314 : v_global = &v;
3295 : }
3296 52640 : else if (v.type() == GHOSTED)
3297 : {
3298 1504 : v_local = &v;
3299 1504 : v_global = &v;
3300 : }
3301 : else // unknown v.type()
3302 0 : libmesh_error_msg("ERROR: Unknown v.type() == " << v.type());
3303 :
3304 : // We should never hit these asserts because we should error-out in
3305 : // else clause above. Just to be sure we don't try to use v_local
3306 : // and v_global uninitialized...
3307 1818 : libmesh_assert(v_local);
3308 1818 : libmesh_assert(v_global);
3309 :
3310 : // Do we have any non_homogeneous constraints?
3311 : const AdjointDofConstraintValues::const_iterator
3312 1818 : adjoint_constraint_map_it = _adjoint_constraint_values.find(q);
3313 : const DofConstraintValueMap * constraint_map =
3314 59943 : (adjoint_constraint_map_it == _adjoint_constraint_values.end()) ?
3315 1430 : nullptr : &adjoint_constraint_map_it->second;
3316 :
3317 444326 : for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3318 : {
3319 385813 : if (!this->local_index(constrained_dof))
3320 72566 : continue;
3321 :
3322 28787 : Number exact_value = 0;
3323 309230 : if (constraint_map)
3324 : {
3325 74368 : if (const auto adjoint_constraint_it = constraint_map->find(constrained_dof);
3326 7168 : adjoint_constraint_it != constraint_map->end())
3327 4476 : exact_value = adjoint_constraint_it->second;
3328 : }
3329 :
3330 523838 : for (const auto & j : constraint_row)
3331 214608 : exact_value += j.second * (*v_local)(j.first);
3332 :
3333 309230 : v_global->set(constrained_dof, exact_value);
3334 : }
3335 :
3336 : // If the old vector was serial, we probably need to send our values
3337 : // to other processors
3338 58513 : if (v.type() == SERIAL)
3339 : {
3340 : #ifndef NDEBUG
3341 0 : v_global->close();
3342 : #endif
3343 99 : v_global->localize (v);
3344 : }
3345 58513 : v.close();
3346 54877 : }
3347 :
3348 :
3349 :
3350 : std::pair<Real, Real>
3351 0 : DofMap::max_constraint_error (const System & system,
3352 : NumericVector<Number> * v) const
3353 : {
3354 0 : if (!v)
3355 0 : v = system.solution.get();
3356 0 : NumericVector<Number> & vec = *v;
3357 :
3358 : // We'll assume the vector is closed
3359 0 : libmesh_assert (vec.closed());
3360 :
3361 0 : Real max_absolute_error = 0., max_relative_error = 0.;
3362 :
3363 0 : const MeshBase & mesh = system.get_mesh();
3364 :
3365 0 : libmesh_assert_equal_to (this, &(system.get_dof_map()));
3366 :
3367 : // indices on each element
3368 0 : std::vector<dof_id_type> local_dof_indices;
3369 :
3370 0 : for (const auto & elem : mesh.active_local_element_ptr_range())
3371 : {
3372 0 : this->dof_indices(elem, local_dof_indices);
3373 0 : std::vector<dof_id_type> raw_dof_indices = local_dof_indices;
3374 :
3375 : // Constraint matrix for each element
3376 0 : DenseMatrix<Number> C;
3377 :
3378 0 : this->build_constraint_matrix (C, local_dof_indices);
3379 :
3380 : // Continue if the element is unconstrained
3381 0 : if (!C.m())
3382 0 : continue;
3383 :
3384 0 : libmesh_assert_equal_to (C.m(), raw_dof_indices.size());
3385 0 : libmesh_assert_equal_to (C.n(), local_dof_indices.size());
3386 :
3387 0 : for (auto i : make_range(C.m()))
3388 : {
3389 : // Recalculate any constrained dof owned by this processor
3390 0 : dof_id_type global_dof = raw_dof_indices[i];
3391 0 : if (this->is_constrained_dof(global_dof) &&
3392 0 : global_dof >= vec.first_local_index() &&
3393 0 : global_dof < vec.last_local_index())
3394 : {
3395 : #ifndef NDEBUG
3396 : DofConstraints::const_iterator
3397 0 : pos = _dof_constraints.find(global_dof);
3398 :
3399 0 : libmesh_assert (pos != _dof_constraints.end());
3400 : #endif
3401 :
3402 0 : Number exact_value = 0;
3403 : DofConstraintValueMap::const_iterator rhsit =
3404 0 : _primal_constraint_values.find(global_dof);
3405 0 : if (rhsit != _primal_constraint_values.end())
3406 0 : exact_value = rhsit->second;
3407 :
3408 0 : for (auto j : make_range(C.n()))
3409 : {
3410 0 : if (local_dof_indices[j] != global_dof)
3411 0 : exact_value += C(i,j) *
3412 0 : vec(local_dof_indices[j]);
3413 : }
3414 :
3415 0 : max_absolute_error = std::max(max_absolute_error,
3416 0 : std::abs(vec(global_dof) - exact_value));
3417 0 : max_relative_error = std::max(max_relative_error,
3418 0 : std::abs(vec(global_dof) - exact_value)
3419 0 : / std::abs(exact_value));
3420 : }
3421 : }
3422 0 : }
3423 :
3424 0 : return std::pair<Real, Real>(max_absolute_error, max_relative_error);
3425 : }
3426 :
3427 :
3428 :
3429 65774545 : void DofMap::build_constraint_matrix (DenseMatrix<Number> & C,
3430 : std::vector<dof_id_type> & elem_dofs,
3431 : const bool called_recursively) const
3432 : {
3433 15088600 : LOG_SCOPE_IF("build_constraint_matrix()", "DofMap", !called_recursively);
3434 :
3435 : // Create a set containing the DOFs we already depend on
3436 : typedef std::set<dof_id_type> RCSet;
3437 7351614 : RCSet dof_set;
3438 :
3439 7351614 : bool we_have_constraints = false;
3440 :
3441 : // Next insert any other dofs the current dofs might be constrained
3442 : // in terms of. Note that in this case we may not be done: Those
3443 : // may in turn depend on others. So, we need to repeat this process
3444 : // in that case until the system depends only on unconstrained
3445 : // degrees of freedom.
3446 436054483 : for (const auto & dof : elem_dofs)
3447 370279938 : if (this->is_constrained_dof(dof))
3448 : {
3449 2656360 : we_have_constraints = true;
3450 :
3451 : // If the DOF is constrained
3452 : DofConstraints::const_iterator
3453 2656360 : pos = _dof_constraints.find(dof);
3454 :
3455 2656360 : libmesh_assert (pos != _dof_constraints.end());
3456 :
3457 2656360 : const DofConstraintRow & constraint_row = pos->second;
3458 :
3459 : // Constraint rows in p refinement may be empty
3460 : //libmesh_assert (!constraint_row.empty());
3461 :
3462 64201690 : for (const auto & item : constraint_row)
3463 35054530 : dof_set.insert (item.first);
3464 : }
3465 :
3466 : // May be safe to return at this point
3467 : // (but remember to stop the perflog)
3468 65774545 : if (!we_have_constraints)
3469 6756734 : return;
3470 :
3471 80552325 : for (const auto & dof : elem_dofs)
3472 6635662 : dof_set.erase (dof);
3473 :
3474 : // If we added any DOFS then we need to do this recursively.
3475 : // It is possible that we just added a DOF that is also
3476 : // constrained!
3477 : //
3478 : // Also, we need to handle the special case of an element having DOFs
3479 : // constrained in terms of other, local DOFs
3480 6999341 : if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
3481 499739 : !called_recursively) // case 2: constrained in terms of our own DOFs
3482 : {
3483 : const unsigned int old_size =
3484 594676 : cast_int<unsigned int>(elem_dofs.size());
3485 :
3486 : // Add new dependency dofs to the end of the current dof set
3487 2952565 : elem_dofs.insert(elem_dofs.end(),
3488 892116 : dof_set.begin(), dof_set.end());
3489 :
3490 : // Now we can build the constraint matrix.
3491 : // Note that resize also zeros for a DenseMatrix<Number>.
3492 3249801 : C.resize (old_size,
3493 : cast_int<unsigned int>(elem_dofs.size()));
3494 :
3495 : // Create the C constraint matrix.
3496 37103010 : for (unsigned int i=0; i != old_size; i++)
3497 36861047 : if (this->is_constrained_dof(elem_dofs[i]))
3498 : {
3499 : // If the DOF is constrained
3500 : DofConstraints::const_iterator
3501 1328180 : pos = _dof_constraints.find(elem_dofs[i]);
3502 :
3503 1328180 : libmesh_assert (pos != _dof_constraints.end());
3504 :
3505 1328180 : const DofConstraintRow & constraint_row = pos->second;
3506 :
3507 : // p refinement creates empty constraint rows
3508 : // libmesh_assert (!constraint_row.empty());
3509 :
3510 32100845 : for (const auto & item : constraint_row)
3511 888952910 : for (unsigned int j=0,
3512 3379092 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
3513 892332002 : j != n_elem_dofs; j++)
3514 960941471 : if (elem_dofs[j] == item.first)
3515 19216607 : C(i,j) = item.second;
3516 : }
3517 : else
3518 : {
3519 18185434 : C(i,i) = 1.;
3520 : }
3521 :
3522 : // May need to do this recursively. It is possible
3523 : // that we just replaced a constrained DOF with another
3524 : // constrained DOF.
3525 3844477 : DenseMatrix<Number> Cnew;
3526 :
3527 3249801 : this->build_constraint_matrix (Cnew, elem_dofs, true);
3528 :
3529 3249801 : if ((C.n() == Cnew.m()) &&
3530 0 : (Cnew.n() == elem_dofs.size())) // If the constraint matrix
3531 0 : C.right_multiply(Cnew); // is constrained...
3532 :
3533 297440 : libmesh_assert_equal_to (C.n(), elem_dofs.size());
3534 2655125 : }
3535 : }
3536 :
3537 :
3538 :
3539 837096 : void DofMap::build_constraint_matrix_and_vector (DenseMatrix<Number> & C,
3540 : DenseVector<Number> & H,
3541 : std::vector<dof_id_type> & elem_dofs,
3542 : int qoi_index,
3543 : const bool called_recursively) const
3544 : {
3545 201168 : LOG_SCOPE_IF("build_constraint_matrix_and_vector()", "DofMap", !called_recursively);
3546 :
3547 : // Create a set containing the DOFs we already depend on
3548 : typedef std::set<dof_id_type> RCSet;
3549 83834 : RCSet dof_set;
3550 :
3551 83834 : bool we_have_constraints = false;
3552 :
3553 : // Next insert any other dofs the current dofs might be constrained
3554 : // in terms of. Note that in this case we may not be done: Those
3555 : // may in turn depend on others. So, we need to repeat this process
3556 : // in that case until the system depends only on unconstrained
3557 : // degrees of freedom.
3558 17099880 : for (const auto & dof : elem_dofs)
3559 16262784 : if (this->is_constrained_dof(dof))
3560 : {
3561 177422 : we_have_constraints = true;
3562 :
3563 : // If the DOF is constrained
3564 : DofConstraints::const_iterator
3565 177422 : pos = _dof_constraints.find(dof);
3566 :
3567 177422 : libmesh_assert (pos != _dof_constraints.end());
3568 :
3569 177422 : const DofConstraintRow & constraint_row = pos->second;
3570 :
3571 : // Constraint rows in p refinement may be empty
3572 : //libmesh_assert (!constraint_row.empty());
3573 :
3574 2095906 : for (const auto & item : constraint_row)
3575 258382 : dof_set.insert (item.first);
3576 : }
3577 :
3578 : // May be safe to return at this point
3579 : // (but remember to stop the perflog)
3580 837096 : if (!we_have_constraints)
3581 50334 : return;
3582 :
3583 6371490 : for (const auto & dof : elem_dofs)
3584 588695 : dof_set.erase (dof);
3585 :
3586 : // If we added any DOFS then we need to do this recursively.
3587 : // It is possible that we just added a DOF that is also
3588 : // constrained!
3589 : //
3590 : // Also, we need to handle the special case of an element having DOFs
3591 : // constrained in terms of other, local DOFs
3592 386305 : if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
3593 30719 : !called_recursively) // case 2: constrained in terms of our own DOFs
3594 : {
3595 16750 : const DofConstraintValueMap * rhs_values = nullptr;
3596 177793 : if (qoi_index < 0)
3597 176113 : rhs_values = &_primal_constraint_values;
3598 1820 : else if (auto it = _adjoint_constraint_values.find(qoi_index);
3599 140 : it != _adjoint_constraint_values.end())
3600 1680 : rhs_values = &it->second;
3601 :
3602 : const unsigned int old_size =
3603 33500 : cast_int<unsigned int>(elem_dofs.size());
3604 :
3605 : // Add new dependency dofs to the end of the current dof set
3606 161043 : elem_dofs.insert(elem_dofs.end(),
3607 50250 : dof_set.begin(), dof_set.end());
3608 :
3609 : // Now we can build the constraint matrix and vector.
3610 : // Note that resize also zeros for a DenseMatrix and DenseVector
3611 177793 : C.resize (old_size,
3612 : cast_int<unsigned int>(elem_dofs.size()));
3613 161043 : H.resize (old_size);
3614 :
3615 : // Create the C constraint matrix.
3616 3150234 : for (unsigned int i=0; i != old_size; i++)
3617 3263886 : if (this->is_constrained_dof(elem_dofs[i]))
3618 : {
3619 : // If the DOF is constrained
3620 : DofConstraints::const_iterator
3621 88711 : pos = _dof_constraints.find(elem_dofs[i]);
3622 :
3623 88711 : libmesh_assert (pos != _dof_constraints.end());
3624 :
3625 88711 : const DofConstraintRow & constraint_row = pos->second;
3626 :
3627 : // p refinement creates empty constraint rows
3628 : // libmesh_assert (!constraint_row.empty());
3629 :
3630 1047953 : for (const auto & item : constraint_row)
3631 1837682 : for (unsigned int j=0,
3632 21240 : n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
3633 1858922 : j != n_elem_dofs; j++)
3634 1872664 : if (elem_dofs[j] == item.first)
3635 139811 : C(i,j) = item.second;
3636 :
3637 918762 : if (rhs_values)
3638 : {
3639 918762 : if (const auto rhsit = rhs_values->find(elem_dofs[i]);
3640 88711 : rhsit != rhs_values->end())
3641 297079 : H(i) = rhsit->second;
3642 : }
3643 : }
3644 : else
3645 : {
3646 1910916 : C(i,i) = 1.;
3647 : }
3648 :
3649 : // May need to do this recursively. It is possible
3650 : // that we just replaced a constrained DOF with another
3651 : // constrained DOF.
3652 211293 : DenseMatrix<Number> Cnew;
3653 177793 : DenseVector<Number> Hnew;
3654 :
3655 177793 : this->build_constraint_matrix_and_vector (Cnew, Hnew, elem_dofs,
3656 : qoi_index, true);
3657 :
3658 177793 : if ((C.n() == Cnew.m()) && // If the constraint matrix
3659 0 : (Cnew.n() == elem_dofs.size())) // is constrained...
3660 : {
3661 : // If x = Cy + h and y = Dz + g
3662 : // Then x = (CD)z + (Cg + h)
3663 0 : C.vector_mult_add(H, 1, Hnew);
3664 :
3665 0 : C.right_multiply(Cnew);
3666 : }
3667 :
3668 16750 : libmesh_assert_equal_to (C.n(), elem_dofs.size());
3669 144293 : }
3670 : }
3671 :
3672 :
3673 335611 : void DofMap::allgather_recursive_constraints(MeshBase & mesh)
3674 : {
3675 : // This function must be run on all processors at once
3676 9620 : parallel_object_only();
3677 :
3678 : // Return immediately if there's nothing to gather
3679 345231 : if (this->n_processors() == 1)
3680 303036 : return;
3681 :
3682 : // We might get to return immediately if none of the processors
3683 : // found any constraints
3684 316740 : unsigned int has_constraints = !_dof_constraints.empty()
3685 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3686 19523 : || !_node_constraints.empty()
3687 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3688 : ;
3689 316740 : this->comm().max(has_constraints);
3690 326360 : if (!has_constraints)
3691 8266 : return;
3692 :
3693 : // If we have heterogeneous adjoint constraints we need to
3694 : // communicate those too.
3695 : const unsigned int max_qoi_num =
3696 32575 : _adjoint_constraint_values.empty() ?
3697 674 : 0 : _adjoint_constraint_values.rbegin()->first+1;
3698 :
3699 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3700 : // We may need to send nodes ahead of data about them
3701 4062 : std::vector<Parallel::Request> packed_range_sends;
3702 :
3703 : // We may be receiving packed_range sends out of order with
3704 : // parallel_sync tags, so make sure they're received correctly.
3705 5416 : Parallel::MessageTag range_tag = this->comm().get_unique_tag();
3706 :
3707 : // We only need to do these sends on a distributed mesh
3708 2708 : const bool dist_mesh = !mesh.is_serial();
3709 : #endif
3710 :
3711 : // We might have calculated constraints for constrained dofs
3712 : // which have support on other processors.
3713 : // Push these out first.
3714 : {
3715 2708 : std::map<processor_id_type, std::set<dof_id_type>> pushed_ids;
3716 :
3717 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3718 2708 : std::map<processor_id_type, std::set<dof_id_type>> pushed_node_ids;
3719 : #endif
3720 :
3721 2708 : const unsigned int sys_num = this->sys_number();
3722 :
3723 : // Collect the constraints to push to each processor
3724 322642 : for (auto & elem : as_range(mesh.active_not_local_elements_begin(),
3725 12538700 : mesh.active_not_local_elements_end()))
3726 : {
3727 6334977 : const unsigned short n_nodes = elem->n_nodes();
3728 :
3729 : // Just checking dof_indices on the foreign element isn't
3730 : // enough. Consider a central hanging node between a coarse
3731 : // Q2/Q1 element and its finer neighbors on a higher-ranked
3732 : // processor. The coarse element's processor will own the node,
3733 : // and will thereby own the pressure dof on that node, despite
3734 : // the fact that that pressure dof doesn't directly exist on the
3735 : // coarse element!
3736 : //
3737 : // So, we loop through dofs manually.
3738 :
3739 : {
3740 6334977 : const unsigned int n_vars = elem->n_vars(sys_num);
3741 14436327 : for (unsigned int v=0; v != n_vars; ++v)
3742 : {
3743 8101350 : const unsigned int n_comp = elem->n_comp(sys_num,v);
3744 10894759 : for (unsigned int c=0; c != n_comp; ++c)
3745 : {
3746 : const unsigned int id =
3747 2793409 : elem->dof_number(sys_num,v,c);
3748 2728684 : if (this->is_constrained_dof(id))
3749 0 : pushed_ids[elem->processor_id()].insert(id);
3750 : }
3751 : }
3752 : }
3753 :
3754 36319528 : for (unsigned short n = 0; n != n_nodes; ++n)
3755 : {
3756 29984551 : const Node & node = elem->node_ref(n);
3757 29984551 : const unsigned int n_vars = node.n_vars(sys_num);
3758 70703611 : for (unsigned int v=0; v != n_vars; ++v)
3759 : {
3760 40719060 : const unsigned int n_comp = node.n_comp(sys_num,v);
3761 84759053 : for (unsigned int c=0; c != n_comp; ++c)
3762 : {
3763 : const unsigned int id =
3764 44039993 : node.dof_number(sys_num,v,c);
3765 42086896 : if (this->is_constrained_dof(id))
3766 255289 : pushed_ids[elem->processor_id()].insert(id);
3767 : }
3768 : }
3769 : }
3770 :
3771 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3772 3339544 : for (unsigned short n = 0; n != n_nodes; ++n)
3773 4222297 : if (this->is_constrained_node(elem->node_ptr(n)))
3774 31443 : pushed_node_ids[elem->processor_id()].insert(elem->node_id(n));
3775 : #endif
3776 29867 : }
3777 :
3778 : // Rewrite those id sets as vectors for sending and receiving,
3779 : // then find the corresponding data for each id, then push it all.
3780 : std::map<processor_id_type, std::vector<dof_id_type>>
3781 2708 : pushed_id_vecs, received_id_vecs;
3782 72775 : for (auto & p : pushed_ids)
3783 40200 : pushed_id_vecs[p.first].assign(p.second.begin(), p.second.end());
3784 :
3785 : std::map<processor_id_type, std::vector<std::vector<std::pair<dof_id_type,Real>>>>
3786 2708 : pushed_keys_vals, received_keys_vals;
3787 2708 : std::map<processor_id_type, std::vector<std::vector<Number>>> pushed_rhss, received_rhss;
3788 72775 : for (auto & p : pushed_id_vecs)
3789 : {
3790 40200 : auto & keys_vals = pushed_keys_vals[p.first];
3791 40868 : keys_vals.reserve(p.second.size());
3792 :
3793 40200 : auto & rhss = pushed_rhss[p.first];
3794 40868 : rhss.reserve(p.second.size());
3795 199031 : for (auto & pushed_id : p.second)
3796 : {
3797 158831 : const DofConstraintRow & row = _dof_constraints[pushed_id];
3798 162575 : keys_vals.emplace_back(row.begin(), row.end());
3799 :
3800 313918 : rhss.push_back(std::vector<Number>(max_qoi_num+1));
3801 3744 : std::vector<Number> & rhs = rhss.back();
3802 : DofConstraintValueMap::const_iterator rhsit =
3803 3744 : _primal_constraint_values.find(pushed_id);
3804 162575 : rhs[max_qoi_num] =
3805 161409 : (rhsit == _primal_constraint_values.end()) ?
3806 2578 : 0 : rhsit->second;
3807 160737 : for (unsigned int q = 0; q != max_qoi_num; ++q)
3808 : {
3809 : AdjointDofConstraintValues::const_iterator adjoint_map_it =
3810 124 : _adjoint_constraint_values.find(q);
3811 :
3812 1906 : if (adjoint_map_it == _adjoint_constraint_values.end())
3813 0 : continue;
3814 :
3815 : const DofConstraintValueMap & constraint_map =
3816 124 : adjoint_map_it->second;
3817 :
3818 : DofConstraintValueMap::const_iterator adj_rhsit =
3819 124 : constraint_map.find(pushed_id);
3820 :
3821 2030 : rhs[q] =
3822 2124 : (adj_rhsit == constraint_map.end()) ?
3823 94 : 0 : adj_rhsit->second;
3824 : }
3825 : }
3826 : }
3827 :
3828 : auto ids_action_functor =
3829 38864 : [& received_id_vecs]
3830 : (processor_id_type pid,
3831 39532 : const std::vector<dof_id_type> & data)
3832 : {
3833 40200 : received_id_vecs[pid] = data;
3834 33243 : };
3835 :
3836 : Parallel::push_parallel_vector_data
3837 32575 : (this->comm(), pushed_id_vecs, ids_action_functor);
3838 :
3839 : auto keys_vals_action_functor =
3840 38864 : [& received_keys_vals]
3841 : (processor_id_type pid,
3842 39532 : const std::vector<std::vector<std::pair<dof_id_type,Real>>> & data)
3843 : {
3844 40200 : received_keys_vals[pid] = data;
3845 33243 : };
3846 :
3847 : Parallel::push_parallel_vector_data
3848 32575 : (this->comm(), pushed_keys_vals, keys_vals_action_functor);
3849 :
3850 : auto rhss_action_functor =
3851 38864 : [& received_rhss]
3852 : (processor_id_type pid,
3853 39532 : const std::vector<std::vector<Number>> & data)
3854 : {
3855 40200 : received_rhss[pid] = data;
3856 33243 : };
3857 :
3858 : Parallel::push_parallel_vector_data
3859 32575 : (this->comm(), pushed_rhss, rhss_action_functor);
3860 :
3861 : // Now we have all the DofConstraint rows and rhs values received
3862 : // from others, so add the DoF constraints that we've been sent
3863 :
3864 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3865 : std::map<processor_id_type, std::vector<dof_id_type>>
3866 2708 : pushed_node_id_vecs, received_node_id_vecs;
3867 4230 : for (auto & p : pushed_node_ids)
3868 1522 : pushed_node_id_vecs[p.first].assign(p.second.begin(), p.second.end());
3869 :
3870 : std::map<processor_id_type, std::vector<std::vector<std::pair<dof_id_type,Real>>>>
3871 2708 : pushed_node_keys_vals, received_node_keys_vals;
3872 2708 : std::map<processor_id_type, std::vector<Point>> pushed_offsets, received_offsets;
3873 :
3874 4230 : for (auto & p : pushed_node_id_vecs)
3875 : {
3876 1522 : const processor_id_type pid = p.first;
3877 :
3878 : // FIXME - this could be an unordered set, given a
3879 : // hash<pointers> specialization
3880 1522 : std::set<const Node *> nodes_requested;
3881 :
3882 1522 : auto & node_keys_vals = pushed_node_keys_vals[pid];
3883 2283 : node_keys_vals.reserve(p.second.size());
3884 :
3885 1522 : auto & offsets = pushed_offsets[pid];
3886 2283 : offsets.reserve(p.second.size());
3887 :
3888 12272 : for (auto & pushed_node_id : p.second)
3889 : {
3890 10750 : const Node * node = mesh.node_ptr(pushed_node_id);
3891 10750 : NodeConstraintRow & row = _node_constraints[node].first;
3892 5375 : const std::size_t row_size = row.size();
3893 : node_keys_vals.push_back
3894 10750 : (std::vector<std::pair<dof_id_type,Real>>());
3895 : std::vector<std::pair<dof_id_type,Real>> & this_node_kv =
3896 5375 : node_keys_vals.back();
3897 10750 : this_node_kv.reserve(row_size);
3898 46032 : for (const auto & j : row)
3899 : {
3900 35282 : this_node_kv.emplace_back(j.first->id(), j.second);
3901 :
3902 : // If we're not sure whether our send
3903 : // destination already has this node, let's give
3904 : // it a copy.
3905 35282 : if (j.first->processor_id() != pid && dist_mesh)
3906 0 : nodes_requested.insert(j.first);
3907 : }
3908 :
3909 10750 : offsets.push_back(_node_constraints[node].second);
3910 :
3911 : }
3912 :
3913 : // Constraining nodes might not even exist on our
3914 : // correspondant's subset of a distributed mesh, so let's
3915 : // make them exist.
3916 1522 : if (dist_mesh)
3917 : {
3918 0 : packed_range_sends.push_back(Parallel::Request());
3919 0 : this->comm().send_packed_range
3920 0 : (pid, &mesh, nodes_requested.begin(), nodes_requested.end(),
3921 0 : packed_range_sends.back(), range_tag);
3922 : }
3923 : }
3924 :
3925 : auto node_ids_action_functor =
3926 : [& received_node_id_vecs]
3927 : (processor_id_type pid,
3928 761 : const std::vector<dof_id_type> & data)
3929 : {
3930 1522 : received_node_id_vecs[pid] = data;
3931 3469 : };
3932 :
3933 : Parallel::push_parallel_vector_data
3934 2708 : (this->comm(), pushed_node_id_vecs, node_ids_action_functor);
3935 :
3936 : auto node_keys_vals_action_functor =
3937 : [& received_node_keys_vals]
3938 : (processor_id_type pid,
3939 761 : const std::vector<std::vector<std::pair<dof_id_type,Real>>> & data)
3940 : {
3941 1522 : received_node_keys_vals[pid] = data;
3942 3469 : };
3943 :
3944 : Parallel::push_parallel_vector_data
3945 2708 : (this->comm(), pushed_node_keys_vals,
3946 : node_keys_vals_action_functor);
3947 :
3948 : auto node_offsets_action_functor =
3949 : [& received_offsets]
3950 : (processor_id_type pid,
3951 761 : const std::vector<Point> & data)
3952 : {
3953 1522 : received_offsets[pid] = data;
3954 3469 : };
3955 :
3956 : Parallel::push_parallel_vector_data
3957 2708 : (this->comm(), pushed_offsets, node_offsets_action_functor);
3958 :
3959 : #endif
3960 :
3961 : // Add all the dof constraints that I've been sent
3962 72775 : for (auto & [pid, pushed_ids_to_me] : received_id_vecs)
3963 : {
3964 668 : libmesh_assert(received_keys_vals.count(pid));
3965 668 : libmesh_assert(received_rhss.count(pid));
3966 40200 : const auto & pushed_keys_vals_to_me = received_keys_vals.at(pid);
3967 40200 : const auto & pushed_rhss_to_me = received_rhss.at(pid);
3968 :
3969 668 : libmesh_assert_equal_to (pushed_ids_to_me.size(),
3970 : pushed_keys_vals_to_me.size());
3971 668 : libmesh_assert_equal_to (pushed_ids_to_me.size(),
3972 : pushed_rhss_to_me.size());
3973 :
3974 199031 : for (auto i : index_range(pushed_ids_to_me))
3975 : {
3976 162575 : dof_id_type constrained = pushed_ids_to_me[i];
3977 :
3978 : // If we don't already have a constraint for this dof,
3979 : // add the one we were sent
3980 104355 : if (!this->is_constrained_dof(constrained))
3981 : {
3982 52594 : DofConstraintRow & row = _dof_constraints[constrained];
3983 190277 : for (auto & kv : pushed_keys_vals_to_me[i])
3984 : {
3985 1678 : libmesh_assert_less(kv.first, this->n_dofs());
3986 136752 : row[kv.first] = kv.second;
3987 : }
3988 :
3989 53525 : const Number primal_rhs = pushed_rhss_to_me[i][max_qoi_num];
3990 :
3991 931 : if (libmesh_isnan(primal_rhs))
3992 0 : libmesh_assert(pushed_keys_vals_to_me[i].empty());
3993 :
3994 52594 : if (primal_rhs != Number(0))
3995 11164 : _primal_constraint_values[constrained] = primal_rhs;
3996 : else
3997 541 : _primal_constraint_values.erase(constrained);
3998 :
3999 52594 : for (unsigned int q = 0; q != max_qoi_num; ++q)
4000 : {
4001 : AdjointDofConstraintValues::iterator adjoint_map_it =
4002 0 : _adjoint_constraint_values.find(q);
4003 :
4004 0 : const Number adj_rhs = pushed_rhss_to_me[i][q];
4005 :
4006 0 : if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
4007 : adj_rhs == Number(0))
4008 0 : continue;
4009 :
4010 0 : if (adjoint_map_it == _adjoint_constraint_values.end())
4011 0 : adjoint_map_it = _adjoint_constraint_values.emplace
4012 0 : (q, DofConstraintValueMap()).first;
4013 :
4014 : DofConstraintValueMap & constraint_map =
4015 0 : adjoint_map_it->second;
4016 :
4017 0 : if (adj_rhs != Number(0))
4018 0 : constraint_map[constrained] = adj_rhs;
4019 : else
4020 0 : constraint_map.erase(constrained);
4021 : }
4022 : }
4023 : }
4024 : }
4025 :
4026 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4027 : // Add all the node constraints that I've been sent
4028 4230 : for (auto & [pid, pushed_node_ids_to_me] : received_node_id_vecs)
4029 : {
4030 : // Before we act on any new constraint rows, we may need to
4031 : // make sure we have all the nodes involved!
4032 1522 : if (dist_mesh)
4033 0 : this->comm().receive_packed_range
4034 0 : (pid, &mesh, null_output_iterator<Node>(),
4035 : (Node**)nullptr, range_tag);
4036 :
4037 761 : libmesh_assert(received_node_keys_vals.count(pid));
4038 761 : libmesh_assert(received_offsets.count(pid));
4039 1522 : const auto & pushed_node_keys_vals_to_me = received_node_keys_vals.at(pid);
4040 1522 : const auto & pushed_offsets_to_me = received_offsets.at(pid);
4041 :
4042 761 : libmesh_assert_equal_to (pushed_node_ids_to_me.size(),
4043 : pushed_node_keys_vals_to_me.size());
4044 761 : libmesh_assert_equal_to (pushed_node_ids_to_me.size(),
4045 : pushed_offsets_to_me.size());
4046 :
4047 12272 : for (auto i : index_range(pushed_node_ids_to_me))
4048 : {
4049 10750 : dof_id_type constrained_id = pushed_node_ids_to_me[i];
4050 :
4051 : // If we don't already have a constraint for this node,
4052 : // add the one we were sent
4053 10750 : const Node * constrained = mesh.node_ptr(constrained_id);
4054 9025 : if (!this->is_constrained_node(constrained))
4055 : {
4056 3450 : NodeConstraintRow & row = _node_constraints[constrained].first;
4057 18086 : for (auto & kv : pushed_node_keys_vals_to_me[i])
4058 : {
4059 12911 : const Node * key_node = mesh.node_ptr(kv.first);
4060 5601 : libmesh_assert(key_node);
4061 12911 : row[key_node] = kv.second;
4062 : }
4063 5175 : _node_constraints[constrained].second = pushed_offsets_to_me[i];
4064 : }
4065 : }
4066 : }
4067 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4068 : }
4069 :
4070 : // Now start checking for any other constraints we need
4071 : // to know about, requesting them recursively.
4072 :
4073 : // Create sets containing the DOFs and nodes we already depend on
4074 : typedef std::set<dof_id_type> DoF_RCSet;
4075 2708 : DoF_RCSet unexpanded_dofs;
4076 :
4077 1043056 : for (const auto & i : _dof_constraints)
4078 1010481 : unexpanded_dofs.insert(i.first);
4079 :
4080 : // Gather all the dof constraints we need
4081 32575 : this->gather_constraints(mesh, unexpanded_dofs, false);
4082 :
4083 : // Gather all the node constraints we need
4084 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4085 : typedef std::set<const Node *> Node_RCSet;
4086 2708 : Node_RCSet unexpanded_nodes;
4087 :
4088 177026 : for (const auto & i : _node_constraints)
4089 174318 : unexpanded_nodes.insert(i.first);
4090 :
4091 : // We have to keep recursing while the unexpanded set is
4092 : // nonempty on *any* processor
4093 2708 : bool unexpanded_set_nonempty = !unexpanded_nodes.empty();
4094 2708 : this->comm().max(unexpanded_set_nonempty);
4095 :
4096 21092 : while (unexpanded_set_nonempty)
4097 : {
4098 : // Let's make sure we don't lose sync in this loop.
4099 9192 : parallel_object_only();
4100 :
4101 : // Request sets
4102 18384 : Node_RCSet node_request_set;
4103 :
4104 : // Request sets to send to each processor
4105 : std::map<processor_id_type, std::vector<dof_id_type>>
4106 18384 : requested_node_ids;
4107 :
4108 : // And the sizes of each
4109 18384 : std::map<processor_id_type, dof_id_type> node_ids_on_proc;
4110 :
4111 : // Fill (and thereby sort and uniq!) the main request sets
4112 235224 : for (const auto & i : unexpanded_nodes)
4113 : {
4114 216840 : NodeConstraintRow & row = _node_constraints[i].first;
4115 959675 : for (const auto & j : row)
4116 : {
4117 742835 : const Node * const node = j.first;
4118 340586 : libmesh_assert(node);
4119 :
4120 : // If it's non-local and we haven't already got a
4121 : // constraint for it, we might need to ask for one
4122 1223401 : if ((node->processor_id() != this->processor_id()) &&
4123 78317 : !_node_constraints.count(node))
4124 25977 : node_request_set.insert(node);
4125 : }
4126 : }
4127 :
4128 : // Clear the unexpanded constraint sets; we're about to expand
4129 : // them
4130 9192 : unexpanded_nodes.clear();
4131 :
4132 : // Count requests by processor
4133 61692 : for (const auto & node : node_request_set)
4134 : {
4135 21654 : libmesh_assert(node);
4136 21654 : libmesh_assert_less (node->processor_id(), this->n_processors());
4137 43308 : node_ids_on_proc[node->processor_id()]++;
4138 : }
4139 :
4140 31502 : for (auto pair : node_ids_on_proc)
4141 13118 : requested_node_ids[pair.first].reserve(pair.second);
4142 :
4143 : // Prepare each processor's request set
4144 61692 : for (const auto & node : node_request_set)
4145 43308 : requested_node_ids[node->processor_id()].push_back(node->id());
4146 :
4147 : typedef std::vector<std::pair<dof_id_type, Real>> row_datum;
4148 :
4149 : auto node_row_gather_functor =
4150 : [this,
4151 : & mesh,
4152 : dist_mesh,
4153 : & packed_range_sends,
4154 : & range_tag]
4155 : (processor_id_type pid,
4156 : const std::vector<dof_id_type> & ids,
4157 163286 : std::vector<row_datum> & data)
4158 : {
4159 : // FIXME - this could be an unordered set, given a
4160 : // hash<pointers> specialization
4161 13118 : std::set<const Node *> nodes_requested;
4162 :
4163 : // Fill those requests
4164 13118 : const std::size_t query_size = ids.size();
4165 :
4166 13118 : data.resize(query_size);
4167 56426 : for (std::size_t i=0; i != query_size; ++i)
4168 : {
4169 43308 : dof_id_type constrained_id = ids[i];
4170 43308 : const Node * constrained_node = mesh.node_ptr(constrained_id);
4171 21654 : if (_node_constraints.count(constrained_node))
4172 : {
4173 42522 : const NodeConstraintRow & row = _node_constraints[constrained_node].first;
4174 21261 : std::size_t row_size = row.size();
4175 63783 : data[i].reserve(row_size);
4176 185359 : for (const auto & j : row)
4177 : {
4178 142837 : const Node * node = j.first;
4179 224558 : data[i].emplace_back(node->id(), j.second);
4180 :
4181 : // If we're not sure whether our send
4182 : // destination already has this node, let's give
4183 : // it a copy.
4184 142837 : if (node->processor_id() != pid && dist_mesh)
4185 0 : nodes_requested.insert(node);
4186 :
4187 : // We can have 0 nodal constraint
4188 : // coefficients, where no Lagrange constraint
4189 : // exists but non-Lagrange basis constraints
4190 : // might.
4191 : // libmesh_assert(j.second);
4192 : }
4193 : }
4194 : else
4195 : {
4196 : // We have to distinguish "constraint with no
4197 : // constraining nodes" (e.g. due to user node
4198 : // constraint equations) from "no constraint".
4199 : // We'll use invalid_id for the latter.
4200 1179 : data[i].emplace_back(DofObject::invalid_id, Real(0));
4201 : }
4202 : }
4203 :
4204 : // Constraining nodes might not even exist on our
4205 : // correspondant's subset of a distributed mesh, so let's
4206 : // make them exist.
4207 13118 : if (dist_mesh)
4208 : {
4209 0 : packed_range_sends.push_back(Parallel::Request());
4210 0 : this->comm().send_packed_range
4211 0 : (pid, &mesh, nodes_requested.begin(), nodes_requested.end(),
4212 0 : packed_range_sends.back(), range_tag);
4213 : }
4214 13118 : };
4215 :
4216 : typedef Point node_rhs_datum;
4217 :
4218 : auto node_rhs_gather_functor =
4219 : [this,
4220 : & mesh]
4221 : (processor_id_type,
4222 : const std::vector<dof_id_type> & ids,
4223 163910 : std::vector<node_rhs_datum> & data)
4224 : {
4225 : // Fill those requests
4226 13118 : const std::size_t query_size = ids.size();
4227 :
4228 13118 : data.resize(query_size);
4229 56426 : for (std::size_t i=0; i != query_size; ++i)
4230 : {
4231 43308 : dof_id_type constrained_id = ids[i];
4232 43308 : const Node * constrained_node = mesh.node_ptr(constrained_id);
4233 21654 : if (_node_constraints.count(constrained_node))
4234 42522 : data[i] = _node_constraints[constrained_node].second;
4235 : else
4236 1179 : data[i](0) = std::numeric_limits<Real>::quiet_NaN();
4237 : }
4238 22310 : };
4239 :
4240 : auto node_row_action_functor =
4241 : [this,
4242 : & mesh,
4243 : dist_mesh,
4244 : & range_tag,
4245 : & unexpanded_nodes]
4246 : (processor_id_type pid,
4247 : const std::vector<dof_id_type> & ids,
4248 235096 : const std::vector<row_datum> & data)
4249 : {
4250 : // Before we act on any new constraint rows, we may need to
4251 : // make sure we have all the nodes involved!
4252 13118 : if (dist_mesh)
4253 0 : this->comm().receive_packed_range
4254 0 : (pid, &mesh, null_output_iterator<Node>(),
4255 : (Node**)nullptr, range_tag);
4256 :
4257 : // Add any new constraint rows we've found
4258 13118 : const std::size_t query_size = ids.size();
4259 :
4260 56426 : for (std::size_t i=0; i != query_size; ++i)
4261 : {
4262 43308 : const dof_id_type constrained_id = ids[i];
4263 :
4264 : // An empty row is an constraint with an empty row; for
4265 : // no constraint we use a "no row" placeholder
4266 64962 : if (data[i].empty())
4267 : {
4268 0 : const Node * constrained_node = mesh.node_ptr(constrained_id);
4269 0 : NodeConstraintRow & row = _node_constraints[constrained_node].first;
4270 0 : row.clear();
4271 : }
4272 43308 : else if (data[i][0].first != DofObject::invalid_id)
4273 : {
4274 42522 : const Node * constrained_node = mesh.node_ptr(constrained_id);
4275 42522 : NodeConstraintRow & row = _node_constraints[constrained_node].first;
4276 21261 : row.clear();
4277 206620 : for (auto & pair : data[i])
4278 : {
4279 : const Node * key_node =
4280 142837 : mesh.node_ptr(pair.first);
4281 61116 : libmesh_assert(key_node);
4282 142837 : row[key_node] = pair.second;
4283 : }
4284 :
4285 : // And prepare to check for more recursive constraints
4286 21261 : unexpanded_nodes.insert(constrained_node);
4287 : }
4288 : }
4289 13118 : };
4290 :
4291 : auto node_rhs_action_functor =
4292 : [this,
4293 : & mesh]
4294 : (processor_id_type,
4295 : const std::vector<dof_id_type> & ids,
4296 71914 : const std::vector<node_rhs_datum> & data)
4297 : {
4298 : // Add rhs data for any new node constraint rows we've found
4299 13118 : const std::size_t query_size = ids.size();
4300 :
4301 56426 : for (std::size_t i=0; i != query_size; ++i)
4302 : {
4303 43308 : dof_id_type constrained_id = ids[i];
4304 43308 : const Node * constrained_node = mesh.node_ptr(constrained_id);
4305 :
4306 64962 : if (!libmesh_isnan(data[i](0)))
4307 42522 : _node_constraints[constrained_node].second = data[i];
4308 : else
4309 393 : _node_constraints.erase(constrained_node);
4310 : }
4311 22310 : };
4312 :
4313 : // Now request node constraint rows from other processors
4314 9192 : row_datum * node_row_ex = nullptr;
4315 : Parallel::pull_parallel_vector_data
4316 18384 : (this->comm(), requested_node_ids, node_row_gather_functor,
4317 : node_row_action_functor, node_row_ex);
4318 :
4319 : // And request node constraint right hand sides from other procesors
4320 9192 : node_rhs_datum * node_rhs_ex = nullptr;
4321 : Parallel::pull_parallel_vector_data
4322 18384 : (this->comm(), requested_node_ids, node_rhs_gather_functor,
4323 : node_rhs_action_functor, node_rhs_ex);
4324 :
4325 :
4326 : // We have to keep recursing while the unexpanded set is
4327 : // nonempty on *any* processor
4328 18384 : unexpanded_set_nonempty = !unexpanded_nodes.empty();
4329 18384 : this->comm().max(unexpanded_set_nonempty);
4330 : }
4331 2708 : Parallel::wait(packed_range_sends);
4332 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4333 : }
4334 :
4335 :
4336 :
4337 335611 : void DofMap::process_constraints (MeshBase & mesh)
4338 : {
4339 : // We've computed our local constraints, but they may depend on
4340 : // non-local constraints that we'll need to take into account.
4341 335611 : this->allgather_recursive_constraints(mesh);
4342 :
4343 335611 : if (_error_on_constraint_loop)
4344 : {
4345 : // Optionally check for constraint loops and throw an error
4346 : // if they're detected. We always do this check below in dbg/devel
4347 : // mode but here we optionally do it in opt mode as well.
4348 71 : check_for_constraint_loops();
4349 : }
4350 :
4351 : // Adjoints will be constrained where the primal is
4352 : // Therefore, we will expand the adjoint_constraint_values
4353 : // map whenever the primal_constraint_values map is expanded
4354 :
4355 : // First, figure out the total number of QoIs
4356 : const unsigned int max_qoi_num =
4357 335540 : _adjoint_constraint_values.empty() ?
4358 712 : 0 : _adjoint_constraint_values.rbegin()->first+1;
4359 :
4360 : // Create a set containing the DOFs we already depend on
4361 : typedef std::set<dof_id_type> RCSet;
4362 19236 : RCSet unexpanded_set;
4363 :
4364 1499456 : for (const auto & i : _dof_constraints)
4365 1163916 : unexpanded_set.insert(i.first);
4366 :
4367 362776 : while (!unexpanded_set.empty())
4368 232952 : for (RCSet::iterator i = unexpanded_set.begin();
4369 1248944 : i != unexpanded_set.end(); /* nothing */)
4370 : {
4371 : // If the DOF is constrained
4372 : DofConstraints::iterator
4373 103470 : pos = _dof_constraints.find(*i);
4374 :
4375 103470 : libmesh_assert (pos != _dof_constraints.end());
4376 :
4377 1221708 : DofConstraintRow & constraint_row = pos->second;
4378 :
4379 : DofConstraintValueMap::iterator rhsit =
4380 103470 : _primal_constraint_values.find(*i);
4381 1270490 : Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
4382 89707 : 0 : rhsit->second;
4383 :
4384 : // A vector of DofConstraintValueMaps for each adjoint variable
4385 206940 : std::vector<DofConstraintValueMap::iterator> adjoint_rhs_iterators;
4386 1221708 : adjoint_rhs_iterators.resize(max_qoi_num);
4387 :
4388 : // Another to hold the adjoint constraint rhs
4389 1325178 : std::vector<Number> adjoint_constraint_rhs(max_qoi_num, 0.0);
4390 :
4391 : // Find and gather recursive constraints for each adjoint variable
4392 1242189 : for (auto & adjoint_map : _adjoint_constraint_values)
4393 : {
4394 20481 : const std::size_t q = adjoint_map.first;
4395 20481 : adjoint_rhs_iterators[q] = adjoint_map.second.find(*i);
4396 :
4397 22911 : adjoint_constraint_rhs[q] =
4398 23661 : (adjoint_rhs_iterators[q] == adjoint_map.second.end()) ?
4399 750 : 0 : adjoint_rhs_iterators[q]->second;
4400 : }
4401 :
4402 206940 : std::vector<dof_id_type> constraints_to_expand;
4403 :
4404 3141191 : for (const auto & item : constraint_row)
4405 3452631 : if (item.first != *i && this->is_constrained_dof(item.first))
4406 : {
4407 67909 : unexpanded_set.insert(item.first);
4408 67909 : constraints_to_expand.push_back(item.first);
4409 : }
4410 :
4411 1289617 : for (const auto & expandable : constraints_to_expand)
4412 : {
4413 67909 : const Real this_coef = constraint_row[expandable];
4414 :
4415 : DofConstraints::const_iterator
4416 5529 : subpos = _dof_constraints.find(expandable);
4417 :
4418 5529 : libmesh_assert (subpos != _dof_constraints.end());
4419 :
4420 5529 : const DofConstraintRow & subconstraint_row = subpos->second;
4421 :
4422 72617 : for (const auto & item : subconstraint_row)
4423 : {
4424 : // Assert that the constraint does not form a cycle.
4425 607 : libmesh_assert(item.first != expandable);
4426 4708 : constraint_row[item.first] += item.second * this_coef;
4427 : }
4428 :
4429 67909 : if (auto subrhsit = _primal_constraint_values.find(expandable);
4430 5529 : subrhsit != _primal_constraint_values.end())
4431 19326 : constraint_rhs += subrhsit->second * this_coef;
4432 :
4433 : // Find and gather recursive constraints for each adjoint variable
4434 67909 : for (const auto & adjoint_map : _adjoint_constraint_values)
4435 : {
4436 0 : if (auto adjoint_subrhsit = adjoint_map.second.find(expandable);
4437 0 : adjoint_subrhsit != adjoint_map.second.end())
4438 0 : adjoint_constraint_rhs[adjoint_map.first] += adjoint_subrhsit->second * this_coef;
4439 : }
4440 :
4441 5529 : constraint_row.erase(expandable);
4442 : }
4443 :
4444 1221708 : if (rhsit == _primal_constraint_values.end())
4445 : {
4446 584445 : if (constraint_rhs != Number(0))
4447 11050 : _primal_constraint_values[*i] = constraint_rhs;
4448 : else
4449 53918 : _primal_constraint_values.erase(*i);
4450 : }
4451 : else
4452 : {
4453 568659 : if (constraint_rhs != Number(0))
4454 139496 : rhsit->second = constraint_rhs;
4455 : else
4456 431633 : _primal_constraint_values.erase(rhsit);
4457 : }
4458 :
4459 : // Finally fill in the adjoint constraints for each adjoint variable if possible
4460 1242189 : for (auto & adjoint_map : _adjoint_constraint_values)
4461 : {
4462 20481 : const std::size_t q = adjoint_map.first;
4463 :
4464 22911 : if(adjoint_rhs_iterators[q] == adjoint_map.second.end())
4465 : {
4466 11390 : if (adjoint_constraint_rhs[q] != Number(0))
4467 0 : (adjoint_map.second)[*i] = adjoint_constraint_rhs[q];
4468 : else
4469 1680 : adjoint_map.second.erase(*i);
4470 : }
4471 : else
4472 : {
4473 9281 : if (adjoint_constraint_rhs[q] != Number(0))
4474 6810 : adjoint_rhs_iterators[q]->second = adjoint_constraint_rhs[q];
4475 : else
4476 2106 : adjoint_map.second.erase(adjoint_rhs_iterators[q]);
4477 : }
4478 : }
4479 :
4480 1221708 : if (constraints_to_expand.empty())
4481 1082779 : i = unexpanded_set.erase(i);
4482 : else
4483 3225 : ++i;
4484 : }
4485 :
4486 : // In parallel we can't guarantee that nodes/dofs which constrain
4487 : // others are on processors which are aware of that constraint, yet
4488 : // we need such awareness for sparsity pattern generation. So send
4489 : // other processors any constraints they might need to know about.
4490 335540 : this->scatter_constraints(mesh);
4491 :
4492 : // Now that we have our root constraint dependencies sorted out, add
4493 : // them to the send_list
4494 335540 : this->add_constraints_to_send_list(mesh);
4495 335540 : }
4496 :
4497 :
4498 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
4499 0 : void DofMap::check_for_cyclic_constraints()
4500 : {
4501 : // Eventually make this officially libmesh_deprecated();
4502 0 : check_for_constraint_loops();
4503 0 : }
4504 :
4505 71 : void DofMap::check_for_constraint_loops()
4506 : {
4507 : // Create a set containing the DOFs we already depend on
4508 : typedef std::set<dof_id_type> RCSet;
4509 4 : RCSet unexpanded_set;
4510 :
4511 : // Use dof_constraints_copy in this method so that we don't
4512 : // mess with _dof_constraints.
4513 4 : DofConstraints dof_constraints_copy = _dof_constraints;
4514 :
4515 213 : for (const auto & i : dof_constraints_copy)
4516 142 : unexpanded_set.insert(i.first);
4517 :
4518 71 : while (!unexpanded_set.empty())
4519 73 : for (RCSet::iterator i = unexpanded_set.begin();
4520 142 : i != unexpanded_set.end(); /* nothing */)
4521 : {
4522 : // If the DOF is constrained
4523 : DofConstraints::iterator
4524 4 : pos = dof_constraints_copy.find(*i);
4525 :
4526 4 : libmesh_assert (pos != dof_constraints_copy.end());
4527 :
4528 142 : DofConstraintRow & constraint_row = pos->second;
4529 :
4530 : // Comment out "rhs" parts of this method copied from process_constraints
4531 : // DofConstraintValueMap::iterator rhsit =
4532 : // _primal_constraint_values.find(*i);
4533 : // Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
4534 : // 0 : rhsit->second;
4535 :
4536 8 : std::vector<dof_id_type> constraints_to_expand;
4537 :
4538 284 : for (const auto & item : constraint_row)
4539 142 : if (item.first != *i && this->is_constrained_dof(item.first))
4540 : {
4541 142 : unexpanded_set.insert(item.first);
4542 142 : constraints_to_expand.push_back(item.first);
4543 : }
4544 :
4545 213 : for (const auto & expandable : constraints_to_expand)
4546 : {
4547 142 : const Real this_coef = constraint_row[expandable];
4548 :
4549 : DofConstraints::const_iterator
4550 4 : subpos = dof_constraints_copy.find(expandable);
4551 :
4552 4 : libmesh_assert (subpos != dof_constraints_copy.end());
4553 :
4554 4 : const DofConstraintRow & subconstraint_row = subpos->second;
4555 :
4556 213 : for (const auto & item : subconstraint_row)
4557 : {
4558 213 : libmesh_error_msg_if(item.first == expandable, "Constraint loop detected");
4559 :
4560 71 : constraint_row[item.first] += item.second * this_coef;
4561 : }
4562 :
4563 : // Comment out "rhs" parts of this method copied from process_constraints
4564 : // DofConstraintValueMap::const_iterator subrhsit =
4565 : // _primal_constraint_values.find(expandable);
4566 : // if (subrhsit != _primal_constraint_values.end())
4567 : // constraint_rhs += subrhsit->second * this_coef;
4568 :
4569 2 : constraint_row.erase(expandable);
4570 : }
4571 :
4572 : // Comment out "rhs" parts of this method copied from process_constraints
4573 : // if (rhsit == _primal_constraint_values.end())
4574 : // {
4575 : // if (constraint_rhs != Number(0))
4576 : // _primal_constraint_values[*i] = constraint_rhs;
4577 : // else
4578 : // _primal_constraint_values.erase(*i);
4579 : // }
4580 : // else
4581 : // {
4582 : // if (constraint_rhs != Number(0))
4583 : // rhsit->second = constraint_rhs;
4584 : // else
4585 : // _primal_constraint_values.erase(rhsit);
4586 : // }
4587 :
4588 71 : if (constraints_to_expand.empty())
4589 0 : i = unexpanded_set.erase(i);
4590 : else
4591 2 : ++i;
4592 : }
4593 0 : }
4594 : #else
4595 : void DofMap::check_for_constraint_loops() {}
4596 : void DofMap::check_for_cyclic_constraints()
4597 : {
4598 : // Do nothing
4599 : }
4600 : #endif
4601 :
4602 :
4603 335540 : void DofMap::scatter_constraints(MeshBase & mesh)
4604 : {
4605 : // At this point each processor with a constrained node knows
4606 : // the corresponding constraint row, but we also need each processor
4607 : // with a constrainer node to know the corresponding row(s).
4608 :
4609 : // This function must be run on all processors at once
4610 9618 : parallel_object_only();
4611 :
4612 : // Return immediately if there's nothing to gather
4613 345158 : if (this->n_processors() == 1)
4614 303034 : return;
4615 :
4616 : // We might get to return immediately if none of the processors
4617 : // found any constraints
4618 316673 : unsigned int has_constraints = !_dof_constraints.empty()
4619 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4620 19531 : || !_node_constraints.empty()
4621 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4622 : ;
4623 316673 : this->comm().max(has_constraints);
4624 326291 : if (!has_constraints)
4625 8266 : return;
4626 :
4627 : // We may be receiving packed_range sends out of order with
4628 : // parallel_sync tags, so make sure they're received correctly.
4629 35210 : Parallel::MessageTag range_tag = this->comm().get_unique_tag();
4630 :
4631 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4632 2704 : std::map<processor_id_type, std::set<dof_id_type>> pushed_node_ids;
4633 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4634 :
4635 2704 : std::map<processor_id_type, std::set<dof_id_type>> pushed_ids;
4636 :
4637 : // Collect the dof constraints I need to push to each processor
4638 1352 : dof_id_type constrained_proc_id = 0;
4639 1044701 : for (const auto & [constrained, row] : _dof_constraints)
4640 : {
4641 1201918 : while (constrained >= _end_df[constrained_proc_id])
4642 90472 : constrained_proc_id++;
4643 :
4644 1110937 : if (constrained_proc_id != this->processor_id())
4645 100725 : continue;
4646 :
4647 2293157 : for (auto & j : row)
4648 : {
4649 1384268 : const dof_id_type constraining = j.first;
4650 :
4651 1384268 : processor_id_type constraining_proc_id = 0;
4652 5312279 : while (constraining >= _end_df[constraining_proc_id])
4653 3703428 : constraining_proc_id++;
4654 :
4655 1562741 : if (constraining_proc_id != this->processor_id() &&
4656 27053 : constraining_proc_id != constrained_proc_id)
4657 213325 : pushed_ids[constraining_proc_id].insert(constrained);
4658 : }
4659 : }
4660 :
4661 : // Pack the dof constraint rows and rhs's to push
4662 :
4663 : std::map<processor_id_type,
4664 : std::vector<std::vector<std::pair<dof_id_type, Real>>>>
4665 2704 : pushed_keys_vals, pushed_keys_vals_to_me;
4666 :
4667 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, Number>>>
4668 2704 : pushed_ids_rhss, pushed_ids_rhss_to_me;
4669 :
4670 : auto gather_ids =
4671 59604 : [this,
4672 : & pushed_ids,
4673 : & pushed_keys_vals,
4674 : & pushed_ids_rhss]
4675 298421 : ()
4676 : {
4677 103967 : for (const auto & [pid, pid_ids] : pushed_ids)
4678 : {
4679 684 : const std::size_t ids_size = pid_ids.size();
4680 : std::vector<std::vector<std::pair<dof_id_type, Real>>> &
4681 38955 : keys_vals = pushed_keys_vals[pid];
4682 : std::vector<std::pair<dof_id_type,Number>> &
4683 38955 : ids_rhss = pushed_ids_rhss[pid];
4684 38955 : keys_vals.resize(ids_size);
4685 38955 : ids_rhss.resize(ids_size);
4686 :
4687 : std::size_t push_i;
4688 684 : std::set<dof_id_type>::const_iterator it;
4689 80303 : for (push_i = 0, it = pid_ids.begin();
4690 291272 : it != pid_ids.end(); ++push_i, ++it)
4691 : {
4692 252317 : const dof_id_type constrained = *it;
4693 252317 : DofConstraintRow & row = _dof_constraints[constrained];
4694 252317 : keys_vals[push_i].assign(row.begin(), row.end());
4695 :
4696 : DofConstraintValueMap::const_iterator rhsit =
4697 21016 : _primal_constraint_values.find(constrained);
4698 252317 : ids_rhss[push_i].first = constrained;
4699 252317 : ids_rhss[push_i].second =
4700 273666 : (rhsit == _primal_constraint_values.end()) ?
4701 333 : 0 : rhsit->second;
4702 : }
4703 : }
4704 65012 : };
4705 :
4706 32506 : gather_ids();
4707 :
4708 : auto ids_rhss_action_functor =
4709 37587 : [& pushed_ids_rhss_to_me]
4710 : (processor_id_type pid,
4711 38271 : const std::vector<std::pair<dof_id_type, Number>> & data)
4712 : {
4713 38955 : pushed_ids_rhss_to_me[pid] = data;
4714 70777 : };
4715 :
4716 : auto keys_vals_action_functor =
4717 37587 : [& pushed_keys_vals_to_me]
4718 : (processor_id_type pid,
4719 38271 : const std::vector<std::vector<std::pair<dof_id_type, Real>>> & data)
4720 : {
4721 38955 : pushed_keys_vals_to_me[pid] = data;
4722 33190 : };
4723 :
4724 : Parallel::push_parallel_vector_data
4725 32506 : (this->comm(), pushed_ids_rhss, ids_rhss_action_functor);
4726 : Parallel::push_parallel_vector_data
4727 32506 : (this->comm(), pushed_keys_vals, keys_vals_action_functor);
4728 :
4729 : // Now work on traded dof constraint rows
4730 : auto receive_dof_constraints =
4731 59604 : [this,
4732 : & pushed_ids_rhss_to_me,
4733 : & pushed_keys_vals_to_me]
4734 265381 : ()
4735 : {
4736 103967 : for (const auto & [pid, ids_rhss] : pushed_ids_rhss_to_me)
4737 : {
4738 38955 : const auto & keys_vals = pushed_keys_vals_to_me[pid];
4739 :
4740 684 : libmesh_assert_equal_to
4741 : (ids_rhss.size(), keys_vals.size());
4742 :
4743 : // Add the dof constraints that I've been sent
4744 291272 : for (auto i : index_range(ids_rhss))
4745 : {
4746 273333 : dof_id_type constrained = ids_rhss[i].first;
4747 :
4748 : // If we don't already have a constraint for this dof,
4749 : // add the one we were sent
4750 89657 : if (!this->is_constrained_dof(constrained))
4751 : {
4752 181036 : DofConstraintRow & row = _dof_constraints[constrained];
4753 702501 : for (auto & key_val : keys_vals[i])
4754 : {
4755 48033 : libmesh_assert_less(key_val.first, this->n_dofs());
4756 501769 : row[key_val.first] = key_val.second;
4757 : }
4758 200732 : if (ids_rhss[i].second != Number(0))
4759 7896 : _primal_constraint_values[constrained] =
4760 108 : ids_rhss[i].second;
4761 : else
4762 19588 : _primal_constraint_values.erase(constrained);
4763 : }
4764 : }
4765 : }
4766 66364 : };
4767 :
4768 32506 : receive_dof_constraints();
4769 :
4770 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4771 : // Collect the node constraints to push to each processor
4772 219544 : for (auto & i : _node_constraints)
4773 : {
4774 216840 : const Node * constrained = i.first;
4775 :
4776 325260 : if (constrained->processor_id() != this->processor_id())
4777 24811 : continue;
4778 :
4779 83609 : NodeConstraintRow & row = i.second.first;
4780 743150 : for (auto & j : row)
4781 : {
4782 575932 : const Node * constraining = j.first;
4783 :
4784 897246 : if (constraining->processor_id() != this->processor_id() &&
4785 28686 : constraining->processor_id() != constrained->processor_id())
4786 28686 : pushed_node_ids[constraining->processor_id()].insert(constrained->id());
4787 : }
4788 : }
4789 :
4790 : // Pack the node constraint rows and rhss to push
4791 : std::map<processor_id_type,
4792 : std::vector<std::vector<std::pair<dof_id_type,Real>>>>
4793 2704 : pushed_node_keys_vals, pushed_node_keys_vals_to_me;
4794 : std::map<processor_id_type, std::vector<std::pair<dof_id_type, Point>>>
4795 2704 : pushed_node_ids_offsets, pushed_node_ids_offsets_to_me;
4796 2704 : std::map<processor_id_type, std::vector<const Node *>> pushed_node_vecs;
4797 :
4798 4134 : for (const auto & [pid, pid_ids]: pushed_node_ids)
4799 : {
4800 715 : const std::size_t ids_size = pid_ids.size();
4801 : std::vector<std::vector<std::pair<dof_id_type,Real>>> &
4802 1430 : keys_vals = pushed_node_keys_vals[pid];
4803 : std::vector<std::pair<dof_id_type, Point>> &
4804 1430 : ids_offsets = pushed_node_ids_offsets[pid];
4805 1430 : keys_vals.resize(ids_size);
4806 1430 : ids_offsets.resize(ids_size);
4807 1430 : std::set<Node *> nodes;
4808 :
4809 : std::size_t push_i;
4810 715 : std::set<dof_id_type>::const_iterator it;
4811 17455 : for (push_i = 0, it = pid_ids.begin();
4812 18170 : it != pid_ids.end(); ++push_i, ++it)
4813 : {
4814 16740 : Node * constrained = mesh.node_ptr(*it);
4815 :
4816 16740 : if (constrained->processor_id() != pid)
4817 8370 : nodes.insert(constrained);
4818 :
4819 16740 : NodeConstraintRow & row = _node_constraints[constrained].first;
4820 8370 : std::size_t row_size = row.size();
4821 25110 : keys_vals[push_i].reserve(row_size);
4822 80026 : for (const auto & j : row)
4823 : {
4824 63286 : Node * constraining = const_cast<Node *>(j.first);
4825 :
4826 96617 : keys_vals[push_i].emplace_back(constraining->id(), j.second);
4827 :
4828 63286 : if (constraining->processor_id() != pid)
4829 15612 : nodes.insert(constraining);
4830 : }
4831 :
4832 16740 : ids_offsets[push_i].first = *it;
4833 16740 : ids_offsets[push_i].second = _node_constraints[constrained].second;
4834 : }
4835 :
4836 1430 : if (!mesh.is_serial())
4837 : {
4838 0 : auto & pid_nodes = pushed_node_vecs[pid];
4839 0 : pid_nodes.assign(nodes.begin(), nodes.end());
4840 : }
4841 : }
4842 :
4843 : auto node_ids_offsets_action_functor =
4844 : [& pushed_node_ids_offsets_to_me]
4845 : (processor_id_type pid,
4846 715 : const std::vector<std::pair<dof_id_type, Point>> & data)
4847 : {
4848 1430 : pushed_node_ids_offsets_to_me[pid] = data;
4849 3419 : };
4850 :
4851 : auto node_keys_vals_action_functor =
4852 : [& pushed_node_keys_vals_to_me]
4853 : (processor_id_type pid,
4854 715 : const std::vector<std::vector<std::pair<dof_id_type, Real>>> & data)
4855 : {
4856 1430 : pushed_node_keys_vals_to_me[pid] = data;
4857 3419 : };
4858 :
4859 : // Trade pushed node constraint rows
4860 : Parallel::push_parallel_vector_data
4861 2704 : (this->comm(), pushed_node_ids_offsets, node_ids_offsets_action_functor);
4862 : Parallel::push_parallel_vector_data
4863 2704 : (this->comm(), pushed_node_keys_vals, node_keys_vals_action_functor);
4864 :
4865 : // Constraining nodes might not even exist on our subset of a
4866 : // distributed mesh, so let's make them exist.
4867 :
4868 : // Node unpack() now automatically adds them to the context mesh
4869 0 : auto null_node_functor = [](processor_id_type, const std::vector<const Node *> &){};
4870 :
4871 2704 : if (!mesh.is_serial())
4872 : Parallel::push_parallel_packed_range
4873 36 : (this->comm(), pushed_node_vecs, &mesh, null_node_functor);
4874 :
4875 4134 : for (const auto & [pid, ids_offsets] : pushed_node_ids_offsets_to_me)
4876 : {
4877 1430 : const auto & keys_vals = pushed_node_keys_vals_to_me[pid];
4878 :
4879 715 : libmesh_assert_equal_to
4880 : (ids_offsets.size(), keys_vals.size());
4881 :
4882 : // Add the node constraints that I've been sent
4883 18170 : for (auto i : index_range(ids_offsets))
4884 : {
4885 16740 : dof_id_type constrained_id = ids_offsets[i].first;
4886 :
4887 : // If we don't already have a constraint for this node,
4888 : // add the one we were sent
4889 16740 : const Node * constrained = mesh.node_ptr(constrained_id);
4890 11766 : if (!this->is_constrained_node(constrained))
4891 : {
4892 9948 : NodeConstraintRow & row = _node_constraints[constrained].first;
4893 53308 : for (auto & key_val : keys_vals[i])
4894 : {
4895 38386 : const Node * key_node = mesh.node_ptr(key_val.first);
4896 38386 : row[key_node] = key_val.second;
4897 : }
4898 9948 : _node_constraints[constrained].second =
4899 9948 : ids_offsets[i].second;
4900 : }
4901 : }
4902 : }
4903 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4904 :
4905 : // Next we need to push constraints to processors which don't own
4906 : // the constrained dof, don't own the constraining dof, but own an
4907 : // element supporting the constraining dof.
4908 : //
4909 : // We need to be able to quickly look up constrained dof ids by what
4910 : // constrains them, so that we can handle the case where we see a
4911 : // foreign element containing one of our constraining DoF ids and we
4912 : // need to push that constraint.
4913 : //
4914 : // Getting distributed adaptive sparsity patterns right is hard.
4915 :
4916 : typedef std::map<dof_id_type, std::set<dof_id_type>> DofConstrainsMap;
4917 2704 : DofConstrainsMap dof_id_constrains;
4918 :
4919 1158354 : for (const auto & [constrained, row] : _dof_constraints)
4920 : {
4921 2996751 : for (const auto & j : row)
4922 : {
4923 1870903 : const dof_id_type constraining = j.first;
4924 :
4925 195306 : dof_id_type constraining_proc_id = 0;
4926 7361501 : while (constraining >= _end_df[constraining_proc_id])
4927 5196092 : constraining_proc_id++;
4928 :
4929 2066209 : if (constraining_proc_id == this->processor_id())
4930 1384268 : dof_id_constrains[constraining].insert(constrained);
4931 : }
4932 : }
4933 :
4934 : // Loop over all foreign elements, find any supporting our
4935 : // constrained dof indices.
4936 1352 : pushed_ids.clear();
4937 :
4938 322490 : for (const auto & elem : as_range(mesh.active_not_local_elements_begin(),
4939 12537800 : mesh.active_not_local_elements_end()))
4940 : {
4941 517660 : std::vector<dof_id_type> my_dof_indices;
4942 6334655 : this->dof_indices (elem, my_dof_indices);
4943 :
4944 51556245 : for (const auto & dof : my_dof_indices)
4945 : {
4946 45221590 : if (auto dcmi = dof_id_constrains.find(dof);
4947 1754347 : dcmi != dof_id_constrains.end())
4948 : {
4949 705213 : for (const auto & constrained : dcmi->second)
4950 : {
4951 23562 : dof_id_type the_constrained_proc_id = 0;
4952 2631383 : while (constrained >= _end_df[the_constrained_proc_id])
4953 2049196 : the_constrained_proc_id++;
4954 :
4955 549216 : const processor_id_type elemproc = elem->processor_id();
4956 549216 : if (elemproc != the_constrained_proc_id)
4957 351240 : pushed_ids[elemproc].insert(constrained);
4958 : }
4959 : }
4960 : }
4961 29802 : }
4962 :
4963 1352 : pushed_ids_rhss.clear();
4964 1352 : pushed_ids_rhss_to_me.clear();
4965 1352 : pushed_keys_vals.clear();
4966 1352 : pushed_keys_vals_to_me.clear();
4967 :
4968 32506 : gather_ids();
4969 :
4970 : // Trade pushed dof constraint rows
4971 : Parallel::push_parallel_vector_data
4972 32506 : (this->comm(), pushed_ids_rhss, ids_rhss_action_functor);
4973 : Parallel::push_parallel_vector_data
4974 32506 : (this->comm(), pushed_keys_vals, keys_vals_action_functor);
4975 :
4976 32506 : receive_dof_constraints();
4977 :
4978 : // Finally, we need to handle the case of remote dof coupling. If a
4979 : // processor's element is coupled to a ghost element, then the
4980 : // processor needs to know about all constraints which affect the
4981 : // dofs on that ghost element, so we'll have to query the ghost
4982 : // element's owner.
4983 :
4984 2704 : GhostingFunctor::map_type elements_to_couple;
4985 2704 : DofMap::CouplingMatricesSet temporary_coupling_matrices;
4986 :
4987 : this->merge_ghost_functor_outputs
4988 98870 : (elements_to_couple,
4989 : temporary_coupling_matrices,
4990 62308 : this->coupling_functors_begin(),
4991 33858 : this->coupling_functors_end(),
4992 65012 : mesh.active_local_elements_begin(),
4993 65012 : mesh.active_local_elements_end(),
4994 : this->processor_id());
4995 :
4996 : // Each ghost-coupled element's owner should get a request for its dofs
4997 2704 : std::set<dof_id_type> requested_dofs;
4998 :
4999 143849 : for (const auto & pr : elements_to_couple)
5000 : {
5001 111343 : const Elem * elem = pr.first;
5002 :
5003 : // FIXME - optimize for the non-fully-coupled case?
5004 20234 : std::vector<dof_id_type> element_dofs;
5005 111343 : this->dof_indices(elem, element_dofs);
5006 :
5007 774634 : for (auto dof : element_dofs)
5008 627031 : requested_dofs.insert(dof);
5009 : }
5010 :
5011 32506 : this->gather_constraints(mesh, requested_dofs, false);
5012 29802 : }
5013 :
5014 :
5015 65081 : void DofMap::gather_constraints (MeshBase & /*mesh*/,
5016 : std::set<dof_id_type> & unexpanded_dofs,
5017 : bool /*look_for_constrainees*/)
5018 : {
5019 : typedef std::set<dof_id_type> DoF_RCSet;
5020 :
5021 : // If we have heterogeneous adjoint constraints we need to
5022 : // communicate those too.
5023 : const unsigned int max_qoi_num =
5024 65081 : _adjoint_constraint_values.empty() ?
5025 1348 : 0 : _adjoint_constraint_values.rbegin()->first+1;
5026 :
5027 : // We have to keep recursing while the unexpanded set is
5028 : // nonempty on *any* processor
5029 65081 : bool unexpanded_set_nonempty = !unexpanded_dofs.empty();
5030 65081 : this->comm().max(unexpanded_set_nonempty);
5031 :
5032 104898 : while (unexpanded_set_nonempty)
5033 : {
5034 : // Let's make sure we don't lose sync in this loop.
5035 1612 : parallel_object_only();
5036 :
5037 : // Request sets
5038 3224 : DoF_RCSet dof_request_set;
5039 :
5040 : // Request sets to send to each processor
5041 : std::map<processor_id_type, std::vector<dof_id_type>>
5042 3224 : requested_dof_ids;
5043 :
5044 : // And the sizes of each
5045 : std::map<processor_id_type, dof_id_type>
5046 3224 : dof_ids_on_proc;
5047 :
5048 : // Fill (and thereby sort and uniq!) the main request sets
5049 1561863 : for (const auto & unexpanded_dof : unexpanded_dofs)
5050 : {
5051 : // If we were asked for a DoF and we don't already have a
5052 : // constraint for it, then we need to check for one.
5053 1522046 : if (auto pos = _dof_constraints.find(unexpanded_dof);
5054 124549 : pos == _dof_constraints.end())
5055 : {
5056 488970 : if (!this->local_index(unexpanded_dof) &&
5057 19720 : !_dof_constraints.count(unexpanded_dof) )
5058 368687 : dof_request_set.insert(unexpanded_dof);
5059 : }
5060 : // If we were asked for a DoF and we already have a
5061 : // constraint for it, then we need to check if the
5062 : // constraint is recursive.
5063 : else
5064 : {
5065 100870 : const DofConstraintRow & row = pos->second;
5066 2749834 : for (const auto & j : row)
5067 : {
5068 1697038 : const dof_id_type constraining_dof = j.first;
5069 :
5070 : // If it's non-local and we haven't already got a
5071 : // constraint for it, we might need to ask for one
5072 1566446 : if (!this->local_index(constraining_dof) &&
5073 32786 : !_dof_constraints.count(constraining_dof))
5074 379819 : dof_request_set.insert(constraining_dof);
5075 : }
5076 : }
5077 : }
5078 :
5079 : // Clear the unexpanded constraint set; we're about to expand it
5080 1612 : unexpanded_dofs.clear();
5081 :
5082 : // Count requests by processor
5083 39817 : processor_id_type proc_id = 0;
5084 584272 : for (const auto & i : dof_request_set)
5085 : {
5086 667564 : while (i >= _end_df[proc_id])
5087 94945 : proc_id++;
5088 544455 : dof_ids_on_proc[proc_id]++;
5089 : }
5090 :
5091 86147 : for (auto & pair : dof_ids_on_proc)
5092 : {
5093 46330 : requested_dof_ids[pair.first].reserve(pair.second);
5094 : }
5095 :
5096 : // Prepare each processor's request set
5097 39817 : proc_id = 0;
5098 584272 : for (const auto & i : dof_request_set)
5099 : {
5100 667564 : while (i >= _end_df[proc_id])
5101 94945 : proc_id++;
5102 544455 : requested_dof_ids[proc_id].push_back(i);
5103 : }
5104 :
5105 : typedef std::vector<std::pair<dof_id_type, Real>> row_datum;
5106 :
5107 : typedef std::vector<Number> rhss_datum;
5108 :
5109 : auto row_gather_functor =
5110 44466 : [this]
5111 : (processor_id_type,
5112 : const std::vector<dof_id_type> & ids,
5113 547225 : std::vector<row_datum> & data)
5114 : {
5115 : // Fill those requests
5116 1864 : const std::size_t query_size = ids.size();
5117 :
5118 46330 : data.resize(query_size);
5119 590785 : for (std::size_t i=0; i != query_size; ++i)
5120 : {
5121 572218 : dof_id_type constrained = ids[i];
5122 27763 : if (_dof_constraints.count(constrained))
5123 : {
5124 17192 : DofConstraintRow & row = _dof_constraints[constrained];
5125 906 : std::size_t row_size = row.size();
5126 18098 : data[i].reserve(row_size);
5127 41588 : for (const auto & j : row)
5128 : {
5129 25518 : data[i].push_back(j);
5130 :
5131 : // We should never have an invalid constraining
5132 : // dof id
5133 1122 : libmesh_assert(j.first != DofObject::invalid_id);
5134 :
5135 : // We should never have a 0 constraint
5136 : // coefficient; that's implicit via sparse
5137 : // constraint storage
5138 : //
5139 : // But we can't easily control how users add
5140 : // constraints, so we can't safely assert that
5141 : // we're being efficient here.
5142 : //
5143 : // libmesh_assert(j.second);
5144 : }
5145 : }
5146 : else
5147 : {
5148 : // We have to distinguish "constraint with no
5149 : // constraining dofs" (e.g. due to Dirichlet
5150 : // constraint equations) from "no constraint".
5151 : // We'll use invalid_id for the latter.
5152 554120 : data[i].emplace_back(DofObject::invalid_id, Real(0));
5153 : }
5154 : }
5155 84535 : };
5156 :
5157 : auto rhss_gather_functor =
5158 44466 : [this,
5159 : max_qoi_num]
5160 : (processor_id_type,
5161 : const std::vector<dof_id_type> & ids,
5162 548105 : std::vector<rhss_datum> & data)
5163 : {
5164 : // Fill those requests
5165 1864 : const std::size_t query_size = ids.size();
5166 :
5167 46330 : data.resize(query_size);
5168 590785 : for (std::size_t i=0; i != query_size; ++i)
5169 : {
5170 544455 : dof_id_type constrained = ids[i];
5171 55526 : data[i].clear();
5172 27763 : if (_dof_constraints.count(constrained))
5173 : {
5174 : DofConstraintValueMap::const_iterator rhsit =
5175 906 : _primal_constraint_values.find(constrained);
5176 906 : data[i].push_back
5177 17369 : ((rhsit == _primal_constraint_values.end()) ?
5178 177 : 0 : rhsit->second);
5179 :
5180 17192 : for (unsigned int q = 0; q != max_qoi_num; ++q)
5181 : {
5182 : AdjointDofConstraintValues::const_iterator adjoint_map_it =
5183 0 : _adjoint_constraint_values.find(q);
5184 :
5185 0 : if (adjoint_map_it == _adjoint_constraint_values.end())
5186 : {
5187 0 : data[i].push_back(0);
5188 0 : continue;
5189 : }
5190 :
5191 : const DofConstraintValueMap & constraint_map =
5192 0 : adjoint_map_it->second;
5193 :
5194 : DofConstraintValueMap::const_iterator adj_rhsit =
5195 0 : constraint_map.find(constrained);
5196 0 : data[i].push_back
5197 0 : ((adj_rhsit == constraint_map.end()) ?
5198 0 : 0 : adj_rhsit->second);
5199 : }
5200 : }
5201 : }
5202 47942 : };
5203 :
5204 : auto row_action_functor =
5205 44466 : [this,
5206 : & unexpanded_dofs]
5207 : (processor_id_type,
5208 : const std::vector<dof_id_type> & ids,
5209 33690 : const std::vector<row_datum> & data)
5210 : {
5211 : // Add any new constraint rows we've found
5212 1864 : const std::size_t query_size = ids.size();
5213 :
5214 590785 : for (std::size_t i=0; i != query_size; ++i)
5215 : {
5216 544455 : const dof_id_type constrained = ids[i];
5217 :
5218 : // An empty row is an constraint with an empty row; for
5219 : // no constraint we use a "no row" placeholder
5220 572218 : if (data[i].empty())
5221 : {
5222 1842 : DofConstraintRow & row = _dof_constraints[constrained];
5223 182 : row.clear();
5224 : }
5225 542613 : else if (data[i][0].first != DofObject::invalid_id)
5226 : {
5227 15350 : DofConstraintRow & row = _dof_constraints[constrained];
5228 724 : row.clear();
5229 40470 : for (auto & pair : data[i])
5230 : {
5231 1122 : libmesh_assert_less(pair.first, this->n_dofs());
5232 24396 : row[pair.first] = pair.second;
5233 : }
5234 :
5235 : // And prepare to check for more recursive constraints
5236 14626 : unexpanded_dofs.insert(constrained);
5237 : }
5238 : }
5239 84535 : };
5240 :
5241 : auto rhss_action_functor =
5242 44466 : [this,
5243 : max_qoi_num]
5244 : (processor_id_type,
5245 : const std::vector<dof_id_type> & ids,
5246 18522 : const std::vector<rhss_datum> & data)
5247 : {
5248 : // Add rhs data for any new constraint rows we've found
5249 1864 : const std::size_t query_size = ids.size();
5250 :
5251 590785 : for (std::size_t i=0; i != query_size; ++i)
5252 : {
5253 572218 : if (!data[i].empty())
5254 : {
5255 17192 : dof_id_type constrained = ids[i];
5256 17192 : if (data[i][0] != Number(0))
5257 519 : _primal_constraint_values[constrained] = data[i][0];
5258 : else
5259 895 : _primal_constraint_values.erase(constrained);
5260 :
5261 17192 : for (unsigned int q = 0; q != max_qoi_num; ++q)
5262 : {
5263 : AdjointDofConstraintValues::iterator adjoint_map_it =
5264 0 : _adjoint_constraint_values.find(q);
5265 :
5266 0 : if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
5267 0 : data[i][q+1] == Number(0))
5268 0 : continue;
5269 :
5270 0 : if (adjoint_map_it == _adjoint_constraint_values.end())
5271 0 : adjoint_map_it = _adjoint_constraint_values.emplace
5272 0 : (q, DofConstraintValueMap()).first;
5273 :
5274 : DofConstraintValueMap & constraint_map =
5275 0 : adjoint_map_it->second;
5276 :
5277 0 : if (data[i][q+1] != Number(0))
5278 0 : constraint_map[constrained] =
5279 0 : data[i][q+1];
5280 : else
5281 0 : constraint_map.erase(constrained);
5282 : }
5283 : }
5284 : }
5285 :
5286 84535 : };
5287 :
5288 : // Now request constraint rows from other processors
5289 1612 : row_datum * row_ex = nullptr;
5290 : Parallel::pull_parallel_vector_data
5291 39817 : (this->comm(), requested_dof_ids, row_gather_functor,
5292 : row_action_functor, row_ex);
5293 :
5294 : // And request constraint right hand sides from other procesors
5295 1612 : rhss_datum * rhs_ex = nullptr;
5296 : Parallel::pull_parallel_vector_data
5297 39817 : (this->comm(), requested_dof_ids, rhss_gather_functor,
5298 : rhss_action_functor, rhs_ex);
5299 :
5300 : // We have to keep recursing while the unexpanded set is
5301 : // nonempty on *any* processor
5302 39817 : unexpanded_set_nonempty = !unexpanded_dofs.empty();
5303 39817 : this->comm().max(unexpanded_set_nonempty);
5304 : }
5305 65081 : }
5306 :
5307 335540 : void DofMap::add_constraints_to_send_list (const MeshBase & mesh)
5308 : {
5309 : // This function must be run on all processors at once
5310 9618 : parallel_object_only();
5311 :
5312 : // Return immediately if there's nothing to gather
5313 345158 : if (this->n_processors() == 1)
5314 331814 : return;
5315 :
5316 : // We might get to return immediately if none of the processors
5317 : // found any constraints
5318 326291 : unsigned int has_constraints = !_dof_constraints.empty();
5319 316673 : this->comm().max(has_constraints);
5320 326291 : if (!has_constraints)
5321 8554 : return;
5322 :
5323 1094518 : auto add_row = [this](const DofConstraintRow & constraint_row) {
5324 2318063 : for (const auto & j : constraint_row)
5325 : {
5326 1402949 : dof_id_type constraint_dependency = j.first;
5327 :
5328 : // No point in adding one of our own dofs to the send_list
5329 1251048 : if (this->local_index(constraint_dependency))
5330 1172163 : continue;
5331 :
5332 230786 : _send_list.push_back(constraint_dependency);
5333 : }
5334 223510 : };
5335 :
5336 : // We usually only need dependencies of our own constrained dofs
5337 1240501 : for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
5338 1208571 : if (this->local_index(constrained_dof))
5339 908889 : add_row(constraint_row);
5340 :
5341 : // If we only need constraint DoFs constraining DoFs which are
5342 : // algebraically local, we're done.
5343 31930 : if (!this->has_static_condensation() && !_need_ghost_constraints)
5344 956 : return;
5345 :
5346 : // If we use StaticCondensation, though, we may need constraint DoFs
5347 : // constraining DoFs which are not local (they're on someone else's
5348 : // node) but which are supported on local elements. Let's get those
5349 : // too if we have to.
5350 : //
5351 : // Kokkos-MOOSE also requires that, because it cannot use VecSetValues
5352 : // and MatSetValues which can dynamically cache remote entries. Instead,
5353 : // it allocates ghost entries in advance based on libMesh send_list and
5354 : // accumulates to those, and does the assembly at once. When ghost DOFs
5355 : // are constrained, it has to add residual/Jacobian contributions to the
5356 : // dependencies of ghost DOFs locally as well.
5357 : //
5358 : // We'll potentially be hitting the same constrained DoFs from
5359 : // multiple directions.
5360 303192 : for (auto & elem : mesh.active_local_element_ptr_range())
5361 : {
5362 30036 : std::vector<dof_id_type> di;
5363 162942 : this->dof_indices (elem, di);
5364 2350646 : for (const auto & dof_id : di)
5365 2187704 : if (!this->local_index(dof_id))
5366 135715 : if (auto pos = _dof_constraints.find(dof_id);
5367 4042 : pos != _dof_constraints.end())
5368 6225 : add_row(pos->second);
5369 3510 : }
5370 : }
5371 :
5372 :
5373 :
5374 : #endif // LIBMESH_ENABLE_CONSTRAINTS
5375 :
5376 :
5377 : #ifdef LIBMESH_ENABLE_AMR
5378 :
5379 576 : void DofMap::constrain_p_dofs (unsigned int var,
5380 : const Elem * elem,
5381 : unsigned int s,
5382 : unsigned int p)
5383 : {
5384 : // We're constraining dofs on elem which correspond to p refinement
5385 : // levels above p - this only makes sense if elem's p refinement
5386 : // level is above p.
5387 48 : libmesh_assert_greater (elem->p_level(), p);
5388 48 : libmesh_assert_less (s, elem->n_sides());
5389 :
5390 96 : const unsigned int sys_num = this->sys_number();
5391 576 : FEType fe_type = this->variable_type(var);
5392 :
5393 576 : const unsigned int n_nodes = elem->n_nodes();
5394 5760 : for (unsigned int n = 0; n != n_nodes; ++n)
5395 5184 : if (elem->is_node_on_side(n, s))
5396 : {
5397 144 : const Node & node = elem->node_ref(n);
5398 : const unsigned int low_nc =
5399 1728 : FEInterface::n_dofs_at_node (fe_type, p, elem, n);
5400 : const unsigned int high_nc =
5401 1728 : FEInterface::n_dofs_at_node (fe_type, elem, n);
5402 :
5403 : // since we may be running this method concurrently
5404 : // on multiple threads we need to acquire a lock
5405 : // before modifying the _dof_constraints object.
5406 288 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
5407 :
5408 1728 : if (elem->is_vertex(n))
5409 : {
5410 : // Add "this is zero" constraint rows for high p vertex
5411 : // dofs
5412 1152 : for (unsigned int i = low_nc; i != high_nc; ++i)
5413 : {
5414 0 : _dof_constraints[node.dof_number(sys_num,var,i)].clear();
5415 0 : _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
5416 : }
5417 : }
5418 : else
5419 : {
5420 576 : const unsigned int total_dofs = node.n_comp(sys_num, var);
5421 48 : libmesh_assert_greater_equal (total_dofs, high_nc);
5422 : // Add "this is zero" constraint rows for high p
5423 : // non-vertex dofs, which are numbered in reverse
5424 1152 : for (unsigned int j = low_nc; j != high_nc; ++j)
5425 : {
5426 576 : const unsigned int i = total_dofs - j - 1;
5427 576 : _dof_constraints[node.dof_number(sys_num,var,i)].clear();
5428 576 : _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
5429 : }
5430 : }
5431 : }
5432 576 : }
5433 :
5434 : #endif // LIBMESH_ENABLE_AMR
5435 :
5436 :
5437 : #ifdef LIBMESH_ENABLE_DIRICHLET
5438 11174 : void DofMap::add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary)
5439 : {
5440 21716 : _dirichlet_boundaries->push_back(std::make_unique<DirichletBoundary>(dirichlet_boundary));
5441 11174 : }
5442 :
5443 :
5444 2311 : void DofMap::add_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
5445 : unsigned int qoi_index)
5446 : {
5447 : unsigned int old_size = cast_int<unsigned int>
5448 132 : (_adjoint_dirichlet_boundaries.size());
5449 3502 : for (unsigned int i = old_size; i <= qoi_index; ++i)
5450 2348 : _adjoint_dirichlet_boundaries.push_back(std::make_unique<DirichletBoundaries>());
5451 :
5452 : // Make copy of DirichletBoundary, owned by _adjoint_dirichlet_boundaries
5453 2311 : _adjoint_dirichlet_boundaries[qoi_index]->push_back
5454 4490 : (std::make_unique<DirichletBoundary>(dirichlet_boundary));
5455 2311 : }
5456 :
5457 :
5458 132476 : bool DofMap::has_adjoint_dirichlet_boundaries(unsigned int q) const
5459 : {
5460 136260 : if (_adjoint_dirichlet_boundaries.size() > q)
5461 130772 : return true;
5462 :
5463 48 : return false;
5464 : }
5465 :
5466 :
5467 : const DirichletBoundaries *
5468 0 : DofMap::get_adjoint_dirichlet_boundaries(unsigned int q) const
5469 : {
5470 0 : libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),q);
5471 0 : return _adjoint_dirichlet_boundaries[q].get();
5472 : }
5473 :
5474 :
5475 : DirichletBoundaries *
5476 0 : DofMap::get_adjoint_dirichlet_boundaries(unsigned int q)
5477 : {
5478 : unsigned int old_size = cast_int<unsigned int>
5479 0 : (_adjoint_dirichlet_boundaries.size());
5480 0 : for (unsigned int i = old_size; i <= q; ++i)
5481 0 : _adjoint_dirichlet_boundaries.push_back(std::make_unique<DirichletBoundaries>());
5482 :
5483 0 : return _adjoint_dirichlet_boundaries[q].get();
5484 : }
5485 :
5486 :
5487 0 : void DofMap::remove_dirichlet_boundary (const DirichletBoundary & boundary_to_remove)
5488 : {
5489 : // Find a boundary condition matching the one to be removed
5490 0 : auto lam = [&boundary_to_remove](const auto & bdy)
5491 0 : {return bdy->b == boundary_to_remove.b && bdy->variables == boundary_to_remove.variables;};
5492 :
5493 0 : auto it = std::find_if(_dirichlet_boundaries->begin(), _dirichlet_boundaries->end(), lam);
5494 :
5495 : // Assert it was actually found and remove it from the vector
5496 0 : libmesh_assert (it != _dirichlet_boundaries->end());
5497 0 : _dirichlet_boundaries->erase(it);
5498 0 : }
5499 :
5500 :
5501 0 : void DofMap::remove_adjoint_dirichlet_boundary (const DirichletBoundary & boundary_to_remove,
5502 : unsigned int qoi_index)
5503 : {
5504 0 : libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),
5505 : qoi_index);
5506 :
5507 0 : auto lam = [&boundary_to_remove](const auto & bdy)
5508 0 : {return bdy->b == boundary_to_remove.b && bdy->variables == boundary_to_remove.variables;};
5509 :
5510 0 : auto it = std::find_if(_adjoint_dirichlet_boundaries[qoi_index]->begin(),
5511 0 : _adjoint_dirichlet_boundaries[qoi_index]->end(),
5512 0 : lam);
5513 :
5514 : // Assert it was actually found and remove it from the vector
5515 0 : libmesh_assert (it != _adjoint_dirichlet_boundaries[qoi_index]->end());
5516 0 : _adjoint_dirichlet_boundaries[qoi_index]->erase(it);
5517 0 : }
5518 :
5519 :
5520 24611 : void DofMap::check_dirichlet_bcid_consistency (const MeshBase & mesh,
5521 : const DirichletBoundary & boundary) const
5522 : {
5523 : const std::set<boundary_id_type>& mesh_side_bcids =
5524 696 : mesh.get_boundary_info().get_boundary_ids();
5525 : const std::set<boundary_id_type>& mesh_edge_bcids =
5526 696 : mesh.get_boundary_info().get_edge_boundary_ids();
5527 : const std::set<boundary_id_type>& mesh_node_bcids =
5528 696 : mesh.get_boundary_info().get_node_boundary_ids();
5529 696 : const std::set<boundary_id_type>& dbc_bcids = boundary.b;
5530 :
5531 : // DirichletBoundary id sets should be consistent across all ranks
5532 696 : libmesh_assert(mesh.comm().verify(dbc_bcids.size()));
5533 :
5534 71951 : for (const auto & bc_id : dbc_bcids)
5535 : {
5536 : // DirichletBoundary id sets should be consistent across all ranks
5537 1340 : libmesh_assert(mesh.comm().verify(bc_id));
5538 :
5539 19968 : bool found_bcid = (mesh_side_bcids.find(bc_id) != mesh_side_bcids.end() ||
5540 65968 : mesh_edge_bcids.find(bc_id) != mesh_edge_bcids.end() ||
5541 46006 : mesh_node_bcids.find(bc_id) != mesh_node_bcids.end());
5542 :
5543 : // On a distributed mesh, boundary id sets may *not* be
5544 : // consistent across all ranks, since not all ranks see all
5545 : // boundaries
5546 47340 : mesh.comm().max(found_bcid);
5547 :
5548 47340 : libmesh_error_msg_if(!found_bcid,
5549 : "Could not find Dirichlet boundary id " << bc_id << " in mesh!");
5550 : }
5551 24611 : }
5552 :
5553 : #endif // LIBMESH_ENABLE_DIRICHLET
5554 :
5555 :
5556 : #ifdef LIBMESH_ENABLE_PERIODIC
5557 :
5558 513 : void DofMap::add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary)
5559 : {
5560 535 : auto inverse_boundary = periodic_boundary.clone(PeriodicBoundaryBase::INVERSE);
5561 513 : this->add_periodic_boundary(periodic_boundary, *inverse_boundary);
5562 513 : }
5563 :
5564 :
5565 :
5566 513 : void DofMap::add_periodic_boundary (const PeriodicBoundaryBase & boundary,
5567 : const PeriodicBoundaryBase & inverse_boundary)
5568 : {
5569 22 : libmesh_assert_equal_to (boundary.myboundary, inverse_boundary.pairedboundary);
5570 22 : libmesh_assert_equal_to (boundary.pairedboundary, inverse_boundary.myboundary);
5571 22 : libmesh_assert(boundary.get_variables() == inverse_boundary.get_variables());
5572 :
5573 : // See if we already have a periodic boundary associated myboundary...
5574 : PeriodicBoundaryBase * existing_boundary =
5575 513 : _periodic_boundaries->boundary(boundary.myboundary);
5576 :
5577 513 : if (!existing_boundary)
5578 : {
5579 : // ...if not, clone the inputs and add them to the
5580 : // PeriodicBoundaries object.
5581 : // These will be cleaned up automatically in the
5582 : // _periodic_boundaries destructor.
5583 535 : _periodic_boundaries->emplace(boundary.myboundary, boundary.clone());
5584 535 : _periodic_boundaries->emplace(inverse_boundary.myboundary, inverse_boundary.clone());
5585 : }
5586 : else
5587 : {
5588 : // ...otherwise, merge this object's variable IDs with the
5589 : // existing boundary object's.
5590 0 : existing_boundary->merge(boundary);
5591 :
5592 : // Do the same merging process for the inverse boundary. The
5593 : // inverse had better already exist!
5594 : PeriodicBoundaryBase * existing_inverse_boundary =
5595 0 : _periodic_boundaries->boundary(boundary.pairedboundary);
5596 0 : libmesh_assert(existing_inverse_boundary);
5597 0 : existing_inverse_boundary->merge(inverse_boundary);
5598 :
5599 : // If we had to merge different *types* of boundaries then
5600 : // something likely has gone wrong.
5601 : #ifdef LIBMESH_HAVE_RTTI
5602 : // typeid needs to be given references, not just pointers, to
5603 : // return a derived class name.
5604 0 : libmesh_assert(typeid(boundary) == typeid(*existing_boundary));
5605 0 : libmesh_assert(typeid(inverse_boundary) == typeid(*existing_inverse_boundary));
5606 : #endif
5607 : }
5608 513 : }
5609 :
5610 :
5611 : #endif
5612 :
5613 :
5614 : } // namespace libMesh
|
