Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "MultiAppProjectionTransfer.h"
11 :
12 : // MOOSE includes
13 : #include "AddVariableAction.h"
14 : #include "FEProblem.h"
15 : #include "MooseMesh.h"
16 : #include "MooseVariableFE.h"
17 : #include "SystemBase.h"
18 : #include "MooseAppCoordTransform.h"
19 :
20 : #include "libmesh/dof_map.h"
21 : #include "libmesh/linear_implicit_system.h"
22 : #include "libmesh/mesh_function.h"
23 : #include "libmesh/mesh_tools.h"
24 : #include "libmesh/numeric_vector.h"
25 : #include "libmesh/parallel_algebra.h"
26 : #include "libmesh/quadrature_gauss.h"
27 : #include "libmesh/sparse_matrix.h"
28 : #include "libmesh/string_to_enum.h"
29 :
30 : // TIMPI includes
31 : #include "timpi/parallel_sync.h"
32 :
33 : using namespace libMesh;
34 :
35 : void
36 628 : assemble_l2(EquationSystems & es, const std::string & system_name)
37 : {
38 : MultiAppProjectionTransfer * transfer =
39 628 : es.parameters.get<MultiAppProjectionTransfer *>("transfer");
40 628 : transfer->assembleL2(es, system_name);
41 628 : }
42 :
43 : registerMooseObject("MooseApp", MultiAppProjectionTransfer);
44 :
45 : InputParameters
46 14865 : MultiAppProjectionTransfer::validParams()
47 : {
48 14865 : InputParameters params = MultiAppConservativeTransfer::validParams();
49 14865 : params.addClassDescription(
50 : "Perform a projection between a master and sub-application mesh of a field variable.");
51 :
52 14865 : MooseEnum proj_type("l2", "l2");
53 14865 : params.addParam<MooseEnum>("proj_type", proj_type, "The type of the projection.");
54 :
55 44595 : params.addParam<bool>("fixed_meshes",
56 29730 : false,
57 : "Set to true when the meshes are not changing (ie, "
58 : "no movement or adaptivity). This will cache some "
59 : "information to speed up the transfer.");
60 :
61 : // Need one layer of ghosting
62 14865 : params.addRelationshipManager("ElementSideNeighborLayers",
63 : Moose::RelationshipManagerType::GEOMETRIC |
64 : Moose::RelationshipManagerType::ALGEBRAIC);
65 14865 : MultiAppTransfer::addBBoxFactorParam(params);
66 29730 : return params;
67 14865 : }
68 :
69 300 : MultiAppProjectionTransfer::MultiAppProjectionTransfer(const InputParameters & parameters)
70 : : MultiAppConservativeTransfer(parameters),
71 300 : _proj_type(getParam<MooseEnum>("proj_type")),
72 300 : _compute_matrix(true),
73 300 : _fixed_meshes(getParam<bool>("fixed_meshes")),
74 900 : _qps_cached(false)
75 : {
76 300 : if (_to_var_names.size() != 1)
77 0 : paramError("variable", " Support single to-variable only ");
78 :
79 300 : if (_from_var_names.size() != 1)
80 0 : paramError("source_variable", " Support single from-variable only ");
81 300 : }
82 :
83 : void
84 300 : MultiAppProjectionTransfer::initialSetup()
85 : {
86 300 : MultiAppConservativeTransfer::initialSetup();
87 :
88 300 : _proj_sys.resize(_to_problems.size(), NULL);
89 :
90 652 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
91 : {
92 352 : FEProblemBase & to_problem = *_to_problems[i_to];
93 352 : EquationSystems & to_es = to_problem.es();
94 :
95 : // Add the projection system.
96 : FEType fe_type = to_problem
97 352 : .getVariable(0,
98 : _to_var_name,
99 : Moose::VarKindType::VAR_ANY,
100 : Moose::VarFieldType::VAR_FIELD_STANDARD)
101 352 : .feType();
102 :
103 352 : LinearImplicitSystem & proj_sys = to_es.add_system<LinearImplicitSystem>("proj-sys-" + name());
104 :
105 352 : _proj_var_num = proj_sys.add_variable("var", fe_type);
106 352 : proj_sys.attach_assemble_function(assemble_l2);
107 352 : _proj_sys[i_to] = &proj_sys;
108 :
109 : // Prevent the projection system from being written to checkpoint
110 : // files. In the event of a recover or restart, we'll read the checkpoint
111 : // before this initialSetup method is called. As a result, we'll find
112 : // systems in the checkpoint (the projection systems) that we don't know
113 : // what to do with, and there will be a crash. We could fix this by making
114 : // the systems in the constructor, except we don't know how many sub apps
115 : // there are at the time of construction. So instead, we'll just nuke the
116 : // projection system and rebuild it from scratch every recover/restart.
117 352 : proj_sys.hide_output() = true;
118 :
119 : // Reinitialize EquationSystems since we added a system.
120 352 : to_es.reinit();
121 : }
122 300 : }
123 :
124 : void
125 628 : MultiAppProjectionTransfer::assembleL2(EquationSystems & es, const std::string & system_name)
126 : {
127 : // Get the system and mesh from the input arguments.
128 628 : LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);
129 628 : MeshBase & to_mesh = es.get_mesh();
130 :
131 : // Get the meshfunction evaluations and the map that was stashed in the es.
132 628 : std::vector<Real> & final_evals = *es.parameters.get<std::vector<Real> *>("final_evals");
133 : std::map<dof_id_type, unsigned int> & element_map =
134 628 : *es.parameters.get<std::map<dof_id_type, unsigned int> *>("element_map");
135 :
136 : // Setup system vectors and matrices.
137 628 : FEType fe_type = system.variable_type(0);
138 628 : std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
139 628 : QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
140 628 : fe->attach_quadrature_rule(&qrule);
141 628 : const DofMap & dof_map = system.get_dof_map();
142 628 : DenseMatrix<Number> Ke;
143 628 : DenseVector<Number> Fe;
144 628 : std::vector<dof_id_type> dof_indices;
145 628 : const std::vector<Real> & JxW = fe->get_JxW();
146 628 : const std::vector<std::vector<Real>> & phi = fe->get_phi();
147 628 : auto & system_matrix = system.get_system_matrix();
148 :
149 41472 : for (const auto & elem : to_mesh.active_local_element_ptr_range())
150 : {
151 20422 : fe->reinit(elem);
152 :
153 20422 : dof_map.dof_indices(elem, dof_indices);
154 20422 : Ke.resize(dof_indices.size(), dof_indices.size());
155 20422 : Fe.resize(dof_indices.size());
156 :
157 85390 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
158 : {
159 64968 : Real meshfun_eval = 0.;
160 64968 : if (element_map.find(elem->id()) != element_map.end())
161 : {
162 : // We have evaluations for this element.
163 52068 : meshfun_eval = final_evals[element_map[elem->id()] + qp];
164 : }
165 :
166 : // Now compute the element matrix and RHS contributions.
167 309656 : for (unsigned int i = 0; i < phi.size(); i++)
168 : {
169 : // RHS
170 244688 : Fe(i) += JxW[qp] * (meshfun_eval * phi[i][qp]);
171 :
172 244688 : if (_compute_matrix)
173 1229248 : for (unsigned int j = 0; j < phi.size(); j++)
174 : {
175 : // The matrix contribution
176 984560 : Ke(i, j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
177 : }
178 : }
179 64968 : dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
180 :
181 64968 : if (_compute_matrix)
182 64968 : system_matrix.add_matrix(Ke, dof_indices);
183 64968 : system.rhs->add_vector(Fe, dof_indices);
184 : }
185 628 : }
186 628 : }
187 :
188 : void
189 499 : MultiAppProjectionTransfer::execute()
190 : {
191 499 : TIME_SECTION(
192 : "MultiAppProjectionTransfer::execute()", 5, "Transferring variables through projection");
193 :
194 : ////////////////////
195 : // We are going to project the solutions by solving some linear systems. In
196 : // order to assemble the systems, we need to evaluate the "from" domain
197 : // solutions at quadrature points in the "to" domain. Some parallel
198 : // communication is necessary because each processor doesn't necessarily have
199 : // all the "from" information it needs to set its "to" values. We don't want
200 : // to use a bunch of big all-to-all broadcasts, so we'll use bounding boxes to
201 : // figure out which processors have the information we need and only
202 : // communicate with those processors.
203 : //
204 : // Each processor will
205 : // 1. Check its local quadrature points in the "to" domains to see which
206 : // "from" domains they might be in.
207 : // 2. Send quadrature points to the processors with "from" domains that might
208 : // contain those points.
209 : // 3. Recieve quadrature points from other processors, evaluate its mesh
210 : // functions at those points, and send the values back to the proper
211 : // processor
212 : // 4. Recieve mesh function evaluations from all relevant processors and
213 : // decide which one to use at every quadrature point (the lowest global app
214 : // index always wins)
215 : // 5. And use the mesh function evaluations to assemble and solve an L2
216 : // projection system on its local elements.
217 : ////////////////////
218 :
219 : ////////////////////
220 : // For every combination of global "from" problem and local "to" problem, find
221 : // which "from" bounding boxes overlap with which "to" elements. Keep track
222 : // of which processors own bounding boxes that overlap with which elements.
223 : // Build vectors of quadrature points to send to other processors for mesh
224 : // function evaluations.
225 : ////////////////////
226 :
227 : // Get the bounding boxes for the "from" domains.
228 499 : std::vector<BoundingBox> bboxes = getFromBoundingBoxes();
229 :
230 : // Figure out how many "from" domains each processor owns.
231 499 : std::vector<unsigned int> froms_per_proc = getFromsPerProc();
232 :
233 499 : std::map<processor_id_type, std::vector<Point>> outgoing_qps;
234 : std::map<processor_id_type, std::map<std::pair<unsigned int, unsigned int>, unsigned int>>
235 499 : element_index_map;
236 : // element_index_map[i_to, element_id] = index
237 : // outgoing_qps[index] is the first quadrature point in element
238 :
239 499 : if (!_qps_cached)
240 : {
241 1031 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
242 : {
243 : // Indexing into the coordinate transforms vector
244 : const auto to_global_num =
245 580 : _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
246 580 : MeshBase & to_mesh = _to_meshes[i_to]->getMesh();
247 :
248 580 : LinearImplicitSystem & system = *_proj_sys[i_to];
249 :
250 580 : FEType fe_type = system.variable_type(0);
251 580 : std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
252 580 : QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
253 580 : fe->attach_quadrature_rule(&qrule);
254 580 : const std::vector<Point> & xyz = fe->get_xyz();
255 :
256 580 : unsigned int from0 = 0;
257 1382 : for (processor_id_type i_proc = 0; i_proc < n_processors();
258 802 : from0 += froms_per_proc[i_proc], i_proc++)
259 : {
260 802 : for (const auto & elem :
261 51268 : as_range(to_mesh.local_elements_begin(), to_mesh.local_elements_end()))
262 : {
263 24832 : fe->reinit(elem);
264 :
265 24832 : bool qp_hit = false;
266 51680 : for (unsigned int i_from = 0; i_from < froms_per_proc[i_proc] && !qp_hit; i_from++)
267 : {
268 77993 : for (unsigned int qp = 0; qp < qrule.n_points() && !qp_hit; qp++)
269 : {
270 51145 : Point qpt = xyz[qp];
271 51145 : if (bboxes[from0 + i_from].contains_point((*_to_transforms[to_global_num])(qpt)))
272 15125 : qp_hit = true;
273 : }
274 : }
275 :
276 24832 : if (qp_hit)
277 : {
278 : // The selected processor's bounding box contains at least one
279 : // quadrature point from this element. Send all qps from this element
280 : // and remember where they are in the array using the map.
281 15125 : std::pair<unsigned int, unsigned int> key(i_to, elem->id());
282 15125 : element_index_map[i_proc][key] = outgoing_qps[i_proc].size();
283 65091 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
284 : {
285 49966 : Point qpt = xyz[qp];
286 49966 : outgoing_qps[i_proc].push_back((*_to_transforms[to_global_num])(qpt));
287 : }
288 : }
289 802 : }
290 : }
291 580 : }
292 :
293 451 : if (_fixed_meshes)
294 52 : _cached_index_map = element_index_map;
295 : }
296 : else
297 : {
298 48 : element_index_map = _cached_index_map;
299 : }
300 :
301 : ////////////////////
302 : // Request quadrature point evaluations from other processors and handle
303 : // requests sent to this processor.
304 : ////////////////////
305 :
306 : // Get the local bounding boxes.
307 499 : std::vector<BoundingBox> local_bboxes(froms_per_proc[processor_id()]);
308 : {
309 : // Find the index to the first of this processor's local bounding boxes.
310 499 : unsigned int local_start = 0;
311 622 : for (processor_id_type i_proc = 0; i_proc < n_processors() && i_proc != processor_id();
312 : i_proc++)
313 123 : local_start += froms_per_proc[i_proc];
314 :
315 : // Extract the local bounding boxes.
316 1026 : for (unsigned int i_from = 0; i_from < froms_per_proc[processor_id()]; i_from++)
317 527 : local_bboxes[i_from] = bboxes[local_start + i_from];
318 : }
319 :
320 : // Setup the local mesh functions.
321 499 : std::vector<MeshFunction> local_meshfuns;
322 1026 : for (unsigned int i_from = 0; i_from < _from_problems.size(); i_from++)
323 : {
324 527 : FEProblemBase & from_problem = *_from_problems[i_from];
325 527 : MooseVariableFEBase & from_var = from_problem.getVariable(
326 : 0, _from_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
327 527 : System & from_sys = from_var.sys().system();
328 527 : unsigned int from_var_num = from_sys.variable_number(from_var.name());
329 :
330 1054 : local_meshfuns.emplace_back(
331 527 : from_problem.es(), *from_sys.current_local_solution, from_sys.get_dof_map(), from_var_num);
332 527 : local_meshfuns.back().init();
333 527 : local_meshfuns.back().enable_out_of_mesh_mode(OutOfMeshValue);
334 : }
335 :
336 : // Recieve quadrature points from other processors, evaluate mesh frunctions
337 : // at those points, and send the values back.
338 499 : std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> outgoing_evals_ids;
339 :
340 : // If there is no cached data, we need to do communication
341 : // Quadrature points I will receive from remote processors
342 499 : std::map<processor_id_type, std::vector<Point>> incoming_qps;
343 499 : if (!_qps_cached)
344 : {
345 534 : auto qps_action_functor = [&incoming_qps](processor_id_type pid, const std::vector<Point> & qps)
346 : {
347 : // Quadrature points from processor 'pid'
348 534 : auto & incoming_qps_from_pid = incoming_qps[pid];
349 : // Store data for late use
350 534 : incoming_qps_from_pid.reserve(incoming_qps_from_pid.size() + qps.size());
351 534 : std::copy(qps.begin(), qps.end(), std::back_inserter(incoming_qps_from_pid));
352 985 : };
353 :
354 451 : Parallel::push_parallel_vector_data(comm(), outgoing_qps, qps_action_functor);
355 : }
356 :
357 : // Cache data
358 499 : if (!_qps_cached)
359 451 : _cached_qps = incoming_qps;
360 :
361 1081 : for (auto & qps : _cached_qps)
362 : {
363 582 : const processor_id_type pid = qps.first;
364 :
365 1164 : outgoing_evals_ids[pid].resize(qps.second.size(),
366 582 : std::make_pair(OutOfMeshValue, libMesh::invalid_uint));
367 :
368 55048 : for (unsigned int qp = 0; qp < qps.second.size(); qp++)
369 : {
370 54466 : Point qpt = qps.second[qp];
371 :
372 : // Loop until we've found the lowest-ranked app that actually contains
373 : // the quadrature point.
374 110192 : for (unsigned int i_from = 0; i_from < _from_problems.size(); i_from++)
375 : {
376 55726 : if (local_bboxes[i_from].contains_point(qpt))
377 : {
378 : const auto from_global_num =
379 54466 : _current_direction == TO_MULTIAPP ? 0 : _from_local2global_map[i_from];
380 54466 : outgoing_evals_ids[pid][qp].first =
381 54466 : (local_meshfuns[i_from])(_from_transforms[from_global_num]->mapBack(qpt));
382 54466 : if (_current_direction == FROM_MULTIAPP)
383 12338 : outgoing_evals_ids[pid][qp].second = from_global_num;
384 : }
385 : }
386 : }
387 : }
388 :
389 : ////////////////////
390 : // Gather all of the qp evaluations and pick out the best ones for each qp.
391 : ////////////////////
392 :
393 : // Values back from remote processors for my local quadrature points
394 499 : std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> incoming_evals_ids;
395 :
396 : auto evals_action_functor =
397 582 : [&incoming_evals_ids](processor_id_type pid,
398 582 : const std::vector<std::pair<Real, unsigned int>> & evals)
399 : {
400 : // evals for processor 'pid'
401 582 : auto & incoming_evals_ids_for_pid = incoming_evals_ids[pid];
402 : // Copy evals for late use
403 582 : incoming_evals_ids_for_pid.reserve(incoming_evals_ids_for_pid.size() + evals.size());
404 582 : std::copy(evals.begin(), evals.end(), std::back_inserter(incoming_evals_ids_for_pid));
405 1081 : };
406 :
407 499 : Parallel::push_parallel_vector_data(comm(), outgoing_evals_ids, evals_action_functor);
408 :
409 499 : std::vector<std::vector<Real>> final_evals(_to_problems.size());
410 499 : std::vector<std::map<dof_id_type, unsigned int>> trimmed_element_maps(_to_problems.size());
411 :
412 1127 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
413 : {
414 628 : MeshBase & to_mesh = _to_meshes[i_to]->getMesh();
415 628 : LinearImplicitSystem & system = *_proj_sys[i_to];
416 :
417 628 : FEType fe_type = system.variable_type(0);
418 628 : std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
419 628 : QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
420 :
421 41472 : for (const auto & elem : to_mesh.active_local_element_ptr_range())
422 : {
423 20422 : qrule.init(*elem);
424 :
425 20422 : bool element_is_evaled = false;
426 20422 : std::vector<Real> evals(qrule.n_points(), 0.);
427 :
428 85390 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
429 : {
430 64968 : unsigned int lowest_app_rank = libMesh::invalid_uint;
431 143030 : for (auto & values_ids : incoming_evals_ids)
432 : {
433 : // Current processor id
434 78062 : const processor_id_type pid = values_ids.first;
435 :
436 : // Ignore the selected processor if the element wasn't found in it's
437 : // bounding box.
438 : std::map<std::pair<unsigned int, unsigned int>, unsigned int> & map =
439 78062 : element_index_map[pid];
440 78062 : std::pair<unsigned int, unsigned int> key(i_to, elem->id());
441 78062 : if (map.find(key) == map.end())
442 24402 : continue;
443 54466 : unsigned int qp0 = map[key];
444 :
445 : // Ignore the selected processor if it's app has a higher rank than the
446 : // previously found lowest app rank.
447 54466 : if (_current_direction == FROM_MULTIAPP)
448 12338 : if (values_ids.second[qp0 + qp].second >= lowest_app_rank)
449 0 : continue;
450 :
451 : // Ignore the selected processor if the qp was actually outside the
452 : // processor's subapp's mesh.
453 54466 : if (values_ids.second[qp0 + qp].first == OutOfMeshValue)
454 806 : continue;
455 :
456 : // This is the best meshfunction evaluation so far, save it.
457 53660 : element_is_evaled = true;
458 53660 : evals[qp] = values_ids.second[qp0 + qp].first;
459 : }
460 : }
461 :
462 : // If we found good evaluations for any of the qps in this element, save
463 : // those evaluations for later.
464 20422 : if (element_is_evaled)
465 : {
466 16252 : trimmed_element_maps[i_to][elem->id()] = final_evals[i_to].size();
467 68320 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
468 52068 : final_evals[i_to].push_back(evals[qp]);
469 : }
470 21050 : }
471 628 : }
472 :
473 : ////////////////////
474 : // We now have just one or zero mesh function values at all of our local
475 : // quadrature points. Stash those values (and a map linking them to element
476 : // ids) in the equation systems parameters and project the solution.
477 : ////////////////////
478 :
479 1127 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
480 : {
481 628 : _to_es[i_to]->parameters.set<std::vector<Real> *>("final_evals") = &final_evals[i_to];
482 628 : _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>("element_map") =
483 628 : &trimmed_element_maps[i_to];
484 628 : projectSolution(i_to);
485 628 : _to_es[i_to]->parameters.set<std::vector<Real> *>("final_evals") = NULL;
486 628 : _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>("element_map") = NULL;
487 : }
488 :
489 499 : if (_fixed_meshes)
490 100 : _qps_cached = true;
491 :
492 499 : postExecute();
493 499 : }
494 :
495 : void
496 628 : MultiAppProjectionTransfer::projectSolution(unsigned int i_to)
497 : {
498 628 : FEProblemBase & to_problem = *_to_problems[i_to];
499 628 : EquationSystems & proj_es = to_problem.es();
500 628 : LinearImplicitSystem & ls = *_proj_sys[i_to];
501 : // activate the current transfer
502 628 : proj_es.parameters.set<MultiAppProjectionTransfer *>("transfer") = this;
503 :
504 : // TODO: specify solver params in an input file
505 : // solver tolerance
506 628 : Real tol = proj_es.parameters.get<Real>("linear solver tolerance");
507 628 : proj_es.parameters.set<Real>("linear solver tolerance") = 1e-10; // set our tolerance
508 : // solve it
509 628 : ls.solve();
510 628 : proj_es.parameters.set<Real>("linear solver tolerance") = tol; // restore the original tolerance
511 :
512 : // copy projected solution into target es
513 628 : MeshBase & to_mesh = proj_es.get_mesh();
514 :
515 628 : MooseVariableFEBase & to_var = to_problem.getVariable(
516 : 0, _to_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
517 628 : System & to_sys = to_var.sys().system();
518 628 : NumericVector<Number> * to_solution = to_sys.solution.get();
519 :
520 55120 : for (const auto & node : to_mesh.local_node_ptr_range())
521 : {
522 46766 : for (unsigned int comp = 0; comp < node->n_comp(to_sys.number(), to_var.number()); comp++)
523 : {
524 19520 : const dof_id_type proj_index = node->dof_number(ls.number(), _proj_var_num, comp);
525 19520 : const dof_id_type to_index = node->dof_number(to_sys.number(), to_var.number(), comp);
526 19520 : to_solution->set(to_index, (*ls.solution)(proj_index));
527 : }
528 628 : }
529 41472 : for (const auto & elem : to_mesh.active_local_element_ptr_range())
530 26870 : for (unsigned int comp = 0; comp < elem->n_comp(to_sys.number(), to_var.number()); comp++)
531 : {
532 6448 : const dof_id_type proj_index = elem->dof_number(ls.number(), _proj_var_num, comp);
533 6448 : const dof_id_type to_index = elem->dof_number(to_sys.number(), to_var.number(), comp);
534 6448 : to_solution->set(to_index, (*ls.solution)(proj_index));
535 628 : }
536 :
537 628 : to_solution->close();
538 628 : to_sys.update();
539 628 : }
|
