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 592 : FeatureVolumeVectorPostprocessor::validParams()
26 : {
27 592 : InputParameters params = GeneralVectorPostprocessor::validParams();
28 592 : params += BoundaryRestrictable::validParams();
29 :
30 1184 : params.addRequiredParam<UserObjectName>("flood_counter",
31 : "The FeatureFloodCount UserObject to get values from.");
32 1184 : params.addParam<bool>("single_feature_per_element",
33 1184 : 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 1184 : params.addParam<bool>("output_centroids", false, "Set to true to output the feature centroids");
38 592 : 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 592 : params.suppressParameter<bool>("contains_complete_history");
43 :
44 592 : return params;
45 0 : }
46 :
47 296 : FeatureVolumeVectorPostprocessor::FeatureVolumeVectorPostprocessor(
48 296 : const InputParameters & parameters)
49 : : GeneralVectorPostprocessor(parameters),
50 : MooseVariableDependencyInterface(this),
51 : BoundaryRestrictable(this, false),
52 296 : _single_feature_per_elem(getParam<bool>("single_feature_per_element")),
53 592 : _output_centroids(getParam<bool>("output_centroids")),
54 296 : _feature_counter(getUserObject<FeatureFloodCount>("flood_counter")),
55 296 : _var_num(declareVector("var_num")),
56 296 : _feature_volumes(declareVector("feature_volumes")),
57 296 : _intersects_bounds(declareVector("intersects_bounds")),
58 296 : _intersects_specified_bounds(declareVector("intersects_specified_bounds")),
59 296 : _percolated(declareVector("percolated")),
60 296 : _vars(_feature_counter.getFECoupledVars()),
61 296 : _mesh(_subproblem.mesh()),
62 296 : _assembly(_subproblem.assembly(_tid, _sys.number())),
63 296 : _q_point(_assembly.qPoints()),
64 296 : _qrule(_assembly.qRule()),
65 296 : _JxW(_assembly.JxW()),
66 296 : _coord(_assembly.coordTransformation()),
67 296 : _qrule_face(_assembly.qRuleFace()),
68 296 : _JxW_face(_assembly.JxWFace())
69 : {
70 296 : addMooseVariableDependency(_vars);
71 :
72 296 : _is_boundary_restricted = boundaryRestricted();
73 :
74 296 : _coupled_sln.reserve(_vars.size());
75 759 : for (auto & var : _feature_counter.getCoupledVars())
76 463 : _coupled_sln.push_back(&var->sln());
77 :
78 296 : const std::array<std::string, 3> suffix = {{"x", "y", "z"}};
79 296 : if (_output_centroids)
80 372 : for (unsigned int i = 0; i < 3; ++i)
81 279 : _centroid[i] = &declareVector("centroid_" + suffix[i]);
82 296 : }
83 :
84 : void
85 358 : FeatureVolumeVectorPostprocessor::initialize()
86 : {
87 358 : }
88 :
89 : void
90 358 : FeatureVolumeVectorPostprocessor::execute()
91 : {
92 358 : const auto num_features = _feature_counter.getTotalFeatureCount();
93 :
94 : // Reset the variable index and intersect bounds vectors
95 358 : _var_num.assign(num_features, -1); // Invalid
96 358 : _intersects_bounds.assign(num_features, -1); // Invalid
97 358 : _intersects_specified_bounds.assign(num_features, -1); // Invalid
98 358 : _percolated.assign(num_features, -1); // Invalid
99 1241 : for (MooseIndex(num_features) feature_num = 0; feature_num < num_features; ++feature_num)
100 : {
101 883 : auto var_num = _feature_counter.getFeatureVar(feature_num);
102 883 : if (var_num != FeatureFloodCount::invalid_id)
103 752 : _var_num[feature_num] = var_num;
104 :
105 883 : _intersects_bounds[feature_num] =
106 883 : static_cast<unsigned int>(_feature_counter.doesFeatureIntersectBoundary(feature_num));
107 :
108 883 : _intersects_specified_bounds[feature_num] = static_cast<unsigned int>(
109 883 : _feature_counter.doesFeatureIntersectSpecifiedBoundary(feature_num));
110 :
111 883 : _percolated[feature_num] =
112 883 : static_cast<unsigned int>(_feature_counter.isFeaturePercolated(feature_num));
113 : }
114 :
115 358 : if (_output_centroids)
116 : {
117 388 : for (std::size_t i = 0; i < 3; ++i)
118 291 : _centroid[i]->resize(num_features);
119 450 : for (std::size_t feature_num = 0; feature_num < num_features; ++feature_num)
120 : {
121 353 : auto p = _feature_counter.featureCentroid(feature_num);
122 1412 : for (std::size_t i = 0; i < 3; ++i)
123 1059 : (*_centroid[i])[feature_num] = p(i);
124 : }
125 : }
126 :
127 : // Reset the volume vector
128 358 : _feature_volumes.assign(num_features, 0);
129 :
130 : // Calculate coverage of a boundary if one has been supplied in the input file
131 358 : if (_is_boundary_restricted)
132 : {
133 40 : const std::set<BoundaryID> supplied_bnd_ids = BoundaryRestrictable::boundaryIDs();
134 37180 : for (auto elem_it = _mesh.bndElemsBegin(), elem_end = _mesh.bndElemsEnd(); elem_it != elem_end;
135 37100 : ++elem_it)
136 :
137 : // loop over only boundaries supplied by user in boundary param
138 74200 : for (auto & supplied_bnd_id : supplied_bnd_ids)
139 37100 : if (((*elem_it)->_bnd_id) == supplied_bnd_id)
140 : {
141 8500 : const auto & elem = (*elem_it)->_elem;
142 : auto rank = processor_id();
143 :
144 8500 : if (elem->processor_id() == rank)
145 : {
146 4250 : _fe_problem.setCurrentSubdomainID(elem, 0);
147 4250 : _fe_problem.prepare(elem, 0);
148 4250 : _fe_problem.reinitElem(elem, 0);
149 4250 : _fe_problem.reinitElemFace(elem, (*elem_it)->_side, 0);
150 :
151 4250 : const auto & var_to_features = _feature_counter.getVarToFeatureVector(elem->id());
152 :
153 4250 : 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 901538 : for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
159 : {
160 450451 : _fe_problem.setCurrentSubdomainID(elem, 0);
161 450451 : _fe_problem.prepare(elem, 0);
162 450451 : _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 450451 : const auto & var_to_features = _feature_counter.getVarToFeatureVector(elem->id());
171 :
172 450451 : accumulateVolumes(elem, var_to_features, num_features);
173 318 : }
174 358 : }
175 :
176 : void
177 358 : FeatureVolumeVectorPostprocessor::finalize()
178 : {
179 : // Do the parallel sum
180 358 : _communicator.sum(_feature_volumes);
181 358 : }
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 450451 : 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 450451 : unsigned int dominant_feature_id = FeatureFloodCount::invalid_id;
197 : Real max_var_value = std::numeric_limits<Real>::lowest();
198 :
199 1039927 : for (MooseIndex(var_to_features) var_index = 0; var_index < var_to_features.size(); ++var_index)
200 : {
201 : // Only sample "active" variables
202 589476 : 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 204431 : auto integral_value = computeIntegral(var_index);
207 :
208 : // Compute volumes in a simplistic but domain conservative fashion
209 204431 : if (_single_feature_per_elem)
210 : {
211 20520 : 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 183911 : _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 450451 : if (_single_feature_per_elem && dominant_feature_id != FeatureFloodCount::invalid_id)
226 20140 : _feature_volumes[dominant_feature_id] += _assembly.elementVolume(elem);
227 450451 : }
228 :
229 : Real
230 204431 : FeatureVolumeVectorPostprocessor::computeIntegral(std::size_t var_index) const
231 : {
232 : Real sum = 0;
233 :
234 813075 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
235 608644 : sum += _JxW[qp] * _coord[qp] * (*_coupled_sln[var_index])[qp];
236 :
237 204431 : return sum;
238 : }
239 :
240 : void
241 4250 : 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 4250 : unsigned int dominant_feature_id = FeatureFloodCount::invalid_id;
248 : Real max_var_value = std::numeric_limits<Real>::lowest();
249 :
250 8500 : for (MooseIndex(var_to_features) var_index = 0; var_index < var_to_features.size(); ++var_index)
251 : {
252 : // Only sample "active" variables
253 4250 : 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 1280 : auto integral_value = computeFaceIntegral(var_index);
258 :
259 1280 : if (_single_feature_per_elem)
260 : {
261 465 : 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 815 : _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 4250 : if (_single_feature_per_elem && dominant_feature_id != FeatureFloodCount::invalid_id)
276 : {
277 465 : std::unique_ptr<const Elem> side_elem = elem->build_side_ptr(side);
278 465 : _feature_volumes[dominant_feature_id] += _assembly.elementVolume(side_elem.get());
279 465 : }
280 4250 : }
281 :
282 : Real
283 1280 : FeatureVolumeVectorPostprocessor::computeFaceIntegral(std::size_t var_index) const
284 : {
285 : Real sum = 0;
286 5790 : for (unsigned int qp = 0; qp < _qrule_face->n_points(); ++qp)
287 4510 : sum += _JxW_face[qp] * _coord[qp] * (*_coupled_sln[var_index])[qp];
288 :
289 1280 : return sum;
290 : }
|
