www.mooseframework.org
MultiAppMeshFunctionTransfer.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
16 
17 // MOOSE includes
18 #include "DisplacedProblem.h"
19 #include "FEProblem.h"
20 #include "MooseMesh.h"
21 #include "MooseTypes.h"
22 #include "MooseVariable.h"
23 
24 #include "libmesh/meshfree_interpolation.h"
25 #include "libmesh/system.h"
26 #include "libmesh/mesh_function.h"
27 #include "libmesh/mesh_tools.h"
28 #include "libmesh/parallel_algebra.h" // for communicator send and receive stuff
29 
30 template <>
33 {
35  params.addRequiredParam<std::vector<AuxVariableName>>(
36  "variable", "The auxiliary variable to store the transferred values in.");
37  params.addRequiredParam<std::vector<VariableName>>("source_variable",
38  "The variable to transfer from.");
39  params.addParam<bool>("displaced_source_mesh",
40  false,
41  "Whether or not to use the displaced mesh for the source mesh.");
42  params.addParam<bool>("displaced_target_mesh",
43  false,
44  "Whether or not to use the displaced mesh for the target mesh.");
45  params.addParam<bool>(
46  "error_on_miss",
47  false,
48  "Whether or not to error in the case that a target point is not found in the source domain.");
49  return params;
50 }
51 
53  : MultiAppTransfer(parameters),
54  _to_var_name(getParam<std::vector<AuxVariableName>>("variable")),
55  _from_var_name(getParam<std::vector<VariableName>>("source_variable")),
56  _error_on_miss(getParam<bool>("error_on_miss"))
57 {
58  _displaced_source_mesh = getParam<bool>("displaced_source_mesh");
59  _displaced_target_mesh = getParam<bool>("displaced_target_mesh");
60 
61  if (_to_var_name.size() == _from_var_name.size())
62  _var_size = _to_var_name.size();
63  else
64  mooseError("The number of variables to transfer to and from should be equal");
65 }
66 
67 void
69 {
70  for (unsigned int i = 0; i < _var_size; ++i)
71  if (_direction == TO_MULTIAPP)
73  else
75 }
76 
77 void
79 {
80  Moose::out << "Beginning MeshFunctionTransfer " << name() << std::endl;
81 
82  getAppInfo();
83 
84  _send_points.resize(_var_size);
85  _send_evals.resize(_var_size);
86  _send_ids.resize(_var_size);
87  // loop over the vector of variables and make the transfer one by one
88  for (unsigned int i = 0; i < _var_size; ++i)
90 
91  // Make sure all our sends succeeded.
92  for (unsigned int i = 0; i < _var_size; ++i)
93  for (processor_id_type i_proc = 0; i_proc < n_processors(); ++i_proc)
94  {
95  if (i_proc == processor_id())
96  continue;
97  _send_points[i][i_proc].wait();
98  _send_evals[i][i_proc].wait();
100  _send_ids[i][i_proc].wait();
101  }
102 
103  _console << "Finished MeshFunctionTransfer " << name() << std::endl;
104 }
105 
106 void
108 {
109  mooseAssert(i < _var_size, "The variable of index " << i << " does not exist");
110 
119  // Get the bounding boxes for the "from" domains.
120  std::vector<BoundingBox> bboxes = getFromBoundingBoxes();
121 
122  // Figure out how many "from" domains each processor owns.
123  std::vector<unsigned int> froms_per_proc = getFromsPerProc();
124 
125  std::vector<std::vector<Point>> outgoing_points(n_processors());
126  std::vector<std::map<std::pair<unsigned int, unsigned int>, unsigned int>> point_index_map(
127  n_processors());
128  // point_index_map[i_to, element_id] = index
129  // outgoing_points[index] is the first quadrature point in element
130 
131  for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
132  {
133  System * to_sys = find_sys(*_to_es[i_to], _to_var_name[i]);
134  unsigned int sys_num = to_sys->number();
135  unsigned int var_num = to_sys->variable_number(_to_var_name[i]);
136  MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
137  bool is_nodal = to_sys->variable_type(var_num).family == LAGRANGE;
138 
139  if (is_nodal)
140  {
141  MeshBase::const_node_iterator node_it = to_mesh->local_nodes_begin();
142  MeshBase::const_node_iterator node_end = to_mesh->local_nodes_end();
143 
144  for (; node_it != node_end; ++node_it)
145  {
146  Node * node = *node_it;
147 
148  // Skip this node if the variable has no dofs at it.
149  if (node->n_dofs(sys_num, var_num) < 1)
150  continue;
151 
152  // Loop over the "froms" on processor i_proc. If the node is found in
153  // any of the "froms", add that node to the vector that will be sent to
154  // i_proc.
155  unsigned int from0 = 0;
156  for (processor_id_type i_proc = 0; i_proc < n_processors();
157  from0 += froms_per_proc[i_proc], ++i_proc)
158  {
159  bool point_found = false;
160  for (unsigned int i_from = from0; i_from < from0 + froms_per_proc[i_proc] && !point_found;
161  ++i_from)
162  {
163  if (bboxes[i_from].contains_point(*node + _to_positions[i_to]))
164  {
165  std::pair<unsigned int, unsigned int> key(i_to, node->id());
166  point_index_map[i_proc][key] = outgoing_points[i_proc].size();
167  outgoing_points[i_proc].push_back(*node + _to_positions[i_to]);
168  point_found = true;
169  }
170  }
171  }
172  }
173  }
174  else // Elemental
175  {
176  MeshBase::const_element_iterator elem_it = to_mesh->local_elements_begin();
177  MeshBase::const_element_iterator elem_end = to_mesh->local_elements_end();
178 
179  for (; elem_it != elem_end; ++elem_it)
180  {
181  Elem * elem = *elem_it;
182 
183  Point centroid = elem->centroid();
184 
185  // Skip this element if the variable has no dofs at it.
186  if (elem->n_dofs(sys_num, var_num) < 1)
187  continue;
188 
189  // Loop over the "froms" on processor i_proc. If the elem is found in
190  // any of the "froms", add that elem to the vector that will be sent to
191  // i_proc.
192  unsigned int from0 = 0;
193  for (processor_id_type i_proc = 0; i_proc < n_processors();
194  from0 += froms_per_proc[i_proc], ++i_proc)
195  {
196  bool point_found = false;
197  for (unsigned int i_from = from0; i_from < from0 + froms_per_proc[i_proc] && !point_found;
198  ++i_from)
199  {
200  if (bboxes[i_from].contains_point(centroid + _to_positions[i_to]))
201  {
202  std::pair<unsigned int, unsigned int> key(i_to, elem->id());
203  point_index_map[i_proc][key] = outgoing_points[i_proc].size();
204  outgoing_points[i_proc].push_back(centroid + _to_positions[i_to]);
205  point_found = true;
206  }
207  }
208  }
209  }
210  }
211  }
212 
218  // Get the local bounding boxes.
219  std::vector<BoundingBox> local_bboxes(froms_per_proc[processor_id()]);
220  {
221  // Find the index to the first of this processor's local bounding boxes.
222  unsigned int local_start = 0;
223  for (processor_id_type i_proc = 0; i_proc < n_processors() && i_proc != processor_id();
224  ++i_proc)
225  {
226  local_start += froms_per_proc[i_proc];
227  }
228 
229  // Extract the local bounding boxes.
230  for (unsigned int i_from = 0; i_from < froms_per_proc[processor_id()]; ++i_from)
231  {
232  local_bboxes[i_from] = bboxes[local_start + i_from];
233  }
234  }
235 
236  // Setup the local mesh functions.
237  std::vector<std::shared_ptr<MeshFunction>> local_meshfuns;
238  for (unsigned int i_from = 0; i_from < _from_problems.size(); ++i_from)
239  {
240  FEProblemBase & from_problem = *_from_problems[i_from];
241  MooseVariable & from_var = from_problem.getVariable(0, _from_var_name[i]);
242  System & from_sys = from_var.sys().system();
243  unsigned int from_var_num = from_sys.variable_number(from_var.name());
244 
245  std::shared_ptr<MeshFunction> from_func;
246  // TODO: make MultiAppTransfer give me the right es
247  if (_displaced_source_mesh && from_problem.getDisplacedProblem())
248  from_func.reset(new MeshFunction(from_problem.getDisplacedProblem()->es(),
249  *from_sys.current_local_solution,
250  from_sys.get_dof_map(),
251  from_var_num));
252  else
253  from_func.reset(new MeshFunction(from_problem.es(),
254  *from_sys.current_local_solution,
255  from_sys.get_dof_map(),
256  from_var_num));
257  from_func->init(Trees::ELEMENTS);
258  from_func->enable_out_of_mesh_mode(OutOfMeshValue);
259  local_meshfuns.push_back(from_func);
260  }
261 
262  // Send points to other processors.
263  std::vector<std::vector<Real>> incoming_evals(n_processors());
264  std::vector<std::vector<unsigned int>> incoming_app_ids(n_processors());
265  _send_points[i].resize(n_processors());
266  for (processor_id_type i_proc = 0; i_proc < n_processors(); ++i_proc)
267  {
268  if (i_proc == processor_id())
269  continue;
270  _communicator.send(i_proc, outgoing_points[i_proc], _send_points[i][i_proc]);
271  }
272 
273  // Receive points from other processors, evaluate mesh functions at those
274  // points, and send the values back.
275  _send_evals[i].resize(n_processors());
276  _send_ids[i].resize(n_processors());
277 
278  // Create these here so that they live the entire life of this function
279  // and are NOT reused per processor.
280  std::vector<std::vector<Real>> processor_outgoing_evals(n_processors());
281 
282  for (processor_id_type i_proc = 0; i_proc < n_processors(); ++i_proc)
283  {
284  std::vector<Point> incoming_points;
285  if (i_proc == processor_id())
286  incoming_points = outgoing_points[i_proc];
287  else
288  _communicator.receive(i_proc, incoming_points);
289 
290  std::vector<Real> & outgoing_evals = processor_outgoing_evals[i_proc];
291  outgoing_evals.resize(incoming_points.size(), OutOfMeshValue);
292 
293  std::vector<unsigned int> outgoing_ids(incoming_points.size(), -1); // -1 = largest unsigned int
294  for (unsigned int i_pt = 0; i_pt < incoming_points.size(); ++i_pt)
295  {
296  Point pt = incoming_points[i_pt];
297 
298  // Loop until we've found the lowest-ranked app that actually contains
299  // the quadrature point.
300  for (unsigned int i_from = 0;
301  i_from < _from_problems.size() && outgoing_evals[i_pt] == OutOfMeshValue;
302  ++i_from)
303  {
304  if (local_bboxes[i_from].contains_point(pt))
305  {
306  outgoing_evals[i_pt] = (*local_meshfuns[i_from])(pt - _from_positions[i_from]);
307  if (_direction == FROM_MULTIAPP)
308  outgoing_ids[i_pt] = _local2global_map[i_from];
309  }
310  }
311  }
312 
313  if (i_proc == processor_id())
314  {
315  incoming_evals[i_proc] = outgoing_evals;
316  if (_direction == FROM_MULTIAPP)
317  incoming_app_ids[i_proc] = outgoing_ids;
318  }
319  else
320  {
321  _communicator.send(i_proc, outgoing_evals, _send_evals[i][i_proc]);
322  if (_direction == FROM_MULTIAPP)
323  _communicator.send(i_proc, outgoing_ids, _send_ids[i][i_proc]);
324  }
325  }
326 
334  for (processor_id_type i_proc = 0; i_proc < n_processors(); ++i_proc)
335  {
336  if (i_proc == processor_id())
337  continue;
338 
339  _communicator.receive(i_proc, incoming_evals[i_proc]);
340  if (_direction == FROM_MULTIAPP)
341  _communicator.receive(i_proc, incoming_app_ids[i_proc]);
342  }
343 
344  for (unsigned int i_to = 0; i_to < _to_problems.size(); ++i_to)
345  {
346  System * to_sys = find_sys(*_to_es[i_to], _to_var_name[i]);
347 
348  unsigned int sys_num = to_sys->number();
349  unsigned int var_num = to_sys->variable_number(_to_var_name[i]);
350 
351  NumericVector<Real> * solution = nullptr;
352  switch (_direction)
353  {
354  case TO_MULTIAPP:
355  solution = &getTransferVector(i_to, _to_var_name[i]);
356  break;
357  case FROM_MULTIAPP:
358  solution = to_sys->solution.get();
359  break;
360  default:
361  mooseError("Unknown direction");
362  }
363 
364  MeshBase * to_mesh = &_to_meshes[i_to]->getMesh();
365 
366  bool is_nodal = to_sys->variable_type(var_num).family == LAGRANGE;
367 
368  if (is_nodal)
369  {
370  MeshBase::const_node_iterator node_it = to_mesh->local_nodes_begin();
371  MeshBase::const_node_iterator node_end = to_mesh->local_nodes_end();
372 
373  for (; node_it != node_end; ++node_it)
374  {
375  Node * node = *node_it;
376 
377  // Skip this node if the variable has no dofs at it.
378  if (node->n_dofs(sys_num, var_num) < 1)
379  continue;
380 
381  unsigned int lowest_app_rank = libMesh::invalid_uint;
382  Real best_val = 0.;
383  bool point_found = false;
384  for (unsigned int i_proc = 0; i_proc < incoming_evals.size(); ++i_proc)
385  {
386  // Skip this proc if the node wasn't in it's bounding boxes.
387  std::pair<unsigned int, unsigned int> key(i_to, node->id());
388  if (point_index_map[i_proc].find(key) == point_index_map[i_proc].end())
389  continue;
390  unsigned int i_pt = point_index_map[i_proc][key];
391 
392  // Ignore this proc if it's app has a higher rank than the
393  // previously found lowest app rank.
394  if (_direction == FROM_MULTIAPP)
395  {
396  if (incoming_app_ids[i_proc][i_pt] >= lowest_app_rank)
397  continue;
398  }
399 
400  // Ignore this proc if the point was actually outside its meshes.
401  if (incoming_evals[i_proc][i_pt] == OutOfMeshValue)
402  continue;
403 
404  best_val = incoming_evals[i_proc][i_pt];
405  point_found = true;
406  }
407 
408  if (_error_on_miss && !point_found)
409  mooseError("Point not found! ", *node + _to_positions[i_to]);
410 
411  dof_id_type dof = node->dof_number(sys_num, var_num, 0);
412  solution->set(dof, best_val);
413  }
414  }
415  else // Elemental
416  {
417  MeshBase::const_element_iterator elem_it = to_mesh->local_elements_begin();
418  MeshBase::const_element_iterator elem_end = to_mesh->local_elements_end();
419 
420  for (; elem_it != elem_end; ++elem_it)
421  {
422  Elem * elem = *elem_it;
423 
424  // Skip this element if the variable has no dofs at it.
425  if (elem->n_dofs(sys_num, var_num) < 1)
426  continue;
427 
428  unsigned int lowest_app_rank = libMesh::invalid_uint;
429  Real best_val = 0;
430  bool point_found = false;
431  for (unsigned int i_proc = 0; i_proc < incoming_evals.size(); ++i_proc)
432  {
433  // Skip this proc if the elem wasn't in it's bounding boxes.
434  std::pair<unsigned int, unsigned int> key(i_to, elem->id());
435  if (point_index_map[i_proc].find(key) == point_index_map[i_proc].end())
436  continue;
437  unsigned int i_pt = point_index_map[i_proc][key];
438 
439  // Ignore this proc if it's app has a higher rank than the
440  // previously found lowest app rank.
441  if (_direction == FROM_MULTIAPP)
442  {
443  if (incoming_app_ids[i_proc][i_pt] >= lowest_app_rank)
444  continue;
445  }
446 
447  // Ignore this proc if the point was actually outside its meshes.
448  if (incoming_evals[i_proc][i_pt] == OutOfMeshValue)
449  continue;
450 
451  best_val = incoming_evals[i_proc][i_pt];
452  point_found = true;
453  }
454 
455  if (_error_on_miss && !point_found)
456  mooseError("Point not found! ", elem->centroid() + _to_positions[i_to]);
457 
458  dof_id_type dof = elem->dof_number(sys_num, var_num, 0);
459  solution->set(dof, best_val);
460  }
461  }
462  solution->close();
463  to_sys->update();
464  }
465 }
const std::string & name() const
Get the name of the object.
Definition: MooseObject.h:47
void transferVariable(unsigned int i)
Performs the transfer for the variable of index i.
NumericVector< Real > & getTransferVector(unsigned int i_local, std::string var_name)
If we are transferring to a multiapp, return the appropriate solution vector.
std::vector< std::vector< Parallel::Request > > _send_ids
To send app ids to other processors.
void variableIntegrityCheck(const AuxVariableName &var_name) const
Utility to verify that the vEariable in the destination system exists.
unsigned int _var_size
The number of variables to transfer.
Class for stuff related to variables.
Definition: MooseVariable.h:43
MultiAppMeshFunctionTransfer(const InputParameters &parameters)
std::vector< std::vector< Parallel::Request > > _send_points
To send points to other processors.
std::vector< EquationSystems * > _to_es
std::vector< BoundingBox > getFromBoundingBoxes()
Return the bounding boxes of all the "from" domains, including all the domains not local to this proc...
std::vector< FEProblemBase * > _to_problems
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
std::vector< Point > _from_positions
InputParameters validParams< MultiAppMeshFunctionTransfer >()
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
const std::string & name() const
Get the variable name.
virtual EquationSystems & es() override
std::vector< std::vector< Parallel::Request > > _send_evals
To send values to other processors.
std::vector< AuxVariableName > _to_var_name
The vector of variables to transfer to.
virtual std::shared_ptr< DisplacedProblem > getDisplacedProblem()
MooseEnum _direction
Whether we&#39;re transferring to or from the MultiApp.
std::vector< unsigned int > getFromsPerProc()
Return the number of "from" domains that each processor owns.
virtual System & system()=0
Get the reference to the libMesh system.
std::vector< unsigned int > _local2global_map
static const Number OutOfMeshValue
Definition: Transfer.h:81
virtual void initialSetup() override
Method called at the beginning of the simulation for checking integrity or doing one-time setup...
std::vector< VariableName > _from_var_name
The vector of variables to transfer from.
std::vector< Point > _to_positions
InputParameters validParams< MultiAppTransfer >()
Base class for all MultiAppTransfer objects.
virtual void execute() override
Execute the transfer.
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
void mooseError(Args &&...args) const
Definition: MooseObject.h:80
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
static System * find_sys(EquationSystems &es, const std::string &var_name)
Small helper function for finding the system containing the variable.
Definition: Transfer.C:70
void getAppInfo()
This method will fill information into the convenience member variables (_to_problems, _from_meshes, etc.)
SystemBase & sys()
Get the system this variable is part of.
virtual MooseVariable & getVariable(THREAD_ID tid, const std::string &var_name) override
Returns the variable reference for requested variable which may be in any system. ...
std::vector< FEProblemBase * > _from_problems
std::vector< MooseMesh * > _to_meshes