forked from geodynamics/aspect
/
stress.cc
115 lines (95 loc) · 4.65 KB
/
stress.cc
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
/*
Copyright (C) 2011 - 2020 by the authors of the ASPECT code.
This file is part of ASPECT.
ASPECT is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
ASPECT is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ASPECT; see the file LICENSE. If not see
<http://www.gnu.org/licenses/>.
*/
#include <aspect/postprocess/visualization/stress.h>
namespace aspect
{
namespace Postprocess
{
namespace VisualizationPostprocessors
{
template <int dim>
Stress<dim>::
Stress ()
:
DataPostprocessorTensor<dim> ("stress",
update_values | update_gradients | update_quadrature_points)
{}
template <int dim>
void
Stress<dim>::
evaluate_vector_field(const DataPostprocessorInputs::Vector<dim> &input_data,
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = input_data.solution_values.size();
Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
Assert ((computed_quantities[0].size() == Tensor<2,dim>::n_independent_components),
ExcInternalError());
Assert (input_data.solution_values[0].size() == this->introspection().n_components, ExcInternalError());
Assert (input_data.solution_gradients[0].size() == this->introspection().n_components, ExcInternalError());
MaterialModel::MaterialModelInputs<dim> in(input_data,
this->introspection());
MaterialModel::MaterialModelOutputs<dim> out(n_quadrature_points,
this->n_compositional_fields());
// Compute the viscosity...
this->get_material_model().evaluate(in, out);
// ...and use it to compute the stresses
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
const SymmetricTensor<2,dim> compressible_strain_rate
= (this->get_material_model().is_compressible()
?
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);
const double eta = out.viscosities[q];
// Compressive stress is positive in geoscience applications
const SymmetricTensor<2,dim> stress = -2.*eta*compressible_strain_rate +
in.pressure[q] * unit_symmetric_tensor<dim>();
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
computed_quantities[q][Tensor<2,dim>::component_to_unrolled_index(TableIndices<2>(d,e))]
= stress[d][e];
}
// average the values if requested
const auto &viz = this->get_postprocess_manager().template get_matching_postprocessor<Postprocess::Visualization<dim> >();
if (!viz.output_pointwise_stress_and_strain())
average_quantities(computed_quantities);
}
}
}
}
// explicit instantiations
namespace aspect
{
namespace Postprocess
{
namespace VisualizationPostprocessors
{
ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(Stress,
"stress",
"A visualization output object that generates output "
"for the 3 (in 2d) or 6 (in 3d) components of the stress "
"tensor, i.e., for the components of the tensor "
"$-2\\eta\\varepsilon(\\mathbf u)+pI$ "
"in the incompressible case and "
"$-2\\eta\\left[\\varepsilon(\\mathbf u)-"
"\\tfrac 13(\\textrm{tr}\\;\\varepsilon(\\mathbf u))\\mathbf I\\right]+pI$ "
"in the compressible case. Note that the convention of positive "
"compressive stress is followed. ")
}
}
}