Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // C++ includes
21 : #include <vector>
22 : #include <numeric> // std::iota
23 :
24 : // Local includes
25 : #include "libmesh/libmesh_config.h"
26 :
27 : #ifdef LIBMESH_HAVE_METAPHYSICL
28 :
29 : // With quad precision we need the shim function declarations to
30 : // precede the MetaPhysicL use of them
31 : #include "libmesh/libmesh_common.h"
32 : #include "libmesh/compare_types.h"
33 :
34 : // Template specialization declarations in here need to *precede* code
35 : // using them.
36 : #include "metaphysicl/dynamicsparsenumberarray_decl.h"
37 :
38 : using MetaPhysicL::DynamicSparseNumberArray;
39 :
40 : namespace libMesh
41 : {
42 : // From the perspective of libMesh gradient vectors, a DSNA is a
43 : // scalar component
44 : template <typename T, typename IndexType>
45 : struct ScalarTraits<MetaPhysicL::DynamicSparseNumberArray<T,IndexType> >
46 : {
47 : static const bool value = true;
48 : };
49 :
50 : // And although MetaPhysicL knows how to combine DSNA with something
51 : // else, we need to teach libMesh too.
52 : template <typename T, typename IndexType, typename T2>
53 : struct CompareTypes<MetaPhysicL::DynamicSparseNumberArray<T,IndexType>, T2>
54 : {
55 : typedef typename
56 : MetaPhysicL::DynamicSparseNumberArray
57 : <typename CompareTypes<T,T2>::supertype,IndexType> supertype;
58 : };
59 :
60 : template <typename T> struct TypeToSend;
61 :
62 : template <typename T, typename IndexType>
63 : struct TypeToSend<MetaPhysicL::DynamicSparseNumberArray<T,IndexType>> {
64 : typedef std::vector<std::pair<IndexType,T>> type;
65 : };
66 :
67 : template <typename T, typename IndexType>
68 : const std::vector<std::pair<IndexType,T>>
69 0 : convert_to_send(MetaPhysicL::DynamicSparseNumberArray<T,IndexType> & in)
70 : {
71 0 : const std::size_t in_size = in.size();
72 0 : std::vector<std::pair<IndexType,T>> returnval(in_size);
73 :
74 0 : for (std::size_t i=0; i != in_size; ++i)
75 : {
76 0 : returnval[i].first = in.raw_index(i);
77 0 : returnval[i].second = in.raw_at(i);
78 : }
79 0 : return returnval;
80 : }
81 :
82 : template <typename SendT, typename T, typename IndexType>
83 0 : void convert_from_receive (SendT & received,
84 : MetaPhysicL::DynamicSparseNumberArray<T,IndexType> & converted)
85 : {
86 0 : const std::size_t received_size = received.size();
87 0 : converted.resize(received_size);
88 0 : for (std::size_t i=0; i != received_size; ++i)
89 : {
90 0 : converted.raw_index(i) = received[i].first;
91 0 : converted.raw_at(i) = received[i].second;
92 : }
93 0 : }
94 :
95 : }
96 :
97 :
98 : #endif
99 :
100 : #include "libmesh/boundary_info.h"
101 : #include "libmesh/dense_matrix.h"
102 : #include "libmesh/dense_vector.h"
103 : #include "libmesh/dof_map.h"
104 : #include "libmesh/elem.h"
105 : #include "libmesh/fe_base.h"
106 : #include "libmesh/fe_interface.h"
107 : #include "libmesh/generic_projector.h"
108 : #include "libmesh/int_range.h"
109 : #include "libmesh/libmesh_logging.h"
110 : #include "libmesh/linear_solver.h"
111 : #include "libmesh/mesh_base.h"
112 : #include "libmesh/numeric_vector.h"
113 : #include "libmesh/quadrature.h"
114 : #include "libmesh/sparse_matrix.h"
115 : #include "libmesh/system.h"
116 : #include "libmesh/threads.h"
117 : #include "libmesh/wrapped_function.h"
118 : #include "libmesh/wrapped_functor.h"
119 : #include "libmesh/fe_interface.h"
120 :
121 :
122 :
123 : #ifdef LIBMESH_HAVE_METAPHYSICL
124 : // Include MetaPhysicL definitions finally
125 : #include "metaphysicl/dynamicsparsenumberarray.h"
126 :
127 : // And make sure we instantiate the methods we'll need to use on them.
128 : #include "libmesh/dense_matrix_impl.h"
129 :
130 : namespace libMesh {
131 : typedef DynamicSparseNumberArray<Real, dof_id_type> DSNAN;
132 :
133 : template LIBMESH_EXPORT void
134 : DenseMatrix<Real>::cholesky_solve(const DenseVector<DSNAN> &,
135 : DenseVector<DSNAN> &);
136 : template LIBMESH_EXPORT void
137 : DenseMatrix<Real>::_cholesky_back_substitute(const DenseVector<DSNAN> &,
138 : DenseVector<DSNAN> &) const;
139 : }
140 : #endif
141 :
142 :
143 :
144 : namespace libMesh
145 : {
146 :
147 : // ------------------------------------------------------------
148 : // Helper class definitions
149 :
150 : #ifdef LIBMESH_ENABLE_AMR
151 :
152 : /**
153 : * This class builds the send_list of old dof indices
154 : * whose coefficients are needed to perform a projection.
155 : * This may be executed in parallel on multiple threads.
156 : * The end result is a \p send_list vector which is
157 : * unsorted and may contain duplicate elements.
158 : * The \p unique() method can be used to sort and
159 : * create a unique list.
160 : */
161 :
162 53978 : class BuildProjectionList
163 : {
164 : private:
165 : const System & system;
166 :
167 : public:
168 1996 : BuildProjectionList (const System & system_in) :
169 52940 : system(system_in),
170 3992 : send_list()
171 1996 : {}
172 :
173 3882 : BuildProjectionList (BuildProjectionList & other, Threads::split) :
174 3882 : system(other.system),
175 2844 : send_list()
176 1422 : {}
177 :
178 : void unique();
179 : void operator()(const ConstElemRange & range);
180 : void join (const BuildProjectionList & other);
181 : std::vector<dof_id_type> send_list;
182 : };
183 :
184 : #endif // LIBMESH_ENABLE_AMR
185 :
186 :
187 : /**
188 : * This class implements projecting an arbitrary
189 : * boundary function to the current mesh. This
190 : * may be executed in parallel on multiple threads.
191 : */
192 67 : class BoundaryProjectSolution
193 : {
194 : private:
195 : const std::set<boundary_id_type> & b;
196 : const std::vector<unsigned int> & variables;
197 : const System & system;
198 : std::unique_ptr<FunctionBase<Number>> f;
199 : std::unique_ptr<FunctionBase<Gradient>> g;
200 : const Parameters & parameters;
201 : NumericVector<Number> & new_vector;
202 :
203 : public:
204 71 : BoundaryProjectSolution (const std::set<boundary_id_type> & b_in,
205 : const std::vector<unsigned int> & variables_in,
206 : const System & system_in,
207 : FunctionBase<Number> * f_in,
208 : FunctionBase<Gradient> * g_in,
209 : const Parameters & parameters_in,
210 71 : NumericVector<Number> & new_v_in) :
211 67 : b(b_in),
212 67 : variables(variables_in),
213 67 : system(system_in),
214 71 : f(f_in ? f_in->clone() : std::unique_ptr<FunctionBase<Number>>()),
215 67 : g(g_in ? g_in->clone() : std::unique_ptr<FunctionBase<Gradient>>()),
216 67 : parameters(parameters_in),
217 73 : new_vector(new_v_in)
218 : {
219 2 : libmesh_assert(f.get());
220 71 : f->init();
221 71 : if (g.get())
222 0 : g->init();
223 71 : }
224 :
225 : BoundaryProjectSolution (const BoundaryProjectSolution & in) :
226 : b(in.b),
227 : variables(in.variables),
228 : system(in.system),
229 : f(in.f.get() ? in.f->clone() : std::unique_ptr<FunctionBase<Number>>()),
230 : g(in.g.get() ? in.g->clone() : std::unique_ptr<FunctionBase<Gradient>>()),
231 : parameters(in.parameters),
232 : new_vector(in.new_vector)
233 : {
234 : libmesh_assert(f.get());
235 : f->init();
236 : if (g.get())
237 : g->init();
238 : }
239 :
240 : void operator()(const ConstElemRange & range) const;
241 : };
242 :
243 :
244 :
245 : // ------------------------------------------------------------
246 : // System implementation
247 57665 : void System::project_vector (NumericVector<Number> & vector,
248 : int is_adjoint,
249 : std::optional<ConstElemRange> active_local_range,
250 : std::optional<std::vector<unsigned int>> variable_numbers) const
251 : {
252 : // Create a copy of the vector, which currently
253 : // contains the old data.
254 : std::unique_ptr<NumericVector<Number>>
255 57665 : old_vector (vector.clone());
256 :
257 : // Project the old vector to the new vector
258 113334 : this->project_vector (*old_vector, vector, is_adjoint, active_local_range, variable_numbers);
259 57665 : }
260 :
261 :
262 : /**
263 : * This method projects the vector
264 : * via L2 projections or nodal
265 : * interpolations on each element.
266 : */
267 57665 : void System::project_vector (const NumericVector<Number> & old_v,
268 : NumericVector<Number> & new_v,
269 : int is_adjoint,
270 : std::optional<ConstElemRange> active_local_range,
271 : std::optional<std::vector<unsigned int>> variable_numbers) const
272 : {
273 3992 : LOG_SCOPE ("project_vector(old,new)", "System");
274 :
275 : /**
276 : * This method projects a solution from an old mesh to a current, refined
277 : * mesh. The input vector \p old_v gives the solution on the
278 : * old mesh, while the \p new_v gives the solution (to be computed)
279 : * on the new mesh.
280 : */
281 57665 : new_v.clear();
282 :
283 : #ifdef LIBMESH_ENABLE_AMR
284 :
285 : // Resize the new vector and get a serial version.
286 1996 : NumericVector<Number> * new_vector_ptr = nullptr;
287 55669 : std::unique_ptr<NumericVector<Number>> new_vector_built;
288 : NumericVector<Number> * local_old_vector;
289 55669 : std::unique_ptr<NumericVector<Number>> local_old_vector_built;
290 1996 : const NumericVector<Number> * old_vector_ptr = nullptr;
291 :
292 57665 : if (!active_local_range)
293 : {
294 : active_local_range.emplace
295 115330 : (this->get_mesh().active_local_elements_begin(),
296 119322 : this->get_mesh().active_local_elements_end());
297 : }
298 :
299 : // If the old vector was uniprocessor, make the new
300 : // vector uniprocessor
301 57665 : if (old_v.type() == SERIAL)
302 : {
303 733 : new_v.init (this->n_dofs(), false, SERIAL);
304 0 : new_vector_ptr = &new_v;
305 0 : old_vector_ptr = &old_v;
306 : }
307 :
308 : // Otherwise it is a parallel, distributed vector, which
309 : // we need to localize.
310 56932 : else if (old_v.type() == PARALLEL)
311 : {
312 : // Build a send list for efficient localization
313 1256 : BuildProjectionList projection_list(*this);
314 37124 : Threads::parallel_reduce (active_local_range.value(),
315 : projection_list);
316 :
317 : // Create a sorted, unique send_list
318 37124 : projection_list.unique();
319 :
320 37124 : new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
321 71736 : new_vector_built = NumericVector<Number>::build(this->comm());
322 71736 : local_old_vector_built = NumericVector<Number>::build(this->comm());
323 1256 : new_vector_ptr = new_vector_built.get();
324 1256 : local_old_vector = local_old_vector_built.get();
325 37124 : new_vector_ptr->init(this->n_dofs(), this->n_local_dofs(),
326 1256 : this->get_dof_map().get_send_list(), false,
327 2512 : GHOSTED);
328 37124 : local_old_vector->init(old_v.size(), old_v.local_size(),
329 2512 : projection_list.send_list, false, GHOSTED);
330 37124 : old_v.localize(*local_old_vector, projection_list.send_list);
331 37124 : local_old_vector->close();
332 1256 : old_vector_ptr = local_old_vector;
333 : }
334 19808 : else if (old_v.type() == GHOSTED)
335 : {
336 : // Build a send list for efficient localization
337 740 : BuildProjectionList projection_list(*this);
338 19808 : Threads::parallel_reduce (active_local_range.value(),
339 : projection_list);
340 :
341 : // Create a sorted, unique send_list
342 19808 : projection_list.unique();
343 :
344 19808 : new_v.init (this->n_dofs(), this->n_local_dofs(),
345 1480 : this->get_dof_map().get_send_list(), false, GHOSTED);
346 :
347 38136 : local_old_vector_built = NumericVector<Number>::build(this->comm());
348 740 : new_vector_ptr = &new_v;
349 740 : local_old_vector = local_old_vector_built.get();
350 19808 : local_old_vector->init(old_v.size(), old_v.local_size(),
351 1480 : projection_list.send_list, false, GHOSTED);
352 19808 : old_v.localize(*local_old_vector, projection_list.send_list);
353 19808 : local_old_vector->close();
354 740 : old_vector_ptr = local_old_vector;
355 : }
356 : else // unknown old_v.type()
357 0 : libmesh_error_msg("ERROR: Unknown old_v.type() == " << old_v.type());
358 :
359 : // Note that the above will have zeroed the new_vector.
360 : // Just to be sure, assert that new_vector_ptr and old_vector_ptr
361 : // were successfully set before trying to deref them.
362 1996 : libmesh_assert(new_vector_ptr);
363 1996 : libmesh_assert(old_vector_ptr);
364 :
365 1996 : NumericVector<Number> & new_vector = *new_vector_ptr;
366 1996 : const NumericVector<Number> & old_vector = *old_vector_ptr;
367 :
368 57665 : const unsigned int n_variables = this->n_vars();
369 :
370 57665 : if (n_variables)
371 : {
372 3776 : std::vector<unsigned int> vars;
373 53883 : if (variable_numbers)
374 : {
375 0 : vars = *variable_numbers;
376 0 : for (auto v : vars)
377 0 : if (v >= n_variables)
378 0 : libmesh_error_msg("ERROR: variable number " << v <<
379 : " out of range for system with " <<
380 : n_variables << " variables.");
381 : }
382 : else
383 : {
384 53883 : vars.resize(n_variables);
385 1888 : std::iota(vars.begin(), vars.end(), 0);
386 : }
387 :
388 3776 : std::vector<unsigned int> regular_vars, vector_vars;
389 118004 : for (auto var : vars)
390 : {
391 64121 : if (FEInterface::field_type(this->variable_type(var)) == TYPE_SCALAR)
392 58938 : regular_vars.push_back(var);
393 : else
394 5183 : vector_vars.push_back(var);
395 : }
396 :
397 1888 : VectorSetAction<Number> setter(new_vector);
398 :
399 53883 : if (!regular_vars.empty())
400 : {
401 : // Use a typedef to make the calling sequence for parallel_for() a bit more readable
402 : typedef
403 : GenericProjector<OldSolutionValue<Number, &FEMContext::point_value>,
404 : OldSolutionValue<Gradient, &FEMContext::point_gradient>,
405 : Number, VectorSetAction<Number>> FEMProjector;
406 :
407 : OldSolutionValue<Number, &FEMContext::point_value>
408 3484 : f(*this, old_vector, ®ular_vars);
409 : OldSolutionValue<Gradient, &FEMContext::point_gradient>
410 3484 : g(*this, old_vector, ®ular_vars);
411 :
412 52184 : FEMProjector projector(*this, f, &g, setter, regular_vars);
413 48700 : projector.project(active_local_range.value());
414 45216 : }
415 :
416 53883 : if (!vector_vars.empty())
417 : {
418 : typedef
419 : GenericProjector<OldSolutionValue<Gradient, &FEMContext::point_value>,
420 : OldSolutionValue<Tensor, &FEMContext::point_gradient>,
421 : Gradient, VectorSetAction<Number>> FEMVectorProjector;
422 :
423 292 : OldSolutionValue<Gradient, &FEMContext::point_value> f_vector(*this, old_vector, &vector_vars);
424 292 : OldSolutionValue<Tensor, &FEMContext::point_gradient> g_vector(*this, old_vector, &vector_vars);
425 :
426 5475 : FEMVectorProjector vector_projector(*this, f_vector, &g_vector, setter, vector_vars);
427 5183 : vector_projector.project(active_local_range.value());
428 4891 : }
429 :
430 : // Copy the SCALAR dofs from old_vector to new_vector
431 : // Note: We assume that all SCALAR dofs are on the
432 : // processor with highest ID
433 55771 : if (this->processor_id() == (this->n_processors()-1))
434 : {
435 944 : const DofMap & dof_map = this->get_dof_map();
436 20336 : for (auto var : vars)
437 10981 : if (this->variable(var).type().family == SCALAR)
438 : {
439 : // We can just map SCALAR dofs directly across
440 0 : std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
441 0 : dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
442 0 : dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
443 0 : for (auto i : index_range(new_SCALAR_indices))
444 0 : new_vector.set(new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]));
445 : }
446 : }
447 : }
448 :
449 57665 : new_vector.close();
450 :
451 : // If the old vector was serial, we probably need to send our values
452 : // to other processors
453 : //
454 : // FIXME: I'm not sure how to make a NumericVector do that without
455 : // creating a temporary parallel vector to use localize! - RHS
456 57665 : if (old_v.type() == SERIAL)
457 : {
458 733 : std::unique_ptr<NumericVector<Number>> dist_v = NumericVector<Number>::build(this->comm());
459 733 : dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
460 733 : dist_v->close();
461 :
462 625559 : for (auto i : make_range(dist_v->size()))
463 624093 : if (new_vector(i) != 0.0)
464 489127 : dist_v->set(i, new_vector(i));
465 :
466 733 : dist_v->close();
467 :
468 733 : dist_v->localize (new_v, this->get_dof_map().get_send_list());
469 733 : new_v.close();
470 733 : }
471 : // If the old vector was parallel, we need to update it
472 : // and free the localized copies
473 56932 : else if (old_v.type() == PARALLEL)
474 : {
475 : // We may have to set dof values that this processor doesn't
476 : // own in certain special cases, like LAGRANGE FIRST or
477 : // HERMITE THIRD elements on second-order meshes?
478 37124 : new_v = new_vector;
479 37124 : new_v.close();
480 : }
481 :
482 :
483 : // Apply constraints only if we we are asked to
484 57665 : if(this->project_with_constraints)
485 : {
486 53831 : if (is_adjoint == -1)
487 : {
488 50646 : this->get_dof_map().enforce_constraints_exactly(*this, &new_v);
489 : }
490 3185 : else if (is_adjoint >= 0)
491 : {
492 3185 : this->get_dof_map().enforce_adjoint_constraints_exactly(new_v,
493 : is_adjoint);
494 : }
495 : }
496 : #else
497 :
498 : // AMR is disabled: simply copy the vector
499 : new_v = old_v;
500 :
501 : libmesh_ignore(is_adjoint, active_local_range, variable_numbers);
502 :
503 : #endif // #ifdef LIBMESH_ENABLE_AMR
504 57665 : }
505 :
506 :
507 : #ifdef LIBMESH_ENABLE_AMR
508 : #ifdef LIBMESH_HAVE_METAPHYSICL
509 :
510 : template <typename Output>
511 : class DSNAOutput
512 : {
513 : public:
514 : typedef DynamicSparseNumberArray<Output, dof_id_type> type;
515 : };
516 :
517 : template <typename InnerOutput>
518 : class DSNAOutput<VectorValue<InnerOutput> >
519 : {
520 : public:
521 : typedef VectorValue<DynamicSparseNumberArray<InnerOutput, dof_id_type> > type;
522 : };
523 :
524 : /**
525 : * The OldSolutionCoefs input functor class can be used with
526 : * GenericProjector to read solution transfer coefficients on a
527 : * just-refined-and-coarsened mesh.
528 : *
529 : * \author Roy H. Stogner
530 : * \date 2017
531 : */
532 :
533 : template <typename Output,
534 : void (FEMContext::*point_output) (unsigned int,
535 : const Point &,
536 : Output &,
537 : const Real) const>
538 1206 : class OldSolutionCoefs : public OldSolutionBase<Output, point_output>
539 : {
540 : public:
541 : typedef typename DSNAOutput<Output>::type DSNA;
542 : typedef DSNA ValuePushType;
543 : typedef DSNA FunctorValue;
544 :
545 1260 : OldSolutionCoefs(const libMesh::System & sys_in,
546 : const std::vector<unsigned int> * vars) :
547 648 : OldSolutionBase<Output, point_output>(sys_in, vars)
548 : {
549 36 : this->old_context.set_algebraic_type(FEMContext::OLD_DOFS_ONLY);
550 36 : }
551 :
552 3306 : OldSolutionCoefs(const OldSolutionCoefs & in) :
553 3448 : OldSolutionBase<Output, point_output>(in.sys, in.old_context.active_vars())
554 : {
555 142 : this->old_context.set_algebraic_type(FEMContext::OLD_DOFS_ONLY);
556 142 : }
557 :
558 : DSNA eval_at_node (const FEMContext & c,
559 : unsigned int i,
560 : unsigned int elem_dim,
561 : const Node & n,
562 : bool extra_hanging_dofs,
563 : Real /* time */ = 0.);
564 :
565 : DSNA eval_at_point(const FEMContext & c,
566 : unsigned int i,
567 : const Point & p,
568 : Real time,
569 : bool skip_context_check);
570 :
571 0 : void eval_mixed_derivatives (const FEMContext & libmesh_dbg_var(c),
572 : unsigned int i,
573 : unsigned int dim,
574 : const Node & n,
575 : std::vector<DSNA> & derivs)
576 : {
577 0 : LOG_SCOPE ("eval_mixed_derivatives", "OldSolutionCoefs");
578 :
579 : // This should only be called on vertices
580 0 : libmesh_assert_less(c.get_elem().get_node_index(&n),
581 : c.get_elem().n_vertices());
582 :
583 : // Handle offset from non-scalar components in previous variables
584 0 : libmesh_assert_less(i, this->component_to_var.size());
585 0 : unsigned int var = this->component_to_var[i];
586 :
587 : // We have 1 mixed derivative in 2D, 4 in 3D
588 0 : const unsigned int n_mixed = (dim-1) * (dim-1);
589 0 : derivs.resize(n_mixed);
590 :
591 : // Be sure to handle cases where the variable wasn't defined on
592 : // this node (e.g. due to changing subdomain support)
593 0 : const DofObject * old_dof_object = n.get_old_dof_object();
594 0 : if (old_dof_object &&
595 0 : old_dof_object->n_vars(this->sys.number()) &&
596 0 : old_dof_object->n_comp(this->sys.number(), var))
597 : {
598 0 : const dof_id_type first_old_id =
599 0 : old_dof_object->dof_number(this->sys.number(), var, dim);
600 0 : std::vector<dof_id_type> old_ids(n_mixed);
601 0 : std::iota(old_ids.begin(), old_ids.end(), first_old_id);
602 :
603 0 : for (auto d_i : index_range(derivs))
604 : {
605 0 : derivs[d_i].resize(1);
606 0 : derivs[d_i].raw_at(0) = 1;
607 0 : derivs[d_i].raw_index(0) = old_ids[d_i];
608 : }
609 : }
610 : else
611 : {
612 0 : std::fill(derivs.begin(), derivs.end(), 0);
613 : }
614 0 : }
615 :
616 :
617 0 : void eval_old_dofs (const Elem & elem,
618 : unsigned int node_num,
619 : unsigned int var_num,
620 : std::vector<dof_id_type> & indices,
621 : std::vector<DSNA> & values)
622 : {
623 0 : LOG_SCOPE ("eval_old_dofs(node)", "OldSolutionCoefs");
624 :
625 : // We may be reusing a std::vector here, but the following
626 : // dof_indices call appends without first clearing.
627 0 : indices.clear();
628 :
629 0 : this->sys.get_dof_map().dof_indices(elem, node_num, indices, var_num);
630 :
631 0 : std::vector<dof_id_type> old_indices;
632 :
633 0 : this->sys.get_dof_map().old_dof_indices(elem, node_num, old_indices, var_num);
634 :
635 0 : libmesh_assert_equal_to (old_indices.size(), indices.size());
636 :
637 0 : values.resize(old_indices.size());
638 :
639 0 : for (auto i : index_range(values))
640 : {
641 0 : values[i].resize(1);
642 0 : values[i].raw_at(0) = 1;
643 0 : values[i].raw_index(0) = old_indices[i];
644 : }
645 0 : }
646 :
647 :
648 0 : void eval_old_dofs (const Elem & elem,
649 : const FEType & fe_type,
650 : unsigned int sys_num,
651 : unsigned int var_num,
652 : std::vector<dof_id_type> & indices,
653 : std::vector<DSNA> & values)
654 : {
655 0 : LOG_SCOPE ("eval_old_dofs(elem)", "OldSolutionCoefs");
656 :
657 : // We're only to be asked for old dofs on elements that can copy
658 : // them through DO_NOTHING or through refinement.
659 0 : const Elem & old_elem =
660 0 : (elem.refinement_flag() == Elem::JUST_REFINED) ?
661 0 : *elem.parent() : elem;
662 :
663 : // If there are any element-based DOF numbers, get them
664 0 : const unsigned int nc =
665 0 : FEInterface::n_dofs_per_elem(fe_type, &elem);
666 :
667 0 : std::vector<dof_id_type> old_dof_indices(nc);
668 0 : indices.resize(nc);
669 :
670 : // We should never have fewer dofs than necessary on an
671 : // element unless we're getting indices on a parent element,
672 : // and we should never need those indices
673 0 : if (nc != 0)
674 : {
675 0 : const DofObject & old_dof_object = old_elem.get_old_dof_object_ref();
676 :
677 0 : const auto [vg, vig] =
678 0 : elem.var_to_vg_and_offset(sys_num,var_num);
679 :
680 0 : const unsigned int n_comp = elem.n_comp_group(sys_num,vg);
681 0 : libmesh_assert_greater(elem.n_systems(), sys_num);
682 0 : libmesh_assert_greater_equal(n_comp, nc);
683 :
684 0 : for (unsigned int i=0; i<nc; i++)
685 : {
686 0 : const dof_id_type d_old =
687 0 : old_dof_object.dof_number(sys_num, vg, vig, i, n_comp);
688 0 : const dof_id_type d_new =
689 0 : elem.dof_number(sys_num, vg, vig, i, n_comp);
690 0 : libmesh_assert_not_equal_to (d_old, DofObject::invalid_id);
691 0 : libmesh_assert_not_equal_to (d_new, DofObject::invalid_id);
692 :
693 0 : old_dof_indices[i] = d_old;
694 0 : indices[i] = d_new;
695 : }
696 : }
697 :
698 0 : values.resize(old_dof_indices.size());
699 :
700 0 : for (auto i : index_range(values))
701 : {
702 0 : values[i].resize(1);
703 0 : values[i].raw_at(0) = 1;
704 0 : values[i].raw_index(0) = old_dof_indices[i];
705 : }
706 0 : }
707 : };
708 :
709 :
710 :
711 : template<>
712 : inline
713 : DynamicSparseNumberArray<Real, dof_id_type>
714 119439 : OldSolutionCoefs<Real, &FEMContext::point_value>::
715 : eval_at_point(const FEMContext & c,
716 : unsigned int i,
717 : const Point & p,
718 : Real /* time */,
719 : bool skip_context_check)
720 : {
721 20704 : LOG_SCOPE ("eval_at_point()", "OldSolutionCoefs");
722 :
723 119439 : if (!skip_context_check)
724 119439 : if (!this->check_old_context(c, p))
725 0 : return 0;
726 :
727 : // Get finite element object
728 10352 : FEGenericBase<Real> * fe = nullptr;
729 : this->old_context.get_element_fe<Real>
730 10352 : (i, fe, this->old_context.get_elem_dim());
731 :
732 : // Build a FE for calculating phi(p)
733 : FEGenericBase<Real> * fe_new =
734 119439 : this->old_context.build_new_fe(fe, p);
735 :
736 : // Get the values and global indices of the shape functions
737 10352 : const std::vector<std::vector<Real> > & phi = fe_new->get_phi();
738 : const std::vector<dof_id_type> & dof_indices =
739 10352 : this->old_context.get_dof_indices(i);
740 :
741 20704 : const std::size_t n_dofs = phi.size();
742 10352 : libmesh_assert_equal_to(n_dofs, dof_indices.size());
743 :
744 20704 : DynamicSparseNumberArray<Real, dof_id_type> returnval;
745 109087 : returnval.resize(n_dofs);
746 :
747 1116532 : for (auto j : index_range(phi))
748 : {
749 997093 : returnval.raw_at(j) = phi[j][0];
750 1170881 : returnval.raw_index(j) = dof_indices[j];
751 : }
752 :
753 10352 : return returnval;
754 : }
755 :
756 :
757 :
758 : template<>
759 : inline
760 : VectorValue<DynamicSparseNumberArray<Real, dof_id_type> >
761 0 : OldSolutionCoefs<RealGradient, &FEMContext::point_gradient>::
762 : eval_at_point(const FEMContext & c,
763 : unsigned int i,
764 : const Point & p,
765 : Real /* time */,
766 : bool skip_context_check)
767 : {
768 0 : LOG_SCOPE ("eval_at_point()", "OldSolutionCoefs");
769 :
770 0 : if (!skip_context_check)
771 0 : if (!this->check_old_context(c, p))
772 0 : return 0;
773 :
774 : // Get finite element object
775 0 : FEGenericBase<Real> * fe = nullptr;
776 : this->old_context.get_element_fe<Real>
777 0 : (i, fe, this->old_context.get_elem_dim());
778 :
779 : // Build a FE for calculating phi(p)
780 : FEGenericBase<Real> * fe_new =
781 0 : this->old_context.build_new_fe(fe, p);
782 :
783 : // Get the values and global indices of the shape functions
784 0 : const std::vector<std::vector<RealGradient> > & dphi = fe_new->get_dphi();
785 : const std::vector<dof_id_type> & dof_indices =
786 0 : this->old_context.get_dof_indices(i);
787 :
788 0 : const std::size_t n_dofs = dphi.size();
789 0 : libmesh_assert_equal_to(n_dofs, dof_indices.size());
790 :
791 0 : VectorValue<DynamicSparseNumberArray<Real, dof_id_type> > returnval;
792 :
793 0 : for (unsigned int d = 0; d != LIBMESH_DIM; ++d)
794 0 : returnval(d).resize(n_dofs);
795 :
796 0 : for (auto j : index_range(dphi))
797 0 : for (int d = 0; d != LIBMESH_DIM; ++d)
798 : {
799 0 : returnval(d).raw_at(j) = dphi[j][0](d);
800 0 : returnval(d).raw_index(j) = dof_indices[j];
801 : }
802 :
803 0 : return returnval;
804 : }
805 :
806 :
807 : template<>
808 : inline
809 : DynamicSparseNumberArray<Real, dof_id_type>
810 62759 : OldSolutionCoefs<Real, &FEMContext::point_value>::
811 : eval_at_node(const FEMContext & c,
812 : unsigned int i,
813 : unsigned int /* elem_dim */,
814 : const Node & n,
815 : bool extra_hanging_dofs,
816 : Real /* time */)
817 : {
818 10206 : LOG_SCOPE ("Real eval_at_node()", "OldSolutionCoefs");
819 :
820 : // Optimize for the common case, where this node was part of the
821 : // old solution.
822 : //
823 : // Be sure to handle cases where the variable wasn't defined on
824 : // this node (due to changing subdomain support) or where the
825 : // variable has no components on this node (due to Elem order
826 : // exceeding FE order) or where the old_dof_object dofs might
827 : // correspond to non-vertex dofs (due to extra_hanging_dofs and
828 : // refinement)
829 :
830 10206 : const Elem::RefinementState flag = c.get_elem().refinement_flag();
831 :
832 5103 : const DofObject * old_dof_object = n.get_old_dof_object();
833 66707 : if (old_dof_object &&
834 66639 : (!extra_hanging_dofs ||
835 56569 : flag == Elem::JUST_COARSENED ||
836 61604 : flag == Elem::DO_NOTHING) &&
837 185967 : old_dof_object->n_vars(sys.number()) &&
838 61604 : old_dof_object->n_comp(sys.number(), i))
839 : {
840 7658 : DynamicSparseNumberArray<Real, dof_id_type> returnval;
841 : const dof_id_type old_id =
842 46874 : old_dof_object->dof_number(sys.number(), i, 0);
843 43045 : returnval.resize(1);
844 46874 : returnval.raw_at(0) = 1;
845 46874 : returnval.raw_index(0) = old_id;
846 3829 : return returnval;
847 : }
848 :
849 15885 : return this->eval_at_point(c, i, n, 0, false);
850 : }
851 :
852 :
853 :
854 : template<>
855 : inline
856 : VectorValue<DynamicSparseNumberArray<Real, dof_id_type> >
857 0 : OldSolutionCoefs<RealGradient, &FEMContext::point_gradient>::
858 : eval_at_node(const FEMContext & c,
859 : unsigned int i,
860 : unsigned int elem_dim,
861 : const Node & n,
862 : bool extra_hanging_dofs,
863 : Real /* time */)
864 : {
865 0 : LOG_SCOPE ("RealGradient eval_at_node()", "OldSolutionCoefs");
866 :
867 : // Optimize for the common case, where this node was part of the
868 : // old solution.
869 : //
870 : // Be sure to handle cases where the variable wasn't defined on
871 : // this node (due to changing subdomain support) or where the
872 : // variable has no components on this node (due to Elem order
873 : // exceeding FE order) or where the old_dof_object dofs might
874 : // correspond to non-vertex dofs (due to extra_hanging_dofs and
875 : // refinement)
876 :
877 0 : const Elem::RefinementState flag = c.get_elem().refinement_flag();
878 :
879 0 : const DofObject * old_dof_object = n.get_old_dof_object();
880 0 : if (old_dof_object &&
881 0 : (!extra_hanging_dofs ||
882 0 : flag == Elem::JUST_COARSENED ||
883 0 : flag == Elem::DO_NOTHING) &&
884 0 : old_dof_object->n_vars(sys.number()) &&
885 0 : old_dof_object->n_comp(sys.number(), i))
886 : {
887 0 : VectorValue<DynamicSparseNumberArray<Real, dof_id_type> > g;
888 0 : for (unsigned int d = 0; d != elem_dim; ++d)
889 : {
890 : const dof_id_type old_id =
891 0 : old_dof_object->dof_number(sys.number(), i, d+1);
892 0 : g(d).resize(1);
893 0 : g(d).raw_at(0) = 1;
894 0 : g(d).raw_index(0) = old_id;
895 : }
896 0 : return g;
897 : }
898 :
899 0 : return this->eval_at_point(c, i, n, 0, false);
900 : }
901 :
902 :
903 :
904 : /**
905 : * The MatrixFillAction output functor class can be used with
906 : * GenericProjector to write solution transfer coefficients into a
907 : * sparse matrix.
908 : *
909 : * \author Roy H. Stogner
910 : * \date 2017
911 : */
912 : template <typename ValIn, typename ValOut>
913 : class MatrixFillAction
914 : {
915 : public:
916 : typedef DynamicSparseNumberArray<ValIn, dof_id_type> InsertInput;
917 : private:
918 : SparseMatrix<ValOut> & target_matrix;
919 :
920 : public:
921 630 : MatrixFillAction(SparseMatrix<ValOut> & target_mat) :
922 630 : target_matrix(target_mat) {}
923 :
924 166313 : void insert(dof_id_type id,
925 : const DynamicSparseNumberArray<ValIn, dof_id_type> & val)
926 : {
927 : // Lock the target matrix since it is shared among threads.
928 : {
929 28362 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
930 :
931 14181 : const std::size_t dnsa_size = val.size();
932 1210280 : for (unsigned int j = 0; j != dnsa_size; ++j)
933 : {
934 1043967 : const dof_id_type dof_j = val.raw_index(j);
935 1043967 : const ValIn dof_val = val.raw_at(j);
936 1043967 : target_matrix.set(id, dof_j, dof_val);
937 : }
938 : }
939 166313 : }
940 :
941 :
942 0 : void insert(const std::vector<dof_id_type> & dof_indices,
943 : const std::vector<DynamicSparseNumberArray<ValIn, dof_id_type> > & Ue)
944 : {
945 : const numeric_index_type
946 0 : begin_dof = target_matrix.row_start(),
947 0 : end_dof = target_matrix.row_stop();
948 :
949 0 : unsigned int size = Ue.size();
950 :
951 0 : libmesh_assert_equal_to(size, dof_indices.size());
952 :
953 : // Lock the target matrix since it is shared among threads.
954 : {
955 0 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
956 :
957 0 : for (unsigned int i = 0; i != size; ++i)
958 : {
959 0 : const dof_id_type dof_i = dof_indices[i];
960 0 : if ((dof_i >= begin_dof) && (dof_i < end_dof))
961 : {
962 0 : const DynamicSparseNumberArray<ValIn,dof_id_type> & dnsa = Ue[i];
963 0 : const std::size_t dnsa_size = dnsa.size();
964 0 : for (unsigned int j = 0; j != dnsa_size; ++j)
965 : {
966 0 : const dof_id_type dof_j = dnsa.raw_index(j);
967 0 : const ValIn dof_val = dnsa.raw_at(j);
968 0 : target_matrix.set(dof_i, dof_j, dof_val);
969 : }
970 : }
971 : }
972 : }
973 0 : }
974 : };
975 :
976 :
977 :
978 : /**
979 : * This method creates a projection matrix which corresponds to the
980 : * operation of project_vector between old and new solution spaces.
981 : */
982 630 : void System::projection_matrix (SparseMatrix<Number> & proj_mat) const
983 : {
984 36 : LOG_SCOPE ("projection_matrix()", "System");
985 :
986 630 : const unsigned int n_variables = this->n_vars();
987 :
988 630 : if (n_variables)
989 : {
990 : ConstElemRange active_local_elem_range
991 1260 : (this->get_mesh().active_local_elements_begin(),
992 1278 : this->get_mesh().active_local_elements_end());
993 :
994 648 : std::vector<unsigned int> vars(n_variables);
995 18 : std::iota(vars.begin(), vars.end(), 0);
996 :
997 : // Use a typedef to make the calling sequence for parallel_for() a bit more readable
998 : typedef OldSolutionCoefs<Real, &FEMContext::point_value> OldSolutionValueCoefs;
999 : typedef OldSolutionCoefs<RealGradient, &FEMContext::point_gradient> OldSolutionGradientCoefs;
1000 :
1001 : typedef
1002 : GenericProjector<OldSolutionValueCoefs,
1003 : OldSolutionGradientCoefs,
1004 : DynamicSparseNumberArray<Real,dof_id_type>,
1005 : MatrixFillAction<Real, Number> > ProjMatFiller;
1006 :
1007 36 : OldSolutionValueCoefs f(*this, &vars);
1008 36 : OldSolutionGradientCoefs g(*this, &vars);
1009 18 : MatrixFillAction<Real, Number> setter(proj_mat);
1010 :
1011 666 : ProjMatFiller mat_filler(*this, f, &g, setter, vars);
1012 630 : mat_filler.project(active_local_elem_range);
1013 :
1014 : // Set the SCALAR dof transfer entries too.
1015 : // Note: We assume that all SCALAR dofs are on the
1016 : // processor with highest ID
1017 648 : if (this->processor_id() == (this->n_processors()-1))
1018 : {
1019 9 : const DofMap & dof_map = this->get_dof_map();
1020 376 : for (auto var : make_range(this->n_vars()))
1021 187 : if (this->variable(var).type().family == SCALAR)
1022 : {
1023 : // We can just map SCALAR dofs directly across
1024 0 : std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
1025 0 : dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
1026 0 : dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
1027 : const unsigned int new_n_dofs =
1028 0 : cast_int<unsigned int>(new_SCALAR_indices.size());
1029 :
1030 0 : for (unsigned int i=0; i<new_n_dofs; i++)
1031 : {
1032 0 : proj_mat.set( new_SCALAR_indices[i],
1033 0 : old_SCALAR_indices[i], 1);
1034 : }
1035 : }
1036 : }
1037 594 : }
1038 630 : }
1039 : #endif // LIBMESH_HAVE_METAPHYSICL
1040 : #endif // LIBMESH_ENABLE_AMR
1041 :
1042 :
1043 :
1044 : /**
1045 : * This method projects an arbitrary function onto the solution via L2
1046 : * projections and nodal interpolations on each element.
1047 : */
1048 199154 : void System::project_solution (ValueFunctionPointer fptr,
1049 : GradientFunctionPointer gptr,
1050 : const Parameters & function_parameters,
1051 : std::optional<ConstElemRange> active_local_range,
1052 : std::optional<std::vector<unsigned int>> variable_numbers) const
1053 : {
1054 11220 : WrappedFunction<Number> f(*this, fptr, &function_parameters);
1055 5610 : WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
1056 392698 : this->project_solution(&f, &g, active_local_range, variable_numbers);
1057 199154 : }
1058 :
1059 :
1060 : /**
1061 : * This method projects an arbitrary function onto the solution via L2
1062 : * projections and nodal interpolations on each element.
1063 : */
1064 199727 : void System::project_solution (FunctionBase<Number> * f,
1065 : FunctionBase<Gradient> * g,
1066 : std::optional<ConstElemRange> active_local_range,
1067 : std::optional<std::vector<unsigned int>> variable_numbers) const
1068 : {
1069 393895 : this->project_vector(*solution, f, g, /*is_adjoint=*/-1, active_local_range, variable_numbers);
1070 :
1071 205355 : solution->localize(*current_local_solution, _dof_map->get_send_list());
1072 199727 : }
1073 :
1074 :
1075 : /**
1076 : * This method projects an arbitrary function onto the solution via L2
1077 : * projections and nodal interpolations on each element.
1078 : */
1079 213 : void System::project_solution (FEMFunctionBase<Number> * f,
1080 : FEMFunctionBase<Gradient> * g,
1081 : std::optional<ConstElemRange> active_local_range,
1082 : std::optional<std::vector<unsigned int>> variable_numbers) const
1083 : {
1084 420 : this->project_vector(*solution, f, g, /*is_adjoint=*/-1, active_local_range, variable_numbers);
1085 :
1086 219 : solution->localize(*current_local_solution, _dof_map->get_send_list());
1087 213 : }
1088 :
1089 :
1090 : /**
1091 : * This method projects an arbitrary function via L2 projections and
1092 : * nodal interpolations on each element.
1093 : */
1094 840 : void System::project_vector (ValueFunctionPointer fptr,
1095 : GradientFunctionPointer gptr,
1096 : const Parameters & function_parameters,
1097 : NumericVector<Number> & new_vector,
1098 : int is_adjoint,
1099 : std::optional<ConstElemRange> active_local_range,
1100 : std::optional<std::vector<unsigned int>> variable_numbers) const
1101 : {
1102 48 : WrappedFunction<Number> f(*this, fptr, &function_parameters);
1103 24 : WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
1104 1656 : this->project_vector(new_vector, &f, &g, is_adjoint, active_local_range, variable_numbers);
1105 840 : }
1106 :
1107 : /**
1108 : * This method projects an arbitrary function via L2 projections and
1109 : * nodal interpolations on each element.
1110 : */
1111 201277 : void System::project_vector (NumericVector<Number> & new_vector,
1112 : FunctionBase<Number> * f,
1113 : FunctionBase<Gradient> * g,
1114 : int is_adjoint,
1115 : std::optional<ConstElemRange> active_local_range,
1116 : std::optional<std::vector<unsigned int>> variable_numbers) const
1117 : {
1118 11344 : LOG_SCOPE ("project_vector(FunctionBase)", "System");
1119 :
1120 5672 : libmesh_assert(f);
1121 :
1122 11344 : WrappedFunctor<Number> f_fem(*f);
1123 :
1124 201277 : if (g)
1125 : {
1126 5634 : WrappedFunctor<Gradient> g_fem(*g);
1127 :
1128 394423 : this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint, active_local_range, variable_numbers);
1129 : }
1130 : else
1131 2528 : this->project_vector(new_vector, &f_fem, nullptr, is_adjoint, active_local_range, variable_numbers);
1132 201277 : }
1133 :
1134 :
1135 : /**
1136 : * This method projects an arbitrary function via L2 projections and
1137 : * nodal interpolations on each element.
1138 : */
1139 201490 : void System::project_vector (NumericVector<Number> & new_vector,
1140 : FEMFunctionBase<Number> * f,
1141 : FEMFunctionBase<Gradient> * g,
1142 : int is_adjoint,
1143 : std::optional<ConstElemRange> active_local_range,
1144 : std::optional<std::vector<unsigned int>> variable_numbers) const
1145 : {
1146 11356 : LOG_SCOPE ("project_fem_vector()", "System");
1147 :
1148 5678 : libmesh_assert (f);
1149 :
1150 201490 : if (!active_local_range)
1151 : {
1152 : active_local_range.emplace
1153 402838 : (this->get_mesh().active_local_elements_begin(),
1154 414190 : this->get_mesh().active_local_elements_end());
1155 : }
1156 :
1157 5678 : VectorSetAction<Number> setter(new_vector);
1158 :
1159 201490 : const unsigned int n_variables = this->n_vars();
1160 :
1161 11356 : std::vector<unsigned int> vars;
1162 201490 : if (variable_numbers)
1163 : {
1164 71 : vars = *variable_numbers;
1165 142 : for (auto v : vars)
1166 71 : if (v >= n_variables)
1167 0 : libmesh_error_msg("ERROR: variable number " << v <<
1168 : " out of range for system with " <<
1169 : n_variables << " variables.");
1170 : }
1171 : else
1172 : {
1173 201419 : vars.resize(n_variables);
1174 5676 : std::iota(vars.begin(), vars.end(), 0);
1175 : }
1176 :
1177 :
1178 : // Use a typedef to make the calling sequence for parallel_for() a bit more readable
1179 : typedef
1180 : GenericProjector<FEMFunctionWrapper<Number>, FEMFunctionWrapper<Gradient>,
1181 : Number, VectorSetAction<Number>> FEMProjector;
1182 :
1183 11356 : FEMFunctionWrapper<Number> fw(*f);
1184 :
1185 201490 : if (g)
1186 : {
1187 11268 : FEMFunctionWrapper<Gradient> gw(*g);
1188 :
1189 211262 : FEMProjector projector(*this, fw, &gw, setter, vars);
1190 199994 : projector.project(active_local_range.value());
1191 188726 : }
1192 : else
1193 : {
1194 1584 : FEMProjector projector(*this, fw, nullptr, setter, vars);
1195 1496 : projector.project(active_local_range.value());
1196 1408 : }
1197 :
1198 : // Also, load values into the SCALAR dofs
1199 : // Note: We assume that all SCALAR dofs are on the
1200 : // processor with highest ID
1201 207168 : if (this->processor_id() == (this->n_processors()-1))
1202 : {
1203 : // FIXME: Do we want to first check for SCALAR vars before building this? [PB]
1204 39716 : FEMContext context( *this );
1205 :
1206 2839 : const DofMap & dof_map = this->get_dof_map();
1207 68664 : for (auto var : vars)
1208 34626 : if (this->variable(var).type().family == SCALAR)
1209 : {
1210 : // FIXME: We reinit with an arbitrary element in case the user
1211 : // doesn't override FEMFunctionBase::component. Is there
1212 : // any use case we're missing? [PB]
1213 0 : context.pre_fe_reinit(*this, *(this->get_mesh().active_local_elements_begin()));
1214 :
1215 0 : std::vector<dof_id_type> SCALAR_indices;
1216 0 : dof_map.SCALAR_dof_indices (SCALAR_indices, var);
1217 : const unsigned int n_SCALAR_dofs =
1218 0 : cast_int<unsigned int>(SCALAR_indices.size());
1219 :
1220 0 : for (unsigned int i=0; i<n_SCALAR_dofs; i++)
1221 : {
1222 0 : const dof_id_type global_index = SCALAR_indices[i];
1223 : const unsigned int component_index =
1224 0 : this->variable_scalar_number(var,i);
1225 :
1226 0 : new_vector.set(global_index, f->component(context, component_index, Point(), this->time));
1227 : }
1228 : }
1229 28360 : }
1230 :
1231 201490 : new_vector.close();
1232 :
1233 : // Look for spline bases, in which case we need to backtrack
1234 : // to calculate the spline DoF values.
1235 11356 : std::vector<const Variable *> rational_vars;
1236 406459 : for (auto varnum : vars)
1237 : {
1238 204969 : const Variable & var = this->get_dof_map().variable(varnum);
1239 204969 : if (var.type().family == RATIONAL_BERNSTEIN)
1240 5412 : rational_vars.push_back(&var);
1241 : }
1242 :
1243 : // Okay, but are we really using any *spline* bases, or just
1244 : // unconstrained rational bases?
1245 201490 : bool using_spline_bases = false;
1246 201490 : if (!rational_vars.empty())
1247 : {
1248 : // Look for a spline node: a NodeElem with a rational variable
1249 : // on it.
1250 8906 : for (auto & elem : active_local_range.value())
1251 3794 : if (elem->type() == NODEELEM)
1252 300 : for (auto rational_var : rational_vars)
1253 300 : if (rational_var->active_on_subdomain(elem->subdomain_id()))
1254 : {
1255 300 : using_spline_bases = true;
1256 292 : goto checked_on_splines;
1257 : }
1258 : }
1259 :
1260 201186 : checked_on_splines:
1261 :
1262 : // Not every processor may have a NodeElem, especially while
1263 : // we're not partitioning them efficiently yet.
1264 201490 : this->comm().max(using_spline_bases);
1265 :
1266 201490 : if (using_spline_bases)
1267 300 : this->solve_for_unconstrained_dofs(new_vector, is_adjoint);
1268 :
1269 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1270 201490 : if (is_adjoint == -1)
1271 201490 : this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
1272 0 : else if (is_adjoint >= 0)
1273 0 : this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
1274 : is_adjoint);
1275 : #else
1276 : libmesh_ignore(is_adjoint);
1277 : #endif
1278 201490 : }
1279 :
1280 :
1281 : /**
1282 : * This method projects components of an arbitrary boundary function
1283 : * onto the solution via L2 projections and nodal interpolations on
1284 : * each element.
1285 : */
1286 0 : void System::boundary_project_solution (const std::set<boundary_id_type> & b,
1287 : const std::vector<unsigned int> & variables,
1288 : ValueFunctionPointer fptr,
1289 : GradientFunctionPointer gptr,
1290 : const Parameters & function_parameters,
1291 : std::optional<ConstElemRange> active_local_range)
1292 :
1293 : {
1294 0 : WrappedFunction<Number> f(*this, fptr, &function_parameters);
1295 0 : WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
1296 0 : this->boundary_project_solution(b, variables, &f, &g, active_local_range);
1297 0 : }
1298 :
1299 :
1300 : /**
1301 : * This method projects an arbitrary boundary function onto the
1302 : * solution via L2 projections and nodal interpolations on each
1303 : * element.
1304 : */
1305 71 : void System::boundary_project_solution (const std::set<boundary_id_type> & b,
1306 : const std::vector<unsigned int> & variables,
1307 : FunctionBase<Number> * f,
1308 : FunctionBase<Gradient> * g,
1309 : std::optional<ConstElemRange> active_local_range)
1310 : {
1311 140 : this->boundary_project_vector(b, variables, *solution, f, g, -1 /*is_adjoint*/, active_local_range);
1312 :
1313 73 : solution->localize(*current_local_solution);
1314 71 : }
1315 :
1316 :
1317 :
1318 :
1319 :
1320 : /**
1321 : * This method projects an arbitrary boundary function via L2
1322 : * projections and nodal interpolations on each element.
1323 : */
1324 0 : void System::boundary_project_vector (const std::set<boundary_id_type> & b,
1325 : const std::vector<unsigned int> & variables,
1326 : ValueFunctionPointer fptr,
1327 : GradientFunctionPointer gptr,
1328 : const Parameters & function_parameters,
1329 : NumericVector<Number> & new_vector,
1330 : int is_adjoint,
1331 : std::optional<ConstElemRange> active_local_range) const
1332 : {
1333 0 : WrappedFunction<Number> f(*this, fptr, &function_parameters);
1334 0 : WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
1335 0 : this->boundary_project_vector(b, variables, new_vector, &f, &g,
1336 : is_adjoint, active_local_range);
1337 0 : }
1338 :
1339 : /**
1340 : * This method projects an arbitrary function via L2 projections and
1341 : * nodal interpolations on each element.
1342 : */
1343 71 : void System::boundary_project_vector (const std::set<boundary_id_type> & b,
1344 : const std::vector<unsigned int> & variables,
1345 : NumericVector<Number> & new_vector,
1346 : FunctionBase<Number> * f,
1347 : FunctionBase<Gradient> * g,
1348 : int is_adjoint,
1349 : std::optional<ConstElemRange> active_local_range) const
1350 : {
1351 4 : LOG_SCOPE ("boundary_project_vector()", "System");
1352 :
1353 71 : if (!active_local_range)
1354 : {
1355 : active_local_range.emplace
1356 142 : (this->get_mesh().active_local_elements_begin(),
1357 146 : this->get_mesh().active_local_elements_end());
1358 : }
1359 :
1360 : Threads::parallel_for
1361 73 : (active_local_range.value(),
1362 138 : BoundaryProjectSolution(b, variables, *this, f, g,
1363 71 : this->get_equation_systems().parameters,
1364 : new_vector)
1365 : );
1366 :
1367 : // We don't do SCALAR dofs when just projecting the boundary, so
1368 : // we're done here.
1369 :
1370 71 : new_vector.close();
1371 :
1372 : #ifdef LIBMESH_ENABLE_CONSTRAINTS
1373 71 : if (is_adjoint == -1)
1374 71 : this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
1375 0 : else if (is_adjoint >= 0)
1376 0 : this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
1377 : is_adjoint);
1378 : #else
1379 : libmesh_ignore(is_adjoint);
1380 : #endif
1381 71 : }
1382 :
1383 :
1384 :
1385 : #ifdef LIBMESH_ENABLE_AMR
1386 56932 : void BuildProjectionList::unique()
1387 : {
1388 : // Sort the send list. After this duplicated
1389 : // elements will be adjacent in the vector
1390 56932 : std::sort(this->send_list.begin(),
1391 : this->send_list.end());
1392 :
1393 : // Now use std::unique to remove duplicate entries
1394 : std::vector<dof_id_type>::iterator new_end =
1395 52940 : std::unique (this->send_list.begin(),
1396 3992 : this->send_list.end());
1397 :
1398 : // Remove the end of the send_list. Use the "swap trick"
1399 : // from Effective STL
1400 58928 : std::vector<dof_id_type>
1401 1996 : (this->send_list.begin(), new_end).swap (this->send_list);
1402 56932 : }
1403 :
1404 :
1405 :
1406 60814 : void BuildProjectionList::operator()(const ConstElemRange & range)
1407 : {
1408 : // The DofMap for this system
1409 60814 : const DofMap & dof_map = system.get_dof_map();
1410 :
1411 3418 : const dof_id_type first_old_dof = dof_map.first_old_dof();
1412 3418 : const dof_id_type end_old_dof = dof_map.end_old_dof();
1413 :
1414 : // We can handle all the variables at once.
1415 : // The old global DOF indices
1416 6836 : std::vector<dof_id_type> di;
1417 :
1418 : // Iterate over the elements in the range
1419 4299105 : for (const auto & elem : range)
1420 : {
1421 : // If this element doesn't have an old_dof_object with dofs for the
1422 : // current system, then it must be newly added, so the user
1423 : // is responsible for setting the new dofs.
1424 :
1425 : // ... but we need a better way to test for that; the code
1426 : // below breaks on any FE type for which the elem stores no
1427 : // dofs.
1428 : // if (!elem->get_old_dof_object() || !elem->get_old_dof_object()->has_dofs(system.number()))
1429 : // continue;
1430 :
1431 : // Examining refinement flags instead should distinguish
1432 : // between refinement-added and user-added elements lacking
1433 : // old_dof_object
1434 4238291 : const DofObject * old_dof_object = elem->get_old_dof_object();
1435 2119643 : if (!old_dof_object &&
1436 4370073 : elem->refinement_flag() != Elem::JUST_REFINED &&
1437 12070 : elem->refinement_flag() != Elem::JUST_COARSENED)
1438 124482 : continue;
1439 :
1440 920284 : const Elem * parent = elem->parent();
1441 :
1442 4562579 : if (elem->refinement_flag() == Elem::JUST_REFINED)
1443 : {
1444 141898 : libmesh_assert(parent);
1445 :
1446 : // We used to hack_p_level here, but that wasn't thread-safe
1447 : // so now we take p refinement flags into account in
1448 : // old_dof_indices
1449 :
1450 1492916 : dof_map.old_dof_indices (parent, di);
1451 :
1452 11417778 : for (auto & node : elem->node_ref_range())
1453 : {
1454 921438 : const DofObject * old_dofs = node.get_old_dof_object();
1455 :
1456 9782956 : if (old_dofs)
1457 : {
1458 3615272 : const unsigned int sysnum = system.number();
1459 3615272 : const unsigned int nvg = old_dofs->n_var_groups(sysnum);
1460 :
1461 7320054 : for (unsigned int vg=0; vg != nvg; ++vg)
1462 : {
1463 : const unsigned int nvig =
1464 349722 : old_dofs->n_vars(sysnum, vg);
1465 7418448 : for (unsigned int vig=0; vig != nvig; ++vig)
1466 : {
1467 : const unsigned int n_comp =
1468 350494 : old_dofs->n_comp_group(sysnum, vg);
1469 7704831 : for (unsigned int c=0; c != n_comp; ++c)
1470 : {
1471 : const dof_id_type old_id =
1472 3616036 : old_dofs->dof_number(sysnum, vg, vig,
1473 3991165 : c, n_comp);
1474 :
1475 : // We should either have no old id
1476 : // (e.g. on a newly expanded subdomain)
1477 : // or an id from the old system.
1478 375035 : libmesh_assert(old_id < dof_map.n_old_dofs() ||
1479 : old_id == DofObject::invalid_id);
1480 3991165 : di.push_back(old_id);
1481 : }
1482 : }
1483 : }
1484 : }
1485 : }
1486 :
1487 1492916 : std::sort(di.begin(), di.end());
1488 : std::vector<dof_id_type>::iterator new_end =
1489 1492916 : std::unique(di.begin(), di.end());
1490 2843934 : std::vector<dof_id_type>(di.begin(), new_end).swap(di);
1491 : }
1492 2609532 : else if (elem->refinement_flag() == Elem::JUST_COARSENED)
1493 : {
1494 45750 : std::vector<dof_id_type> di_child;
1495 22875 : di.clear();
1496 1054806 : for (auto & child : elem->child_ref_range())
1497 : {
1498 844056 : dof_map.old_dof_indices (&child, di_child);
1499 844056 : di.insert(di.end(), di_child.begin(), di_child.end());
1500 : }
1501 : }
1502 : else
1503 2398782 : dof_map.old_dof_indices (elem, di);
1504 :
1505 33140883 : for (auto di_i : di)
1506 : {
1507 : // If we've just expanded a subdomain for a
1508 : // subdomain-restricted variable, then we may have an
1509 : // old_dof_object that doesn't have an old DoF for every
1510 : // local index.
1511 29038435 : if (di_i == DofObject::invalid_id)
1512 0 : continue;
1513 :
1514 3027853 : libmesh_assert_less(di_i, dof_map.n_old_dofs());
1515 29038435 : if (di_i < first_old_dof || di_i >= end_old_dof)
1516 13243428 : this->send_list.push_back(di_i);
1517 : }
1518 : } // end elem loop
1519 60814 : }
1520 :
1521 :
1522 :
1523 3882 : void BuildProjectionList::join(const BuildProjectionList & other)
1524 : {
1525 : // Joining simply requires I add the dof indices from the other object
1526 2460 : this->send_list.insert(this->send_list.end(),
1527 : other.send_list.begin(),
1528 5688 : other.send_list.end());
1529 3882 : }
1530 : #endif // LIBMESH_ENABLE_AMR
1531 :
1532 :
1533 :
1534 77 : void BoundaryProjectSolution::operator()(const ConstElemRange & range) const
1535 : {
1536 : // We need data to project
1537 4 : libmesh_assert(f.get());
1538 :
1539 : /**
1540 : * This method projects an arbitrary boundary solution to the current
1541 : * mesh. The input function \p f gives the arbitrary solution,
1542 : * while the \p new_vector (which should already be correctly sized)
1543 : * gives the solution (to be computed) on the current mesh.
1544 : */
1545 :
1546 : // The dimensionality of the current mesh
1547 77 : const unsigned int dim = system.get_mesh().mesh_dimension();
1548 :
1549 : // The DofMap for this system
1550 77 : const DofMap & dof_map = system.get_dof_map();
1551 :
1552 : // Boundary info for the current mesh
1553 : const BoundaryInfo & boundary_info =
1554 8 : system.get_mesh().get_boundary_info();
1555 :
1556 : // The element matrix and RHS for projections.
1557 : // Note that Ke is always real-valued, whereas
1558 : // Fe may be complex valued if complex number
1559 : // support is enabled
1560 85 : DenseMatrix<Real> Ke;
1561 77 : DenseVector<Number> Fe;
1562 : // The new element coefficients
1563 77 : DenseVector<Number> Ue;
1564 :
1565 :
1566 : // Loop over all the variables we've been requested to project
1567 154 : for (auto v : make_range(variables.size()))
1568 : {
1569 77 : const unsigned int var = variables[v];
1570 :
1571 77 : const Variable & variable = dof_map.variable(var);
1572 :
1573 4 : const FEType & fe_type = variable.type();
1574 :
1575 77 : if (fe_type.family == SCALAR)
1576 0 : continue;
1577 :
1578 : const unsigned int var_component =
1579 77 : system.variable_scalar_number(var, 0);
1580 :
1581 : // Get FE objects of the appropriate type
1582 81 : std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
1583 :
1584 : // Prepare variables for projection
1585 81 : std::unique_ptr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
1586 81 : std::unique_ptr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
1587 :
1588 : // The values of the shape functions at the quadrature
1589 : // points
1590 4 : const std::vector<std::vector<Real>> & phi = fe->get_phi();
1591 :
1592 : // The gradients of the shape functions at the quadrature
1593 : // points on the child element.
1594 4 : const std::vector<std::vector<RealGradient>> * dphi = nullptr;
1595 :
1596 77 : const FEContinuity cont = fe->get_continuity();
1597 :
1598 77 : if (cont == C_ONE)
1599 : {
1600 : // We'll need gradient data for a C1 projection
1601 0 : libmesh_assert(g.get());
1602 :
1603 : const std::vector<std::vector<RealGradient>> &
1604 0 : ref_dphi = fe->get_dphi();
1605 0 : dphi = &ref_dphi;
1606 : }
1607 :
1608 : // The Jacobian * quadrature weight at the quadrature points
1609 : const std::vector<Real> & JxW =
1610 9 : fe->get_JxW();
1611 :
1612 : // The XYZ locations of the quadrature points
1613 : const std::vector<Point> & xyz_values =
1614 9 : fe->get_xyz();
1615 :
1616 : // The global DOF indices
1617 8 : std::vector<dof_id_type> dof_indices;
1618 : // Side/edge DOF indices
1619 8 : std::vector<unsigned int> side_dofs;
1620 :
1621 : // Container to catch IDs passed back from BoundaryInfo.
1622 8 : std::vector<boundary_id_type> bc_ids;
1623 :
1624 : // Iterate over all the elements in the range
1625 401 : for (const auto & elem : range)
1626 : {
1627 : // Per-subdomain variables don't need to be projected on
1628 : // elements where they're not active
1629 324 : if (!variable.active_on_subdomain(elem->subdomain_id()))
1630 0 : continue;
1631 :
1632 324 : const unsigned short n_nodes = elem->n_nodes();
1633 324 : const unsigned short n_edges = elem->n_edges();
1634 324 : const unsigned short n_sides = elem->n_sides();
1635 :
1636 : // Find out which nodes, edges and sides are on a requested
1637 : // boundary:
1638 351 : std::vector<bool> is_boundary_node(n_nodes, false),
1639 351 : is_boundary_edge(n_edges, false),
1640 351 : is_boundary_side(n_sides, false);
1641 :
1642 : // We also maintain a separate list of nodeset-based boundary nodes
1643 351 : std::vector<bool> is_boundary_nodeset(n_nodes, false);
1644 :
1645 2268 : for (unsigned char s=0; s != n_sides; ++s)
1646 : {
1647 : // First see if this side has been requested
1648 1944 : boundary_info.boundary_ids (elem, s, bc_ids);
1649 162 : bool do_this_side = false;
1650 2484 : for (const auto & bc_id : bc_ids)
1651 648 : if (b.count(bc_id))
1652 : {
1653 9 : do_this_side = true;
1654 9 : break;
1655 : }
1656 1944 : if (!do_this_side)
1657 1683 : continue;
1658 :
1659 18 : is_boundary_side[s] = true;
1660 :
1661 : // Then see what nodes and what edges are on it
1662 972 : for (unsigned int n=0; n != n_nodes; ++n)
1663 864 : if (elem->is_node_on_side(n,s))
1664 72 : is_boundary_node[n] = true;
1665 1404 : for (unsigned int e=0; e != n_edges; ++e)
1666 1296 : if (elem->is_edge_on_side(e,s))
1667 72 : is_boundary_edge[e] = true;
1668 : }
1669 :
1670 : // We can also project on nodes, so we should also independently
1671 : // check whether the nodes have been requested
1672 2916 : for (unsigned int n=0; n != n_nodes; ++n)
1673 : {
1674 2808 : boundary_info.boundary_ids (elem->node_ptr(n), bc_ids);
1675 :
1676 5196 : for (const auto & bc_id : bc_ids)
1677 2604 : if (b.count(bc_id))
1678 : {
1679 74 : is_boundary_node[n] = true;
1680 74 : is_boundary_nodeset[n] = true;
1681 : }
1682 : }
1683 :
1684 : // We can also project on edges, so we should also independently
1685 : // check whether the edges have been requested
1686 4212 : for (unsigned short e=0; e != n_edges; ++e)
1687 : {
1688 3888 : boundary_info.edge_boundary_ids (elem, e, bc_ids);
1689 :
1690 3924 : for (const auto & bc_id : bc_ids)
1691 36 : if (b.count(bc_id))
1692 6 : is_boundary_edge[e] = true;
1693 : }
1694 :
1695 : // Update the DOF indices for this element based on
1696 : // the current mesh
1697 324 : dof_map.dof_indices (elem, dof_indices, var);
1698 :
1699 : // The number of DOFs on the element
1700 : const unsigned int n_dofs =
1701 54 : cast_int<unsigned int>(dof_indices.size());
1702 :
1703 : // Fixed vs. free DoFs on edge/face projections
1704 351 : std::vector<char> dof_is_fixed(n_dofs, false); // bools
1705 351 : std::vector<int> free_dof(n_dofs, 0);
1706 :
1707 : // Zero the interpolated values
1708 297 : Ue.resize (n_dofs); Ue.zero();
1709 :
1710 : // In general, we need a series of
1711 : // projections to ensure a unique and continuous
1712 : // solution. We start by interpolating boundary nodes, then
1713 : // hold those fixed and project boundary edges, then hold
1714 : // those fixed and project boundary faces,
1715 :
1716 : // Interpolate node values first
1717 27 : unsigned int current_dof = 0;
1718 2916 : for (unsigned short n = 0; n != n_nodes; ++n)
1719 : {
1720 : // FIXME: this should go through the DofMap,
1721 : // not duplicate dof_indices code badly!
1722 :
1723 : // This call takes into account elem->p_level() internally.
1724 : const unsigned int nc =
1725 2592 : FEInterface::n_dofs_at_node (fe_type, elem, n);
1726 :
1727 2771 : if ((!elem->is_vertex(n) || !is_boundary_node[n]) &&
1728 574 : !is_boundary_nodeset[n])
1729 : {
1730 2148 : current_dof += nc;
1731 2148 : continue;
1732 : }
1733 444 : if (cont == DISCONTINUOUS)
1734 : {
1735 0 : libmesh_assert_equal_to (nc, 0);
1736 : }
1737 : // Assume that C_ZERO elements have a single nodal
1738 : // value shape function
1739 444 : else if (cont == C_ZERO)
1740 : {
1741 37 : libmesh_assert_equal_to (nc, 1);
1742 851 : Ue(current_dof) = f->component(var_component,
1743 : elem->point(n),
1744 444 : system.time);
1745 444 : dof_is_fixed[current_dof] = true;
1746 444 : current_dof++;
1747 : }
1748 : // The hermite element vertex shape functions are weird
1749 0 : else if (fe_type.family == HERMITE)
1750 : {
1751 0 : Ue(current_dof) = f->component(var_component,
1752 : elem->point(n),
1753 0 : system.time);
1754 0 : dof_is_fixed[current_dof] = true;
1755 0 : current_dof++;
1756 0 : Gradient grad = g->component(var_component,
1757 : elem->point(n),
1758 0 : system.time);
1759 : // x derivative
1760 0 : Ue(current_dof) = grad(0);
1761 0 : dof_is_fixed[current_dof] = true;
1762 0 : current_dof++;
1763 : #if LIBMESH_DIM > 1
1764 0 : if (dim > 1)
1765 : {
1766 : // We'll finite difference mixed derivatives
1767 0 : Point nxminus = elem->point(n),
1768 0 : nxplus = elem->point(n);
1769 0 : nxminus(0) -= TOLERANCE;
1770 0 : nxplus(0) += TOLERANCE;
1771 0 : Gradient gxminus = g->component(var_component,
1772 : nxminus,
1773 0 : system.time);
1774 0 : Gradient gxplus = g->component(var_component,
1775 : nxplus,
1776 0 : system.time);
1777 : // y derivative
1778 0 : Ue(current_dof) = grad(1);
1779 0 : dof_is_fixed[current_dof] = true;
1780 0 : current_dof++;
1781 : // xy derivative
1782 0 : Ue(current_dof) = (gxplus(1) - gxminus(1))
1783 0 : / 2. / TOLERANCE;
1784 0 : dof_is_fixed[current_dof] = true;
1785 0 : current_dof++;
1786 :
1787 : #if LIBMESH_DIM > 2
1788 0 : if (dim > 2)
1789 : {
1790 : // z derivative
1791 0 : Ue(current_dof) = grad(2);
1792 0 : dof_is_fixed[current_dof] = true;
1793 0 : current_dof++;
1794 : // xz derivative
1795 0 : Ue(current_dof) = (gxplus(2) - gxminus(2))
1796 0 : / 2. / TOLERANCE;
1797 0 : dof_is_fixed[current_dof] = true;
1798 0 : current_dof++;
1799 : // We need new points for yz
1800 0 : Point nyminus = elem->point(n),
1801 0 : nyplus = elem->point(n);
1802 0 : nyminus(1) -= TOLERANCE;
1803 0 : nyplus(1) += TOLERANCE;
1804 0 : Gradient gyminus = g->component(var_component,
1805 : nyminus,
1806 0 : system.time);
1807 0 : Gradient gyplus = g->component(var_component,
1808 : nyplus,
1809 0 : system.time);
1810 : // xz derivative
1811 0 : Ue(current_dof) = (gyplus(2) - gyminus(2))
1812 0 : / 2. / TOLERANCE;
1813 0 : dof_is_fixed[current_dof] = true;
1814 0 : current_dof++;
1815 : // Getting a 2nd order xyz is more tedious
1816 0 : Point nxmym = elem->point(n),
1817 0 : nxmyp = elem->point(n),
1818 0 : nxpym = elem->point(n),
1819 0 : nxpyp = elem->point(n);
1820 0 : nxmym(0) -= TOLERANCE;
1821 0 : nxmym(1) -= TOLERANCE;
1822 0 : nxmyp(0) -= TOLERANCE;
1823 0 : nxmyp(1) += TOLERANCE;
1824 0 : nxpym(0) += TOLERANCE;
1825 0 : nxpym(1) -= TOLERANCE;
1826 0 : nxpyp(0) += TOLERANCE;
1827 0 : nxpyp(1) += TOLERANCE;
1828 0 : Gradient gxmym = g->component(var_component,
1829 : nxmym,
1830 0 : system.time);
1831 0 : Gradient gxmyp = g->component(var_component,
1832 : nxmyp,
1833 0 : system.time);
1834 0 : Gradient gxpym = g->component(var_component,
1835 : nxpym,
1836 0 : system.time);
1837 0 : Gradient gxpyp = g->component(var_component,
1838 : nxpyp,
1839 0 : system.time);
1840 0 : Number gxzplus = (gxpyp(2) - gxmyp(2))
1841 0 : / 2. / TOLERANCE;
1842 0 : Number gxzminus = (gxpym(2) - gxmym(2))
1843 0 : / 2. / TOLERANCE;
1844 : // xyz derivative
1845 0 : Ue(current_dof) = (gxzplus - gxzminus)
1846 0 : / 2. / TOLERANCE;
1847 0 : dof_is_fixed[current_dof] = true;
1848 0 : current_dof++;
1849 : }
1850 : #endif // LIBMESH_DIM > 2
1851 : }
1852 : #endif // LIBMESH_DIM > 1
1853 : }
1854 : // Assume that other C_ONE elements have a single nodal
1855 : // value shape function and nodal gradient component
1856 : // shape functions
1857 0 : else if (cont == C_ONE)
1858 : {
1859 0 : libmesh_assert_equal_to (nc, 1 + dim);
1860 0 : Ue(current_dof) = f->component(var_component,
1861 : elem->point(n),
1862 0 : system.time);
1863 0 : dof_is_fixed[current_dof] = true;
1864 0 : current_dof++;
1865 0 : Gradient grad = g->component(var_component,
1866 : elem->point(n),
1867 0 : system.time);
1868 0 : for (unsigned int i=0; i!= dim; ++i)
1869 : {
1870 0 : Ue(current_dof) = grad(i);
1871 0 : dof_is_fixed[current_dof] = true;
1872 0 : current_dof++;
1873 : }
1874 : }
1875 : else
1876 0 : libmesh_error_msg("Unknown continuity " << cont);
1877 : }
1878 :
1879 : // In 3D, project any edge values next
1880 324 : if (dim > 2 && cont != DISCONTINUOUS)
1881 4212 : for (unsigned short e = 0; e != n_edges; ++e)
1882 : {
1883 4212 : if (!is_boundary_edge[e])
1884 3852 : continue;
1885 :
1886 468 : FEInterface::dofs_on_edge(elem, dim, fe_type, e,
1887 : side_dofs);
1888 :
1889 : const unsigned int n_side_dofs =
1890 78 : cast_int<unsigned int>(side_dofs.size());
1891 :
1892 : // Some edge dofs are on nodes and already
1893 : // fixed, others are free to calculate
1894 39 : unsigned int free_dofs = 0;
1895 1404 : for (auto i : make_range(n_side_dofs))
1896 1092 : if (!dof_is_fixed[side_dofs[i]])
1897 78 : free_dof[free_dofs++] = i;
1898 :
1899 : // There may be nothing to project
1900 468 : if (!free_dofs)
1901 396 : continue;
1902 :
1903 33 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1904 33 : Fe.resize (free_dofs); Fe.zero();
1905 : // The new edge coefficients
1906 36 : DenseVector<Number> Uedge(free_dofs);
1907 :
1908 : // Initialize FE data on the edge
1909 39 : fe->attach_quadrature_rule (qedgerule.get());
1910 36 : fe->edge_reinit (elem, e);
1911 3 : const unsigned int n_qp = qedgerule->n_points();
1912 :
1913 : // Loop over the quadrature points
1914 108 : for (unsigned int qp=0; qp<n_qp; qp++)
1915 : {
1916 : // solution at the quadrature point
1917 78 : Number fineval = f->component(var_component,
1918 72 : xyz_values[qp],
1919 72 : system.time);
1920 : // solution grad at the quadrature point
1921 6 : Gradient finegrad;
1922 72 : if (cont == C_ONE)
1923 0 : finegrad = g->component(var_component,
1924 0 : xyz_values[qp],
1925 0 : system.time);
1926 :
1927 : // Form edge projection matrix
1928 204 : for (unsigned int sidei=0, freei=0;
1929 216 : sidei != n_side_dofs; ++sidei)
1930 : {
1931 144 : unsigned int i = side_dofs[sidei];
1932 : // fixed DoFs aren't test functions
1933 156 : if (dof_is_fixed[i])
1934 0 : continue;
1935 300 : for (unsigned int sidej=0, freej=0;
1936 432 : sidej != n_side_dofs; ++sidej)
1937 : {
1938 288 : unsigned int j = side_dofs[sidej];
1939 312 : if (dof_is_fixed[j])
1940 0 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
1941 0 : JxW[qp] * Ue(j);
1942 : else
1943 384 : Ke(freei,freej) += phi[i][qp] *
1944 336 : phi[j][qp] * JxW[qp];
1945 288 : if (cont == C_ONE)
1946 : {
1947 0 : if (dof_is_fixed[j])
1948 0 : Fe(freei) -= ((*dphi)[i][qp] *
1949 0 : (*dphi)[j][qp]) *
1950 0 : JxW[qp] * Ue(j);
1951 : else
1952 0 : Ke(freei,freej) += ((*dphi)[i][qp] *
1953 0 : (*dphi)[j][qp])
1954 0 : * JxW[qp];
1955 : }
1956 312 : if (!dof_is_fixed[j])
1957 288 : freej++;
1958 : }
1959 168 : Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1960 144 : if (cont == C_ONE)
1961 0 : Fe(freei) += (finegrad * (*dphi)[i][qp]) *
1962 0 : JxW[qp];
1963 144 : freei++;
1964 : }
1965 : }
1966 :
1967 36 : Ke.cholesky_solve(Fe, Uedge);
1968 :
1969 : // Transfer new edge solutions to element
1970 108 : for (unsigned int i=0; i != free_dofs; ++i)
1971 : {
1972 84 : Number & ui = Ue(side_dofs[free_dof[i]]);
1973 6 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1974 : std::abs(ui - Uedge(i)) < TOLERANCE);
1975 72 : ui = Uedge(i);
1976 78 : dof_is_fixed[side_dofs[free_dof[i]]] = true;
1977 : }
1978 : }
1979 :
1980 : // Project any side values (edges in 2D, faces in 3D)
1981 324 : if (dim > 1 && cont != DISCONTINUOUS)
1982 2268 : for (unsigned short s = 0; s != n_sides; ++s)
1983 : {
1984 2106 : if (!is_boundary_side[s])
1985 1944 : continue;
1986 :
1987 108 : FEInterface::dofs_on_side(elem, dim, fe_type, s,
1988 : side_dofs);
1989 :
1990 : // Some side dofs are on nodes/edges and already
1991 : // fixed, others are free to calculate
1992 9 : unsigned int free_dofs = 0;
1993 540 : for (auto i : index_range(side_dofs))
1994 468 : if (!dof_is_fixed[side_dofs[i]])
1995 0 : free_dof[free_dofs++] = i;
1996 :
1997 : // There may be nothing to project
1998 108 : if (!free_dofs)
1999 99 : continue;
2000 :
2001 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
2002 0 : Fe.resize (free_dofs); Fe.zero();
2003 : // The new side coefficients
2004 0 : DenseVector<Number> Uside(free_dofs);
2005 :
2006 : // Initialize FE data on the side
2007 0 : fe->attach_quadrature_rule (qsiderule.get());
2008 0 : fe->reinit (elem, s);
2009 0 : const unsigned int n_qp = qsiderule->n_points();
2010 :
2011 : const unsigned int n_side_dofs =
2012 0 : cast_int<unsigned int>(side_dofs.size());
2013 :
2014 : // Loop over the quadrature points
2015 0 : for (unsigned int qp=0; qp<n_qp; qp++)
2016 : {
2017 : // solution at the quadrature point
2018 0 : Number fineval = f->component(var_component,
2019 0 : xyz_values[qp],
2020 0 : system.time);
2021 : // solution grad at the quadrature point
2022 0 : Gradient finegrad;
2023 0 : if (cont == C_ONE)
2024 0 : finegrad = g->component(var_component,
2025 0 : xyz_values[qp],
2026 0 : system.time);
2027 :
2028 : // Form side projection matrix
2029 0 : for (unsigned int sidei=0, freei=0;
2030 0 : sidei != n_side_dofs; ++sidei)
2031 : {
2032 0 : unsigned int i = side_dofs[sidei];
2033 : // fixed DoFs aren't test functions
2034 0 : if (dof_is_fixed[i])
2035 0 : continue;
2036 0 : for (unsigned int sidej=0, freej=0;
2037 0 : sidej != n_side_dofs; ++sidej)
2038 : {
2039 0 : unsigned int j = side_dofs[sidej];
2040 0 : if (dof_is_fixed[j])
2041 0 : Fe(freei) -= phi[i][qp] * phi[j][qp] *
2042 0 : JxW[qp] * Ue(j);
2043 : else
2044 0 : Ke(freei,freej) += phi[i][qp] *
2045 0 : phi[j][qp] * JxW[qp];
2046 0 : if (cont == C_ONE)
2047 : {
2048 0 : if (dof_is_fixed[j])
2049 0 : Fe(freei) -= ((*dphi)[i][qp] *
2050 0 : (*dphi)[j][qp]) *
2051 0 : JxW[qp] * Ue(j);
2052 : else
2053 0 : Ke(freei,freej) += ((*dphi)[i][qp] *
2054 0 : (*dphi)[j][qp])
2055 0 : * JxW[qp];
2056 : }
2057 0 : if (!dof_is_fixed[j])
2058 0 : freej++;
2059 : }
2060 0 : Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
2061 0 : if (cont == C_ONE)
2062 0 : Fe(freei) += (finegrad * (*dphi)[i][qp]) *
2063 0 : JxW[qp];
2064 0 : freei++;
2065 : }
2066 : }
2067 :
2068 0 : Ke.cholesky_solve(Fe, Uside);
2069 :
2070 : // Transfer new side solutions to element
2071 0 : for (unsigned int i=0; i != free_dofs; ++i)
2072 : {
2073 0 : Number & ui = Ue(side_dofs[free_dof[i]]);
2074 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
2075 : std::abs(ui - Uside(i)) < TOLERANCE);
2076 0 : ui = Uside(i);
2077 0 : dof_is_fixed[side_dofs[free_dof[i]]] = true;
2078 : }
2079 : }
2080 :
2081 : const dof_id_type
2082 324 : first = new_vector.first_local_index(),
2083 324 : last = new_vector.last_local_index();
2084 :
2085 : // Lock the new_vector since it is shared among threads.
2086 : {
2087 54 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2088 :
2089 2916 : for (unsigned int i = 0; i < n_dofs; i++)
2090 2592 : if (dof_is_fixed[i] &&
2091 2628 : (dof_indices[i] >= first) &&
2092 36 : (dof_indices[i] < last))
2093 : {
2094 434 : new_vector.set(dof_indices[i], Ue(i));
2095 : }
2096 : }
2097 : } // end elem loop
2098 69 : } // end variables loop
2099 77 : }
2100 :
2101 :
2102 300 : void System::solve_for_unconstrained_dofs(NumericVector<Number> & vec,
2103 : int is_adjoint) const
2104 : {
2105 8 : const DofMap & dof_map = this->get_dof_map();
2106 :
2107 : std::unique_ptr<SparseMatrix<Number>> mat =
2108 308 : SparseMatrix<Number>::build(this->comm());
2109 :
2110 292 : std::unique_ptr<SparsityPattern::Build> sp;
2111 :
2112 300 : if (dof_map.computed_sparsity_already())
2113 0 : dof_map.update_sparsity_pattern(*mat);
2114 : else
2115 : {
2116 300 : mat->attach_dof_map(dof_map);
2117 584 : sp = dof_map.build_sparsity(this->get_mesh());
2118 300 : mat->attach_sparsity_pattern(*sp);
2119 : }
2120 :
2121 300 : mat->init();
2122 :
2123 8 : libmesh_assert_equal_to(vec.size(), dof_map.n_dofs());
2124 8 : libmesh_assert_equal_to(vec.local_size(), dof_map.n_local_dofs());
2125 :
2126 : std::unique_ptr<NumericVector<Number>> rhs =
2127 308 : NumericVector<Number>::build(this->comm());
2128 :
2129 308 : rhs->init(dof_map.n_dofs(), dof_map.n_local_dofs(), false,
2130 16 : PARALLEL);
2131 :
2132 : // Here we start with the unconstrained (and indeterminate) linear
2133 : // system, K*u = f, where K is the identity matrix for constrained
2134 : // DoFs and 0 elsewhere, and f is the current solution values for
2135 : // constrained DoFs and 0 elsewhere.
2136 : // We then apply the usual heterogeneous constraint matrix C and
2137 : // offset h, where u = C*x + h,
2138 : // to get C^T*K*C*x = C^T*f - C^T*K*h
2139 : // - a constrained and no-longer-singular system that finds the
2140 : // closest approximation for the unconstrained degrees of freedom.
2141 : //
2142 : // Here, though "closest" is in an algebraic sense; we're
2143 : // effectively using a pseudoinverse that optimizes in a
2144 : // discretization-dependent norm. That only seems to give ~0.1%
2145 : // excess error even in coarse unit test cases, but at some point it
2146 : // might be reasonable to weight K and f properly.
2147 :
2148 2121 : for (dof_id_type d : IntRange<dof_id_type>(dof_map.first_dof(),
2149 24936 : dof_map.end_dof()))
2150 : {
2151 6799 : if (dof_map.is_constrained_dof(d))
2152 : {
2153 19619 : DenseMatrix<Number> K(1,1);
2154 16891 : DenseVector<Number> F(1);
2155 18255 : std::vector<dof_id_type> dof_indices(1, d);
2156 15527 : K(0,0) = 1;
2157 16891 : F(0) = (*this->solution)(d);
2158 : dof_map.heterogenously_constrain_element_matrix_and_vector
2159 1364 : (K, F, dof_indices, false, is_adjoint);
2160 16891 : mat->add_matrix(K, dof_indices);
2161 15527 : rhs->add_vector(F, dof_indices);
2162 14163 : }
2163 : }
2164 :
2165 : std::unique_ptr<LinearSolver<Number>> linear_solver =
2166 300 : LinearSolver<Number>::build(this->comm());
2167 :
2168 584 : linear_solver->solve(*mat, vec, *rhs,
2169 300 : double(this->get_equation_systems().parameters.get<Real>("linear solver tolerance")),
2170 40 : this->get_equation_systems().parameters.get<unsigned int>("linear solver maximum iterations"));
2171 300 : }
2172 :
2173 :
2174 : } // namespace libMesh
|
