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 "FeatureVolumeVectorPostprocessor.h"
11 :
12 : // MOOSE includes
13 : #include "Assembly.h"
14 : #include "FeatureFloodCount.h"
15 : #include "GrainTrackerInterface.h"
16 : #include "MooseMesh.h"
17 : #include "MooseVariable.h"
18 : #include "SystemBase.h"
19 :
20 : #include "libmesh/quadrature.h"
21 :
22 : registerMooseObject("PhaseFieldApp", FeatureVolumeVectorPostprocessor);
23 :
24 : InputParameters
25 836 : FeatureVolumeVectorPostprocessor::validParams()
26 : {
27 836 : InputParameters params = GeneralVectorPostprocessor::validParams();
28 836 : params += BoundaryRestrictable::validParams();
29 :
30 1672 : params.addRequiredParam<UserObjectName>("flood_counter",
31 : "The FeatureFloodCount UserObject to get values from.");
32 1672 : params.addParam<bool>("single_feature_per_element",
33 1672 : false,
34 : "Set this Boolean if you wish to use an element based volume where"
35 : " the dominant order parameter determines the feature that accumulates the "
36 : "entire element volume");
37 1672 : params.addParam<bool>("output_centroids", false, "Set to true to output the feature centroids");
38 836 : params.addClassDescription("This object is designed to pull information from the data structures "
39 : "of a \"FeatureFloodCount\" or derived object (e.g. individual "
40 : "feature volumes)");
41 :
42 836 : params.suppressParameter<bool>("contains_complete_history");
43 :
44 836 : return params;
45 0 : }
46 :
47 418 : FeatureVolumeVectorPostprocessor::FeatureVolumeVectorPostprocessor(
48 418 : const InputParameters & parameters)
49 : : GeneralVectorPostprocessor(parameters),
50 : MooseVariableDependencyInterface(this),
51 : BoundaryRestrictable(this, false),
52 418 : _single_feature_per_elem(getParam<bool>("single_feature_per_element")),
53 836 : _output_centroids(getParam<bool>("output_centroids")),
54 418 : _feature_counter(getUserObject<FeatureFloodCount>("flood_counter")),
55 418 : _var_num(declareVector("var_num")),
56 418 : _feature_volumes(declareVector("feature_volumes")),
57 418 : _intersects_bounds(declareVector("intersects_bounds")),
58 418 : _intersects_specified_bounds(declareVector("intersects_specified_bounds")),
59 418 : _percolated(declareVector("percolated")),
60 418 : _vars(_feature_counter.getFECoupledVars()),
61 418 : _mesh(_subproblem.mesh()),
62 418 : _assembly(_subproblem.assembly(_tid, _sys.number())),
63 418 : _q_point(_assembly.qPoints()),
64 418 : _qrule(_assembly.qRule()),
65 418 : _JxW(_assembly.JxW()),
66 418 : _coord(_assembly.coordTransformation()),
67 418 : _qrule_face(_assembly.qRuleFace()),
68 418 : _JxW_face(_assembly.JxWFace())
69 : {
70 418 : addMooseVariableDependency(_vars);
71 :
72 418 : _is_boundary_restricted = boundaryRestricted();
73 :
74 418 : _coupled_sln.reserve(_vars.size());
75 1067 : for (auto & var : _feature_counter.getCoupledVars())
76 649 : _coupled_sln.push_back(&var->sln());
77 :
78 418 : const std::array<std::string, 3> suffix = {{"x", "y", "z"}};
79 418 : if (_output_centroids)
80 516 : for (unsigned int i = 0; i < 3; ++i)
81 387 : _centroid[i] = &declareVector("centroid_" + suffix[i]);
82 418 : }
83 :
84 : void
85 454 : FeatureVolumeVectorPostprocessor::initialize()
86 : {
87 454 : }
88 :
89 : void
90 454 : FeatureVolumeVectorPostprocessor::execute()
91 : {
92 454 : const auto num_features = _feature_counter.getTotalFeatureCount();
93 :
94 : // Reset the variable index and intersect bounds vectors
95 454 : _var_num.assign(num_features, -1); // Invalid
96 454 : _intersects_bounds.assign(num_features, -1); // Invalid
97 454 : _intersects_specified_bounds.assign(num_features, -1); // Invalid
98 454 : _percolated.assign(num_features, -1); // Invalid
99 1557 : for (MooseIndex(num_features) feature_num = 0; feature_num < num_features; ++feature_num)
100 : {
101 1103 : auto var_num = _feature_counter.getFeatureVar(feature_num);
102 1103 : if (var_num != FeatureFloodCount::invalid_id)
103 941 : _var_num[feature_num] = var_num;
104 :
105 1103 : _intersects_bounds[feature_num] =
106 1103 : static_cast<unsigned int>(_feature_counter.doesFeatureIntersectBoundary(feature_num));
107 :
108 1103 : _intersects_specified_bounds[feature_num] = static_cast<unsigned int>(
109 1103 : _feature_counter.doesFeatureIntersectSpecifiedBoundary(feature_num));
110 :
111 1103 : _percolated[feature_num] =
112 1103 : static_cast<unsigned int>(_feature_counter.isFeaturePercolated(feature_num));
113 : }
114 :
115 454 : if (_output_centroids)
116 : {
117 484 : for (std::size_t i = 0; i < 3; ++i)
118 363 : _centroid[i]->resize(num_features);
119 558 : for (std::size_t feature_num = 0; feature_num < num_features; ++feature_num)
120 : {
121 437 : auto p = _feature_counter.featureCentroid(feature_num);
122 1748 : for (std::size_t i = 0; i < 3; ++i)
123 1311 : (*_centroid[i])[feature_num] = p(i);
124 : }
125 : }
126 :
127 : // Reset the volume vector
128 454 : _feature_volumes.assign(num_features, 0);
129 :
130 : // Calculate coverage of a boundary if one has been supplied in the input file
131 454 : if (_is_boundary_restricted)
132 : {
133 48 : const std::set<BoundaryID> supplied_bnd_ids = BoundaryRestrictable::boundaryIDs();
134 44616 : for (auto elem_it = _mesh.bndElemsBegin(), elem_end = _mesh.bndElemsEnd(); elem_it != elem_end;
135 44520 : ++elem_it)
136 :
137 : // loop over only boundaries supplied by user in boundary param
138 89040 : for (auto & supplied_bnd_id : supplied_bnd_ids)
139 44520 : if (((*elem_it)->_bnd_id) == supplied_bnd_id)
140 : {
141 10200 : const auto & elem = (*elem_it)->_elem;
142 : auto rank = processor_id();
143 :
144 10200 : if (elem->processor_id() == rank)
145 : {
146 5100 : _fe_problem.setCurrentSubdomainID(elem, 0);
147 5100 : _fe_problem.prepare(elem, 0);
148 5100 : _fe_problem.reinitElem(elem, 0);
149 5100 : _fe_problem.reinitElemFace(elem, (*elem_it)->_side, 0);
150 :
151 5100 : const auto & var_to_features = _feature_counter.getVarToFeatureVector(elem->id());
152 :
153 5100 : accumulateBoundaryFaces(elem, var_to_features, num_features, (*elem_it)->_side);
154 : }
155 : }
156 : }
157 : else // If no boundary is supplied, calculate volumes of features as normal
158 1080530 : for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
159 : {
160 539859 : _fe_problem.setCurrentSubdomainID(elem, 0);
161 539859 : _fe_problem.prepare(elem, 0);
162 539859 : _fe_problem.reinitElem(elem, 0);
163 :
164 : /**
165 : * Here we retrieve the var to features vector on the current element.
166 : * We'll use that information to figure out which variables are non-zero
167 : * (from a threshold perspective) then we can sum those values into
168 : * appropriate grain index locations.
169 : */
170 539859 : const auto & var_to_features = _feature_counter.getVarToFeatureVector(elem->id());
171 :
172 539859 : accumulateVolumes(elem, var_to_features, num_features);
173 406 : }
174 454 : }
175 :
176 : void
177 454 : FeatureVolumeVectorPostprocessor::finalize()
178 : {
179 : // Do the parallel sum
180 454 : _communicator.sum(_feature_volumes);
181 454 : }
182 :
183 : Real
184 0 : FeatureVolumeVectorPostprocessor::getFeatureVolume(unsigned int feature_id) const
185 : {
186 : mooseAssert(feature_id < _feature_volumes.size(), "feature_id is out of range");
187 0 : return _feature_volumes[feature_id];
188 : }
189 :
190 : void
191 539859 : FeatureVolumeVectorPostprocessor::accumulateVolumes(
192 : const Elem * elem,
193 : const std::vector<unsigned int> & var_to_features,
194 : std::size_t libmesh_dbg_var(num_features))
195 : {
196 539859 : unsigned int dominant_feature_id = FeatureFloodCount::invalid_id;
197 : Real max_var_value = std::numeric_limits<Real>::lowest();
198 :
199 1246303 : for (MooseIndex(var_to_features) var_index = 0; var_index < var_to_features.size(); ++var_index)
200 : {
201 : // Only sample "active" variables
202 706444 : if (var_to_features[var_index] != FeatureFloodCount::invalid_id)
203 : {
204 : auto feature_id = var_to_features[var_index];
205 : mooseAssert(feature_id < num_features, "Feature ID out of range");
206 245034 : auto integral_value = computeIntegral(var_index);
207 :
208 : // Compute volumes in a simplistic but domain conservative fashion
209 245034 : if (_single_feature_per_elem)
210 : {
211 24624 : if (integral_value > max_var_value)
212 : {
213 : // Update the current dominant feature and associated value
214 : max_var_value = integral_value;
215 : dominant_feature_id = feature_id;
216 : }
217 : }
218 : // Solution based volume calculation (integral value)
219 : else
220 220410 : _feature_volumes[feature_id] += integral_value;
221 : }
222 : }
223 :
224 : // Accumulate the entire element volume into the dominant feature. Do not use the integral value
225 539859 : if (_single_feature_per_elem && dominant_feature_id != FeatureFloodCount::invalid_id)
226 24168 : _feature_volumes[dominant_feature_id] += _assembly.elementVolume(elem);
227 539859 : }
228 :
229 : Real
230 245034 : FeatureVolumeVectorPostprocessor::computeIntegral(std::size_t var_index) const
231 : {
232 : Real sum = 0;
233 :
234 974114 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
235 729080 : sum += _JxW[qp] * _coord[qp] * (*_coupled_sln[var_index])[qp];
236 :
237 245034 : return sum;
238 : }
239 :
240 : void
241 5100 : FeatureVolumeVectorPostprocessor::accumulateBoundaryFaces(
242 : const Elem * elem,
243 : const std::vector<unsigned int> & var_to_features,
244 : std::size_t libmesh_dbg_var(num_features),
245 : unsigned int side)
246 : {
247 5100 : unsigned int dominant_feature_id = FeatureFloodCount::invalid_id;
248 : Real max_var_value = std::numeric_limits<Real>::lowest();
249 :
250 10200 : for (MooseIndex(var_to_features) var_index = 0; var_index < var_to_features.size(); ++var_index)
251 : {
252 : // Only sample "active" variables
253 5100 : if (var_to_features[var_index] != FeatureFloodCount::invalid_id)
254 : {
255 : auto feature_id = var_to_features[var_index];
256 : mooseAssert(feature_id < num_features, "Feature ID out of range");
257 1536 : auto integral_value = computeFaceIntegral(var_index);
258 :
259 1536 : if (_single_feature_per_elem)
260 : {
261 558 : if (integral_value > max_var_value)
262 : {
263 : // Update the current dominant feature and associated value
264 : max_var_value = integral_value;
265 : dominant_feature_id = feature_id;
266 : }
267 : }
268 : // Solution based boundary area/length calculation (integral value)
269 : else
270 978 : _feature_volumes[feature_id] += integral_value;
271 : }
272 : }
273 :
274 : // Accumulate the boundary area/length into the dominant feature. Do not use the integral value
275 5100 : if (_single_feature_per_elem && dominant_feature_id != FeatureFloodCount::invalid_id)
276 : {
277 558 : std::unique_ptr<const Elem> side_elem = elem->build_side_ptr(side);
278 558 : _feature_volumes[dominant_feature_id] += _assembly.elementVolume(side_elem.get());
279 558 : }
280 5100 : }
281 :
282 : Real
283 1536 : FeatureVolumeVectorPostprocessor::computeFaceIntegral(std::size_t var_index) const
284 : {
285 : Real sum = 0;
286 6948 : for (unsigned int qp = 0; qp < _qrule_face->n_points(); ++qp)
287 5412 : sum += _JxW_face[qp] * _coord[qp] * (*_coupled_sln[var_index])[qp];
288 :
289 1536 : return sum;
290 : }
|
