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