forked from idaholab/moose
/
IntegralReportingIDVPP.C
165 lines (150 loc) · 5.72 KB
/
IntegralReportingIDVPP.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
#include "IntegralReportingIDVPP.h"
#include "MooseVariable.h"
#include "MooseMesh.h"
#include "libmesh/mesh_base.h"
registerMooseObject("ReactorApp", IntegralReportingIDVPP);
InputParameters
IntegralReportingIDVPP::validParams()
{
InputParameters params = ElementVariableVectorPostprocessor::validParams();
params.addRequiredParam<std::vector<ExtraElementIDName>>(
"reporting_id", "The types of reporting id by which to separate integrals.");
params.addClassDescription(
"This IntegralReportingIDVPP source code is to integrate variables based on reporting IDs");
return params;
}
IntegralReportingIDVPP::IntegralReportingIDVPP(const InputParameters & parameters)
: ElementVariableVectorPostprocessor(parameters),
_nvar(coupledComponents("variable")),
_reporting_id(getParam<std::vector<ExtraElementIDName>>("reporting_id")),
_n_reporting_id(_reporting_id.size())
{
// numbers of parsed vpp ids
dof_id_type n_vpp_ids;
// list of extra_integer id
std::set<dof_id_type> ids;
std::set<ExtraElementIDName> _reporting_id_copy(_reporting_id.begin(), _reporting_id.end());
// initial setting for parsing
ExtraElementIDName id_type = *_reporting_id_copy.begin();
_reporting_id_copy.erase(id_type);
unsigned int id_index = getElementIDIndexByName(id_type);
if (blockRestricted())
ids = getElemIDsOnBlocks(id_index, blockIDs());
else
ids = getAllElemIDs(id_index);
n_vpp_ids = ids.size();
for (auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
// initializing vpp_id based on the first reporting id value
auto id = elem->get_extra_integer(id_index);
_unique_vpp_ids[elem->id()] = std::distance(ids.begin(), std::find(ids.begin(), ids.end(), id));
}
// looping over reporting ids
while (!_reporting_id_copy.empty())
{
id_type = *_reporting_id_copy.begin();
_reporting_id_copy.erase(id_type);
id_index = getElementIDIndexByName(id_type);
// parsing reporting ids (find unique combinations of reporting ids)
std::set<std::pair<dof_id_type, dof_id_type>> unique_ids;
for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
auto elem_id = elem->id();
dof_id_type id1 = _unique_vpp_ids[elem_id];
dof_id_type id2 = elem->get_extra_integer(id_index);
auto it = unique_ids.find(std::pair<dof_id_type, dof_id_type>(id1, id2));
if (it == unique_ids.end())
unique_ids.insert(std::pair<dof_id_type, dof_id_type>(id1, id2));
}
comm().set_union(unique_ids);
// update parsed vpp ids
n_vpp_ids = unique_ids.size();
for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
dof_id_type id1 = _unique_vpp_ids[elem->id()];
dof_id_type id2 = elem->get_extra_integer(id_index);
_unique_vpp_ids[elem->id()] = std::distance(
unique_ids.begin(),
std::find(
unique_ids.begin(), unique_ids.end(), std::pair<dof_id_type, dof_id_type>(id1, id2)));
}
}
// set up variable vector containing parsed reporting ids
for (unsigned int i = 0; i < _n_reporting_id; ++i)
{
id_type = _reporting_id[i];
auto & p = declareVector("Level-" + std::to_string(i) + "-" + id_type);
p.resize(n_vpp_ids);
// collect local map of local vpp id to reporting id
std::map<dof_id_type, dof_id_type> reporting_ids;
for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
{
if (!hasBlocks(elem->subdomain_id()))
continue;
auto vpp_id = _unique_vpp_ids[elem->id()];
if (reporting_ids.find(vpp_id) == reporting_ids.end())
reporting_ids[vpp_id] = elem->get_extra_integer(getElementIDIndexByName(id_type));
}
// create global map of local vpp id to reporting id
comm().set_union(reporting_ids);
for (auto it : reporting_ids)
p[it.first] = it.second;
_var_reporting_ids.push_back(&p);
}
// declare vectors containg integral values
for (unsigned int i = 0; i < _nvar; ++i)
{
_vars.push_back(getVar("variable", i));
_var_values.push_back(&coupledValue("variable", i));
auto & p = declareVector(_vars[i]->name());
p.resize(n_vpp_ids);
_var_integrals.push_back(&p);
}
}
void
IntegralReportingIDVPP::initialize()
{
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
for (size_t i = 0; i < (*_var_integrals[ivar]).size(); ++i)
(*_var_integrals[ivar])[i] = 0;
}
void
IntegralReportingIDVPP::execute()
{
if (hasBlocks(_current_elem->subdomain_id()))
{
auto ipos = _unique_vpp_ids[_current_elem->id()];
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
if (_vars[ivar]->hasBlocks(_current_elem->subdomain_id()))
for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
(*_var_integrals[ivar])[ipos] += _JxW[qp] * _coord[qp] * (*_var_values[ivar])[qp];
}
}
void
IntegralReportingIDVPP::finalize()
{
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
gatherSum(*_var_integrals[ivar]);
}
void
IntegralReportingIDVPP::threadJoin(const UserObject & s)
{
const IntegralReportingIDVPP & sibling = static_cast<const IntegralReportingIDVPP &>(s);
for (unsigned int ivar = 0; ivar < _nvar; ++ivar)
for (size_t i = 0; i < (*_var_integrals[ivar]).size(); ++i)
(*_var_integrals[ivar])[i] += (*sibling._var_integrals[ivar])[i];
}