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