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