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 566 : assemble_l2(EquationSystems & es, const std::string & system_name)
37 : {
38 : MultiAppProjectionTransfer * transfer =
39 566 : es.parameters.get<MultiAppProjectionTransfer *>("transfer");
40 566 : transfer->assembleL2(es, system_name);
41 566 : }
42 :
43 : registerMooseObject("MooseApp", MultiAppProjectionTransfer);
44 :
45 : InputParameters
46 3613 : MultiAppProjectionTransfer::validParams()
47 : {
48 3613 : InputParameters params = MultiAppConservativeTransfer::validParams();
49 7226 : params.addClassDescription(
50 : "Perform a projection between a master and sub-application mesh of a field variable.");
51 :
52 14452 : MooseEnum proj_type("l2", "l2");
53 14452 : params.addParam<MooseEnum>("proj_type", proj_type, "The type of the projection.");
54 :
55 7226 : params.addParam<bool>("fixed_meshes",
56 7226 : 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 10839 : params.addRelationshipManager("ElementSideNeighborLayers",
63 : Moose::RelationshipManagerType::GEOMETRIC |
64 : Moose::RelationshipManagerType::ALGEBRAIC);
65 3613 : MultiAppTransfer::addBBoxFactorParam(params);
66 7226 : return params;
67 3613 : }
68 :
69 276 : MultiAppProjectionTransfer::MultiAppProjectionTransfer(const InputParameters & parameters)
70 : : MultiAppConservativeTransfer(parameters),
71 276 : _proj_type(getParam<MooseEnum>("proj_type")),
72 276 : _compute_matrix(true),
73 552 : _fixed_meshes(getParam<bool>("fixed_meshes")),
74 828 : _qps_cached(false)
75 : {
76 276 : if (_to_var_names.size() != 1)
77 0 : paramError("variable", " Support single to-variable only ");
78 :
79 276 : if (_from_var_names.size() != 1)
80 0 : paramError("source_variable", " Support single from-variable only ");
81 276 : }
82 :
83 : void
84 276 : MultiAppProjectionTransfer::initialSetup()
85 : {
86 276 : MultiAppConservativeTransfer::initialSetup();
87 :
88 276 : _proj_sys.resize(_to_problems.size(), NULL);
89 :
90 596 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
91 : {
92 320 : FEProblemBase & to_problem = *_to_problems[i_to];
93 320 : EquationSystems & to_es = to_problem.es();
94 :
95 : // Add the projection system.
96 : FEType fe_type = to_problem
97 320 : .getVariable(0,
98 : _to_var_name,
99 : Moose::VarKindType::VAR_ANY,
100 : Moose::VarFieldType::VAR_FIELD_STANDARD)
101 320 : .feType();
102 :
103 320 : LinearImplicitSystem & proj_sys = to_es.add_system<LinearImplicitSystem>("proj-sys-" + name());
104 :
105 320 : _proj_var_num = proj_sys.add_variable("var", fe_type);
106 320 : proj_sys.attach_assemble_function(assemble_l2);
107 320 : _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 320 : proj_sys.hide_output() = true;
118 :
119 : // Reinitialize EquationSystems since we added a system.
120 320 : to_es.reinit();
121 : }
122 276 : }
123 :
124 : void
125 566 : MultiAppProjectionTransfer::assembleL2(EquationSystems & es, const std::string & system_name)
126 : {
127 : // Get the system and mesh from the input arguments.
128 566 : LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);
129 566 : MeshBase & to_mesh = es.get_mesh();
130 :
131 : // Get the meshfunction evaluations and the map that was stashed in the es.
132 566 : std::vector<Real> & final_evals = *es.parameters.get<std::vector<Real> *>("final_evals");
133 : std::map<dof_id_type, unsigned int> & element_map =
134 566 : *es.parameters.get<std::map<dof_id_type, unsigned int> *>("element_map");
135 :
136 : // Setup system vectors and matrices.
137 566 : FEType fe_type = system.variable_type(0);
138 566 : std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
139 566 : QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
140 566 : fe->attach_quadrature_rule(&qrule);
141 566 : const DofMap & dof_map = system.get_dof_map();
142 566 : DenseMatrix<Number> Ke;
143 566 : DenseVector<Number> Fe;
144 566 : std::vector<dof_id_type> dof_indices;
145 566 : const std::vector<Real> & JxW = fe->get_JxW();
146 566 : const std::vector<std::vector<Real>> & phi = fe->get_phi();
147 566 : auto & system_matrix = system.get_system_matrix();
148 :
149 36870 : for (const auto & elem : to_mesh.active_local_element_ptr_range())
150 : {
151 18152 : fe->reinit(elem);
152 :
153 18152 : dof_map.dof_indices(elem, dof_indices);
154 18152 : Ke.resize(dof_indices.size(), dof_indices.size());
155 18152 : Fe.resize(dof_indices.size());
156 :
157 75794 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
158 : {
159 57642 : Real meshfun_eval = 0.;
160 57642 : if (element_map.find(elem->id()) != element_map.end())
161 : {
162 : // We have evaluations for this element.
163 46172 : meshfun_eval = final_evals[element_map[elem->id()] + qp];
164 : }
165 :
166 : // Now compute the element matrix and RHS contributions.
167 274588 : for (unsigned int i = 0; i < phi.size(); i++)
168 : {
169 : // RHS
170 216946 : Fe(i) += JxW[qp] * (meshfun_eval * phi[i][qp]);
171 :
172 216946 : if (_compute_matrix)
173 1089476 : for (unsigned int j = 0; j < phi.size(); j++)
174 : {
175 : // The matrix contribution
176 872530 : Ke(i, j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
177 : }
178 : }
179 57642 : dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
180 :
181 57642 : if (_compute_matrix)
182 57642 : system_matrix.add_matrix(Ke, dof_indices);
183 57642 : system.rhs->add_vector(Fe, dof_indices);
184 : }
185 566 : }
186 566 : }
187 :
188 : void
189 458 : MultiAppProjectionTransfer::execute()
190 : {
191 2290 : 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 458 : std::vector<BoundingBox> bboxes = getFromBoundingBoxes();
229 :
230 : // Figure out how many "from" domains each processor owns.
231 458 : std::vector<unsigned int> froms_per_proc = getFromsPerProc();
232 :
233 458 : 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 458 : 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 458 : if (!_qps_cached)
240 : {
241 936 : 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 522 : _current_direction == FROM_MULTIAPP ? 0 : _to_local2global_map[i_to];
246 522 : MeshBase & to_mesh = _to_meshes[i_to]->getMesh();
247 :
248 522 : LinearImplicitSystem & system = *_proj_sys[i_to];
249 :
250 522 : FEType fe_type = system.variable_type(0);
251 522 : std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
252 522 : QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
253 522 : fe->attach_quadrature_rule(&qrule);
254 522 : const std::vector<Point> & xyz = fe->get_xyz();
255 :
256 522 : unsigned int from0 = 0;
257 1266 : for (processor_id_type i_proc = 0; i_proc < n_processors();
258 744 : from0 += froms_per_proc[i_proc], i_proc++)
259 : {
260 744 : for (const auto & elem :
261 24250 : as_range(to_mesh.local_elements_begin(), to_mesh.local_elements_end()))
262 : {
263 22762 : fe->reinit(elem);
264 :
265 22762 : bool qp_hit = false;
266 47252 : for (unsigned int i_from = 0; i_from < froms_per_proc[i_proc] && !qp_hit; i_from++)
267 : {
268 71867 : for (unsigned int qp = 0; qp < qrule.n_points() && !qp_hit; qp++)
269 : {
270 47377 : Point qpt = xyz[qp];
271 47377 : if (bboxes[from0 + i_from].contains_point((*_to_transforms[to_global_num])(qpt)))
272 13507 : qp_hit = true;
273 : }
274 : }
275 :
276 22762 : 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 13507 : std::pair<unsigned int, unsigned int> key(i_to, elem->id());
282 13507 : element_index_map[i_proc][key] = outgoing_qps[i_proc].size();
283 58077 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
284 : {
285 44570 : Point qpt = xyz[qp];
286 44570 : outgoing_qps[i_proc].push_back((*_to_transforms[to_global_num])(qpt));
287 : }
288 : }
289 744 : }
290 : }
291 522 : }
292 :
293 414 : if (_fixed_meshes)
294 48 : _cached_index_map = element_index_map;
295 : }
296 : else
297 : {
298 44 : 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 458 : 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 458 : unsigned int local_start = 0;
311 581 : 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 940 : for (unsigned int i_from = 0; i_from < froms_per_proc[processor_id()]; i_from++)
317 482 : local_bboxes[i_from] = bboxes[local_start + i_from];
318 : }
319 :
320 : // Setup the local mesh functions.
321 458 : std::vector<MeshFunction> local_meshfuns;
322 940 : for (unsigned int i_from = 0; i_from < _from_problems.size(); i_from++)
323 : {
324 482 : FEProblemBase & from_problem = *_from_problems[i_from];
325 482 : MooseVariableFEBase & from_var = from_problem.getVariable(
326 : 0, _from_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
327 482 : System & from_sys = from_var.sys().system();
328 482 : unsigned int from_var_num = from_sys.variable_number(from_var.name());
329 :
330 964 : local_meshfuns.emplace_back(
331 482 : from_problem.es(), *from_sys.current_local_solution, from_sys.get_dof_map(), from_var_num);
332 482 : local_meshfuns.back().init();
333 482 : 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 458 : 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 458 : std::map<processor_id_type, std::vector<Point>> incoming_qps;
343 458 : if (!_qps_cached)
344 : {
345 497 : auto qps_action_functor = [&incoming_qps](processor_id_type pid, const std::vector<Point> & qps)
346 : {
347 : // Quadrature points from processor 'pid'
348 497 : auto & incoming_qps_from_pid = incoming_qps[pid];
349 : // Store data for late use
350 497 : incoming_qps_from_pid.reserve(incoming_qps_from_pid.size() + qps.size());
351 497 : std::copy(qps.begin(), qps.end(), std::back_inserter(incoming_qps_from_pid));
352 911 : };
353 :
354 414 : Parallel::push_parallel_vector_data(comm(), outgoing_qps, qps_action_functor);
355 : }
356 :
357 : // Cache data
358 458 : if (!_qps_cached)
359 414 : _cached_qps = incoming_qps;
360 :
361 999 : for (auto & qps : _cached_qps)
362 : {
363 541 : const processor_id_type pid = qps.first;
364 :
365 1082 : outgoing_evals_ids[pid].resize(qps.second.size(),
366 541 : std::make_pair(OutOfMeshValue, libMesh::invalid_uint));
367 :
368 49111 : for (unsigned int qp = 0; qp < qps.second.size(); qp++)
369 : {
370 48570 : Point qpt = qps.second[qp];
371 :
372 : // Loop until we've found the lowest-ranked app that actually contains
373 : // the quadrature point.
374 98220 : for (unsigned int i_from = 0; i_from < _from_problems.size(); i_from++)
375 : {
376 49650 : if (local_bboxes[i_from].contains_point(qpt))
377 : {
378 97140 : outgoing_evals_ids[pid][qp].first = (local_meshfuns[i_from])(
379 145710 : getPointInSourceAppFrame(qpt, i_from, "Projection transfer evaluation"));
380 48570 : if (_current_direction == FROM_MULTIAPP)
381 11022 : outgoing_evals_ids[pid][qp].second = getGlobalSourceAppIndex(i_from);
382 : }
383 : }
384 : }
385 : }
386 :
387 : ////////////////////
388 : // Gather all of the qp evaluations and pick out the best ones for each qp.
389 : ////////////////////
390 :
391 : // Values back from remote processors for my local quadrature points
392 458 : std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> incoming_evals_ids;
393 :
394 : auto evals_action_functor =
395 541 : [&incoming_evals_ids](processor_id_type pid,
396 : const std::vector<std::pair<Real, unsigned int>> & evals)
397 : {
398 : // evals for processor 'pid'
399 541 : auto & incoming_evals_ids_for_pid = incoming_evals_ids[pid];
400 : // Copy evals for late use
401 541 : incoming_evals_ids_for_pid.reserve(incoming_evals_ids_for_pid.size() + evals.size());
402 541 : std::copy(evals.begin(), evals.end(), std::back_inserter(incoming_evals_ids_for_pid));
403 999 : };
404 :
405 458 : Parallel::push_parallel_vector_data(comm(), outgoing_evals_ids, evals_action_functor);
406 :
407 916 : std::vector<std::vector<Real>> final_evals(_to_problems.size());
408 458 : std::vector<std::map<dof_id_type, unsigned int>> trimmed_element_maps(_to_problems.size());
409 :
410 1024 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
411 : {
412 566 : MeshBase & to_mesh = _to_meshes[i_to]->getMesh();
413 566 : LinearImplicitSystem & system = *_proj_sys[i_to];
414 :
415 566 : FEType fe_type = system.variable_type(0);
416 566 : std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
417 566 : QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
418 :
419 18718 : for (const auto & elem : to_mesh.active_local_element_ptr_range())
420 : {
421 18152 : qrule.init(*elem);
422 :
423 18152 : bool element_is_evaled = false;
424 18152 : std::vector<Real> evals(qrule.n_points(), 0.);
425 :
426 75794 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
427 : {
428 57642 : unsigned int lowest_app_rank = libMesh::invalid_uint;
429 128378 : for (auto & values_ids : incoming_evals_ids)
430 : {
431 : // Current processor id
432 70736 : const processor_id_type pid = values_ids.first;
433 :
434 : // Ignore the selected processor if the element wasn't found in it's
435 : // bounding box.
436 : std::map<std::pair<unsigned int, unsigned int>, unsigned int> & map =
437 70736 : element_index_map[pid];
438 70736 : std::pair<unsigned int, unsigned int> key(i_to, elem->id());
439 70736 : if (map.find(key) == map.end())
440 23042 : continue;
441 48570 : unsigned int qp0 = map[key];
442 :
443 : // Ignore the selected processor if it's app has a higher rank than the
444 : // previously found lowest app rank.
445 48570 : if (_current_direction == FROM_MULTIAPP)
446 11022 : if (values_ids.second[qp0 + qp].second >= lowest_app_rank)
447 0 : continue;
448 :
449 : // Ignore the selected processor if the qp was actually outside the
450 : // processor's subapp's mesh.
451 48570 : if (values_ids.second[qp0 + qp].first == OutOfMeshValue)
452 876 : continue;
453 :
454 : // This is the best meshfunction evaluation so far, save it.
455 47694 : element_is_evaled = true;
456 47694 : evals[qp] = values_ids.second[qp0 + qp].first;
457 : }
458 : }
459 :
460 : // If we found good evaluations for any of the qps in this element, save
461 : // those evaluations for later.
462 18152 : if (element_is_evaled)
463 : {
464 14434 : trimmed_element_maps[i_to][elem->id()] = final_evals[i_to].size();
465 60606 : for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
466 46172 : final_evals[i_to].push_back(evals[qp]);
467 : }
468 18718 : }
469 566 : }
470 :
471 : ////////////////////
472 : // We now have just one or zero mesh function values at all of our local
473 : // quadrature points. Stash those values (and a map linking them to element
474 : // ids) in the equation systems parameters and project the solution.
475 : ////////////////////
476 :
477 1024 : for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
478 : {
479 1132 : _to_es[i_to]->parameters.set<std::vector<Real> *>("final_evals") = &final_evals[i_to];
480 1132 : _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>("element_map") =
481 566 : &trimmed_element_maps[i_to];
482 566 : projectSolution(i_to);
483 1132 : _to_es[i_to]->parameters.set<std::vector<Real> *>("final_evals") = NULL;
484 1698 : _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>("element_map") = NULL;
485 : }
486 :
487 458 : if (_fixed_meshes)
488 92 : _qps_cached = true;
489 :
490 458 : postExecute();
491 458 : }
492 :
493 : void
494 566 : MultiAppProjectionTransfer::projectSolution(unsigned int i_to)
495 : {
496 566 : FEProblemBase & to_problem = *_to_problems[i_to];
497 566 : EquationSystems & proj_es = to_problem.es();
498 566 : LinearImplicitSystem & ls = *_proj_sys[i_to];
499 : // activate the current transfer
500 1132 : proj_es.parameters.set<MultiAppProjectionTransfer *>("transfer") = this;
501 :
502 : // TODO: specify solver params in an input file
503 : // solver tolerance
504 566 : Real tol = proj_es.parameters.get<Real>("linear solver tolerance");
505 1132 : proj_es.parameters.set<Real>("linear solver tolerance") = 1e-10; // set our tolerance
506 : // solve it
507 566 : ls.solve();
508 1132 : proj_es.parameters.set<Real>("linear solver tolerance") = tol; // restore the original tolerance
509 :
510 : // copy projected solution into target es
511 566 : MeshBase & to_mesh = proj_es.get_mesh();
512 :
513 566 : MooseVariableFEBase & to_var = to_problem.getVariable(
514 : 0, _to_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
515 566 : System & to_sys = to_var.sys().system();
516 566 : NumericVector<Number> * to_solution = to_sys.solution.get();
517 :
518 24788 : for (const auto & node : to_mesh.local_node_ptr_range())
519 : {
520 41538 : for (unsigned int comp = 0; comp < node->n_comp(to_sys.number(), to_var.number()); comp++)
521 : {
522 17316 : const dof_id_type proj_index = node->dof_number(ls.number(), _proj_var_num, comp);
523 17316 : const dof_id_type to_index = node->dof_number(to_sys.number(), to_var.number(), comp);
524 17316 : to_solution->set(to_index, (*ls.solution)(proj_index));
525 : }
526 566 : }
527 18718 : for (const auto & elem : to_mesh.active_local_element_ptr_range())
528 23906 : for (unsigned int comp = 0; comp < elem->n_comp(to_sys.number(), to_var.number()); comp++)
529 : {
530 5754 : const dof_id_type proj_index = elem->dof_number(ls.number(), _proj_var_num, comp);
531 5754 : const dof_id_type to_index = elem->dof_number(to_sys.number(), to_var.number(), comp);
532 5754 : to_solution->set(to_index, (*ls.solution)(proj_index));
533 566 : }
534 :
535 566 : to_solution->close();
536 566 : to_sys.update();
537 566 : }
|
