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