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