Skip to content

Commit

Permalink
Add weight functions (#21960)
Browse files Browse the repository at this point in the history
  • Loading branch information
dschwen committed Aug 29, 2022
1 parent 712a86c commit d5e6a7b
Show file tree
Hide file tree
Showing 12 changed files with 81 additions and 28 deletions.
10 changes: 8 additions & 2 deletions framework/doc/content/source/userobject/RadialAverage.md
Expand Up @@ -6,9 +6,10 @@ Given a material property and a radius for averaging, the RadialAverage object
computes the spatial average value of the property over the radius.

## Applications

This can be used for nonlocal damage models in the `TensorMechanics` module
where the damage_index that is used for computing the damage stress is average
over a certain radius. This can help alleviate mesh sensitivity in certain
over a certain radius $r_0$. This can help alleviate mesh sensitivity in certain
cases. This can be accomplished by running the RadialAverage object on a local
damage material property. Then using the `NonlocalDamage` model in conjunction
with the `ComputeDamageStress` the damage index used for updating the stress is
Expand All @@ -20,7 +21,7 @@ The RadialAverage user object is derived from `ElementUserObject` and
works in two stages.

1. In the element loop (in the `execute()` method) a list of all quadrature
points, their locations, indices, and selected variable value is compiled.
points, their locations, indices, and selected material property value is compiled.

2. In the `finalize()` method

Expand All @@ -36,6 +37,11 @@ works in two stages.
4. the results from the search are used to spatially averaged


The [!param](/UserObjects/RadialAverage/weights) parameter determines the distance
based weight function to be used in the averaging process. `constant` assigns an equal weight
to all material points, `linear` weights each point by $r_0-r$ (a linear fall-off with distance),
and `cosine` weights each point by $1+\cos (\frac r{r_0}\pi}$.

!syntax parameters /UserObjects/RadialAverage

!syntax inputs /UserObjects/RadialAverage
Expand Down
12 changes: 9 additions & 3 deletions framework/include/userobject/RadialAverage.h
Expand Up @@ -72,10 +72,16 @@ class RadialAverage : public ElementUserObject
protected:
void updateCommunicationLists();

/// material name to get averaged
std::string _averaged_material_name;
/// distance based weight function
enum class WeightsType
{
CONSTANT,
LINEAR,
COSINE
} _weights_type;

/// material to be averaged
const MaterialProperty<Real> & _averaged_material;
const MaterialProperty<Real> & _prop;

/// cut-off radius
const Real _radius;
Expand Down
13 changes: 8 additions & 5 deletions framework/src/userobject/RadialAverage.C
Expand Up @@ -34,8 +34,11 @@ InputParameters
RadialAverage::validParams()
{
InputParameters params = ElementUserObject::validParams();
params.addClassDescription("Perform a radial equal weight average of a material property");
params.addRequiredParam<std::string>("material_name", "Name of the material to average.");
params.addClassDescription("Perform a radial average of a material property");
params.addRequiredParam<MaterialPropertyName>("prop_name",
"Name of the material property to average");
MooseEnum weights_type("constant linear cosine", "linear");
params.addRequiredParam<MooseEnum>("weights", weights_type, "Distance based weight function");
params.addRequiredParam<Real>("radius", "Cut-off radius for the averaging");
params.addRangeCheckedParam<Real>(
"padding",
Expand All @@ -57,8 +60,8 @@ RadialAverage::validParams()

RadialAverage::RadialAverage(const InputParameters & parameters)
: ElementUserObject(parameters),
_averaged_material_name(getParam<std::string>("material_name")),
_averaged_material(getMaterialProperty<Real>(_averaged_material_name)),
_weights_type(getParam<MooseEnum>("weights").getEnum<WeightsType>()),
_prop(getMaterialProperty<Real>("prop_name")),
_radius(getParam<Real>("radius")),
_padding(getParam<Real>("padding")),
_update_communication_lists(false),
Expand Down Expand Up @@ -89,7 +92,7 @@ RadialAverage::execute()

// collect all QP data
for (const auto qp : make_range(_qrule->n_points()))
_qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _averaged_material[qp]);
_qp_data.emplace_back(_q_point[qp], id, qp, _JxW[qp] * _coord[qp], _prop[qp]);

// make sure the result map entry for the current element is sized correctly
auto i = _average.find(id);
Expand Down
20 changes: 18 additions & 2 deletions framework/src/userobject/ThreadedRadialAverageLoop.C
Expand Up @@ -26,6 +26,7 @@ ThreadedRadialAverageLoop::operator()(const QPDataRange & qpdata_range)
const auto radius = _radavg._radius;
const auto & qp_data = _radavg._qp_data;
const auto & kd_tree = _radavg._kd_tree;
const auto & weights_type = _radavg._weights_type;

// tree search data structures
std::vector<std::pair<std::size_t, Real>> ret_matches;
Expand All @@ -52,11 +53,26 @@ ThreadedRadialAverageLoop::operator()(const QPDataRange & qpdata_range)
std::size_t n_result =
kd_tree->radiusSearch(&(local_qp._q_point(0)), radius * radius, ret_matches, search_params);
Real total_vol = 0.0;
Real weight = 1.0;
for (std::size_t j = 0; j < n_result; ++j)
{
const auto & other_qp = qp_data[ret_matches[j].first];
sum += other_qp._value * other_qp._volume;
total_vol += other_qp._volume;
switch (weights_type)
{
case RadialAverage::WeightsType::CONSTANT:
break;

case RadialAverage::WeightsType::LINEAR:
weight = radius - std::sqrt(ret_matches[j].second);
break;

case RadialAverage::WeightsType::COSINE:
weight = std::cos(std::sqrt(ret_matches[j].second) / radius * libMesh::pi) + 1.0;
break;
}

sum += other_qp._value * other_qp._volume * weight;
total_vol += other_qp._volume * weight;
}
sum /= total_vol;
}
Expand Down
Expand Up @@ -60,9 +60,8 @@
[UserObjects]
[ele_avg]
type = RadialAverage
material_name = local_damage_reg
prop_name = local_damage_reg
execute_on = "INITIAL timestep_end"
block = 0
radius = 0.55
[]
[]
Expand Down
Expand Up @@ -59,9 +59,8 @@
[UserObjects]
[ele_avg]
type = RadialAverage
material_name = local_damage
prop_name = local_damage
execute_on = "INITIAL timestep_end"
block = 0
radius = 0.55
[]
[]
Expand Down
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Expand Up @@ -18,13 +18,14 @@
[AuxVariables]
[non_local_material]
family = MONOMIAL
order = SECOND
[]
[]

[AuxKernels]
[non_local]
type = RadialAverageAux
average_UO = ele_avg
average_UO = average
variable = non_local_material
execute_on = 'INITIAL TIMESTEP_BEGIN'
[]
Expand All @@ -46,11 +47,11 @@
[]

[UserObjects]
[ele_avg]
[average]
type = RadialAverage
material_name = local
prop_name = local
weights = constant
execute_on = "INITIAL timestep_end"
block = 0
radius = 0.3
[]
[]
Expand Down
39 changes: 31 additions & 8 deletions test/tests/userobjects/radial_average/tests
@@ -1,13 +1,36 @@
[Tests]
issues = '#21786'
issues = '#21786 #21944 #21963'
design = '/RadialAverage.md'

[radial_average]
type = 'Exodiff'
input = 'time_changing_test.i'
exodiff = 'time_changing_test_out.e'
recover = false
requirement = 'The system shall will compute a radial average of a material that changes over time.'
abs_zero = 1e-10
[weight]
requirement = 'The system shall will compute a radial average of a material property that changes over time with '
[constant]
type = 'Exodiff'
input = 'test.i'
cli_args = 'Outputs/file_base=constant UserObjects/average/weights=constant'
exodiff = 'constant.e'
recover = false
abs_zero = 1e-10
detail = 'an equal weight for all material points'
[]
[linear]
type = 'Exodiff'
input = 'test.i'
cli_args = 'Outputs/file_base=linear UserObjects/average/weights=linear'
exodiff = 'linear.e'
recover = false
abs_zero = 1e-10
detail = 'a weight factor that falls off linearly to zero with distance'
[]
[cosine]
type = 'Exodiff'
input = 'test.i'
cli_args = 'Outputs/file_base=cosine UserObjects/average/weights=cosine'
exodiff = 'cosine.e'
recover = false
abs_zero = 1e-10
detail = 'a weight factor that falls off with a cosine profile with distance'
[]
[]
[]

0 comments on commit d5e6a7b

Please sign in to comment.