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