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 #include "libmesh/default_coupling.h" 35 assemble_l2(EquationSystems & es,
const std::string & system_name)
49 "Perform a projection between a master and sub-application mesh of a field variable.");
52 params.
addParam<
MooseEnum>(
"proj_type", proj_type,
"The type of the projection.");
54 params.
addParam<
bool>(
"fixed_meshes",
56 "Set to true when the meshes are not changing (ie, " 57 "no movement or adaptivity). This will cache some " 58 "information to speed up the transfer.");
70 _proj_type(getParam<
MooseEnum>(
"proj_type")),
71 _compute_matrix(true),
72 _fixed_meshes(getParam<bool>(
"fixed_meshes")),
76 paramError(
"variable",
" Support single to-variable only ");
79 paramError(
"source_variable",
" Support single from-variable only ");
89 for (
unsigned int i_to = 0; i_to <
_to_problems.size(); i_to++)
92 EquationSystems & to_es = to_problem.
es();
95 FEType fe_type = to_problem
102 LinearImplicitSystem & proj_sys = to_es.add_system<LinearImplicitSystem>(
"proj-sys-" +
name());
104 proj_sys.get_dof_map().add_coupling_functor(
105 proj_sys.get_dof_map().default_coupling(),
120 proj_sys.hide_output() =
true;
131 LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);
132 MeshBase & to_mesh = es.get_mesh();
135 std::vector<Real> & final_evals = *es.parameters.get<std::vector<Real> *>(
"final_evals");
136 std::map<dof_id_type, unsigned int> & element_map =
137 *es.parameters.get<std::map<dof_id_type, unsigned int> *>(
"element_map");
140 FEType fe_type = system.variable_type(0);
141 std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
142 QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
143 fe->attach_quadrature_rule(&qrule);
144 const DofMap & dof_map = system.get_dof_map();
147 std::vector<dof_id_type> dof_indices;
148 const std::vector<Real> & JxW = fe->get_JxW();
149 const std::vector<std::vector<Real>> & phi = fe->get_phi();
150 auto & system_matrix = system.get_system_matrix();
152 for (
const auto & elem : to_mesh.active_local_element_ptr_range())
156 dof_map.dof_indices(elem, dof_indices);
157 Ke.
resize(dof_indices.size(), dof_indices.size());
158 Fe.
resize(dof_indices.size());
160 for (
unsigned int qp = 0; qp < qrule.n_points(); qp++)
162 Real meshfun_eval = 0.;
163 if (element_map.find(elem->id()) != element_map.end())
166 meshfun_eval = final_evals[element_map[elem->id()] + qp];
170 for (
unsigned int i = 0; i < phi.size(); i++)
173 Fe(i) += JxW[qp] * (meshfun_eval * phi[i][qp]);
176 for (
unsigned int j = 0; j < phi.size(); j++)
179 Ke(i, j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
182 dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
185 system_matrix.add_matrix(Ke, dof_indices);
186 system.rhs->add_vector(Fe, dof_indices);
195 "MultiAppProjectionTransfer::execute()", 5,
"Transferring variables through projection");
236 std::map<processor_id_type, std::vector<Point>> outgoing_qps;
237 std::map<processor_id_type, std::map<std::pair<unsigned int, unsigned int>,
unsigned int>>
244 for (
unsigned int i_to = 0; i_to <
_to_problems.size(); i_to++)
247 const auto to_global_num =
249 MeshBase & to_mesh =
_to_meshes[i_to]->getMesh();
251 LinearImplicitSystem & system = *
_proj_sys[i_to];
253 FEType fe_type = system.variable_type(0);
254 std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
255 QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
256 fe->attach_quadrature_rule(&qrule);
257 const std::vector<Point> & xyz = fe->get_xyz();
259 unsigned int from0 = 0;
261 from0 += froms_per_proc[i_proc], i_proc++)
263 for (
const auto & elem :
264 as_range(to_mesh.local_elements_begin(), to_mesh.local_elements_end()))
269 for (
unsigned int i_from = 0; i_from < froms_per_proc[i_proc] && !qp_hit; i_from++)
271 for (
unsigned int qp = 0; qp < qrule.n_points() && !qp_hit; qp++)
274 if (bboxes[from0 + i_from].contains_point((*
_to_transforms[to_global_num])(qpt)))
284 std::pair<unsigned int, unsigned int> key(i_to, elem->id());
285 element_index_map[i_proc][key] = outgoing_qps[i_proc].size();
286 for (
unsigned int qp = 0; qp < qrule.n_points(); qp++)
289 outgoing_qps[i_proc].push_back((*
_to_transforms[to_global_num])(qpt));
310 std::vector<BoundingBox> local_bboxes(froms_per_proc[
processor_id()]);
313 unsigned int local_start = 0;
316 local_start += froms_per_proc[i_proc];
319 for (
unsigned int i_from = 0; i_from < froms_per_proc[
processor_id()]; i_from++)
320 local_bboxes[i_from] = bboxes[local_start + i_from];
324 std::vector<MeshFunction> local_meshfuns;
325 for (
unsigned int i_from = 0; i_from <
_from_problems.size(); i_from++)
330 System & from_sys = from_var.
sys().
system();
331 unsigned int from_var_num = from_sys.variable_number(from_var.
name());
333 local_meshfuns.emplace_back(
334 from_problem.
es(), *from_sys.current_local_solution, from_sys.get_dof_map(), from_var_num);
335 local_meshfuns.back().init();
341 std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> outgoing_evals_ids;
345 std::map<processor_id_type, std::vector<Point>> incoming_qps;
348 auto qps_action_functor = [&incoming_qps](
processor_id_type pid,
const std::vector<Point> & qps)
351 auto & incoming_qps_from_pid = incoming_qps[pid];
353 incoming_qps_from_pid.reserve(incoming_qps_from_pid.size() + qps.size());
354 std::copy(qps.begin(), qps.end(), std::back_inserter(incoming_qps_from_pid));
357 Parallel::push_parallel_vector_data(
comm(), outgoing_qps, qps_action_functor);
368 outgoing_evals_ids[pid].resize(qps.second.size(),
371 for (
unsigned int qp = 0; qp < qps.second.size(); qp++)
373 Point qpt = qps.second[qp];
377 for (
unsigned int i_from = 0; i_from <
_from_problems.size(); i_from++)
379 if (local_bboxes[i_from].contains_point(qpt))
381 const auto from_global_num =
383 outgoing_evals_ids[pid][qp].first =
386 outgoing_evals_ids[pid][qp].second = from_global_num;
397 std::map<processor_id_type, std::vector<std::pair<Real, unsigned int>>> incoming_evals_ids;
399 auto evals_action_functor =
401 const std::vector<std::pair<Real, unsigned int>> & evals)
404 auto & incoming_evals_ids_for_pid = incoming_evals_ids[pid];
406 incoming_evals_ids_for_pid.reserve(incoming_evals_ids_for_pid.size() + evals.size());
407 std::copy(evals.begin(), evals.end(), std::back_inserter(incoming_evals_ids_for_pid));
410 Parallel::push_parallel_vector_data(
comm(), outgoing_evals_ids, evals_action_functor);
412 std::vector<std::vector<Real>> final_evals(
_to_problems.size());
413 std::vector<std::map<dof_id_type, unsigned int>> trimmed_element_maps(
_to_problems.size());
415 for (
unsigned int i_to = 0; i_to <
_to_problems.size(); i_to++)
417 MeshBase & to_mesh =
_to_meshes[i_to]->getMesh();
418 LinearImplicitSystem & system = *
_proj_sys[i_to];
420 FEType fe_type = system.variable_type(0);
421 std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
422 QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
423 fe->attach_quadrature_rule(&qrule);
425 for (
const auto & elem : to_mesh.active_local_element_ptr_range())
429 bool element_is_evaled =
false;
430 std::vector<Real> evals(qrule.n_points(), 0.);
432 for (
unsigned int qp = 0; qp < qrule.n_points(); qp++)
435 for (
auto & values_ids : incoming_evals_ids)
442 std::map<std::pair<unsigned int, unsigned int>,
unsigned int> & map =
443 element_index_map[pid];
444 std::pair<unsigned int, unsigned int> key(i_to, elem->id());
445 if (map.find(key) == map.end())
447 unsigned int qp0 = map[key];
452 if (values_ids.second[qp0 + qp].second >= lowest_app_rank)
461 element_is_evaled =
true;
462 evals[qp] = values_ids.second[qp0 + qp].first;
468 if (element_is_evaled)
470 trimmed_element_maps[i_to][elem->id()] = final_evals[i_to].size();
471 for (
unsigned int qp = 0; qp < qrule.n_points(); qp++)
472 final_evals[i_to].push_back(evals[qp]);
483 for (
unsigned int i_to = 0; i_to <
_to_problems.size(); i_to++)
485 _to_es[i_to]->parameters.set<std::vector<Real> *>(
"final_evals") = &final_evals[i_to];
486 _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>(
"element_map") =
487 &trimmed_element_maps[i_to];
489 _to_es[i_to]->parameters.set<std::vector<Real> *>(
"final_evals") = NULL;
490 _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>(
"element_map") = NULL;
503 EquationSystems & proj_es = to_problem.
es();
504 LinearImplicitSystem & ls = *
_proj_sys[i_to];
510 Real tol = proj_es.parameters.get<
Real>(
"linear solver tolerance");
511 proj_es.parameters.set<
Real>(
"linear solver tolerance") = 1e-10;
514 proj_es.parameters.set<
Real>(
"linear solver tolerance") = tol;
517 MeshBase & to_mesh = proj_es.get_mesh();
524 for (
const auto & node : to_mesh.local_node_ptr_range())
526 for (
unsigned int comp = 0; comp < node->n_comp(to_sys.number(), to_var.
number()); comp++)
529 const dof_id_type to_index = node->dof_number(to_sys.number(), to_var.
number(), comp);
530 to_solution->
set(to_index, (*ls.solution)(proj_index));
533 for (
const auto & elem : to_mesh.active_local_element_ptr_range())
534 for (
unsigned int comp = 0; comp < elem->n_comp(to_sys.number(), to_var.
number()); comp++)
537 const dof_id_type to_index = elem->dof_number(to_sys.number(), to_var.
number(), comp);
538 to_solution->
set(to_index, (*ls.solution)(proj_index));
541 to_solution->
close();
std::vector< std::unique_ptr< MultiAppCoordTransform > > _to_transforms
void assemble_l2(EquationSystems &es, const std::string &system_name)
const unsigned int invalid_uint
MooseEnum _current_direction
virtual void get(const std::vector< numeric_index_type > &index, Number *values) const
unsigned int number() const
Get variable number coming from libMesh.
AuxVariableName _to_var_name
std::map< processor_id_type, std::map< std::pair< unsigned int, unsigned int >, unsigned int > > _cached_index_map
std::vector< EquationSystems * > _to_es
void resize(const unsigned int n)
std::map< processor_id_type, std::vector< Point > > _cached_qps
std::vector< BoundingBox > getFromBoundingBoxes()
Return the bounding boxes of all the "from" domains, including all the domains not local to this proc...
const std::vector< VariableName > _from_var_names
Name of variables transferring from.
std::vector< FEProblemBase * > _to_problems
const std::string & name() const override
Get the variable name.
const Parallel::Communicator & comm() const
virtual void postExecute()
Add some extra work if necessary after execute().
This class provides an interface for common operations on field variables of both FE and FV types wit...
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
void assembleL2(EquationSystems &es, const std::string &system_name)
virtual const std::string & name() const
Get the name of the class.
registerMooseObject("MooseApp", MultiAppProjectionTransfer)
const FEType & feType() const
Get the type of finite element object.
virtual const MooseVariableFieldBase & getVariable(const THREAD_ID tid, const std::string &var_name, Moose::VarKindType expected_var_type=Moose::VarKindType::VAR_ANY, Moose::VarFieldType expected_var_field_type=Moose::VarFieldType::VAR_FIELD_ANY) const override
Returns the variable reference for requested variable which must be of the expected_var_type (Nonline...
std::vector< unsigned int > _to_local2global_map
Given local app index, returns global app index.
virtual EquationSystems & es() override
virtual void initialSetup() override
Method called at the beginning of the simulation for checking integrity or doing one-time setup...
uint8_t processor_id_type
processor_id_type n_processors() const
Project values from one domain to another.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
bool _compute_matrix
True, if we need to recompute the projection matrix.
const std::vector< AuxVariableName > _to_var_names
Name of variables transferring to.
std::vector< unsigned int > getFromsPerProc()
Return the number of "from" domains that each processor owns.
Transfers variables on possibly different meshes while conserving a user defined property (Postproces...
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
static InputParameters validParams()
virtual void execute() override
Execute the transfer.
virtual System & system()=0
Get the reference to the libMesh system.
void paramError(const std::string ¶m, Args... args) const
Emits an error prefixed with the file and line number of the given param (from the input file) along ...
static const Number OutOfMeshValue
MultiAppProjectionTransfer(const InputParameters ¶meters)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void resize(const unsigned int new_m, const unsigned int new_n)
std::vector< std::unique_ptr< MultiAppCoordTransform > > _from_transforms
std::vector< unsigned int > _from_local2global_map
Given local app index, returns global app index.
static InputParameters validParams()
static void addBBoxFactorParam(InputParameters ¶ms)
Add the bounding box factor parameter to the supplied input parameters.
const InputParameters & parameters() const
Get the parameters of the object.
void projectSolution(unsigned int to_problem)
virtual void set(const numeric_index_type i, const Number value)=0
std::vector< LinearImplicitSystem * > _proj_sys
virtual void initialSetup() override
Method called at the beginning of the simulation for checking integrity or doing one-time setup...
unsigned int _proj_var_num
Having one projection variable number seems weird, but there is always one variable in every system b...
processor_id_type processor_id() const
SystemBase & sys()
Get the system this variable is part of.
friend void assemble_l2(EquationSystems &es, const std::string &system_name)
std::vector< FEProblemBase * > _from_problems
std::vector< MooseMesh * > _to_meshes
VariableName _from_var_name
This values are used if a derived class only supports one variable.