Hello, Again.
I am trying to post-process an elastic problem for which I have cells with
different material properties. I would like to determine von-mises stress
and output it to the solution file. Presently I save the cell materials
properties in the assemble() function to arrays. I was going to zero the
array index before call to add_data_vector,
ComputeVonMises<dim> vonmises;
....
iCELL = 0;
data_out.add_data_vector (solution, vonmises);
with the class ComputeVonMises as defined below. However it appears that
the function is threaded (nice) but the incrementing iCELL method will not
work. The output of the two iCELL fprintf() statments in vonmises is as
follows:
0 - 1
1 - 2 - 3
3 - 4
4
4 - 5
5 - 6
6 - 7
7 - 8
8 - 9
9 - 10
10 - 11
11 - 12
12 - 13
13 - 14 - 15 - 16
16
16 - 17
There does not appear to be a cell output parameter which would have made
things very easy and I think may be easy to implement. It appears that all
the other parameters are from cell data so I suspect it is readily
available. I have not looked closely at the DataPostprocessorScalar(Vector)
function, but is adding the cell parameter something very easy to do? If I
modify the appropriate code in my dealii-build directory and recompile my
program, will the appropriate function(s) be updated. I could then I think
use user_index or user_pointer. I know I could also do this by writing and
threading my own function but would have to know the output order of the
cells and vertices which is not obvious. Is there another way to do it that
I'm missing? Also I searched all the step-xx looking for material keyword
and none post-process with the different materials.
Thanks again beforehand.
Pete Griffin
/****************************************************************************************
*
* class ComputeVonMises : public DataPostprocessorScalar<dim>
*
* Computes the scalar, Von-Mises Stress, and the location of Max Von-Mises
Stress
*
****************************************************************************************/
template <int dim>
class ComputeVonMises : public DataPostprocessorScalar<dim>
{
public:
ComputeVonMises();
// virtual
void
compute_derived_quantities_vector
(const std::vector<Vector<double> > &uh,
const std::vector<std::vector<Tensor<1, dim> > > &duh,
const std::vector<std::vector<Tensor<2, dim> > > &dduh,
const std::vector<Point<dim> > &normals,
const std::vector<Point<dim> >
&evaluation_points,
std::vector<Vector<double> >
&computed_quantities) const;
};
template <int dim>
ComputeVonMises<dim>::ComputeVonMises ()
:
DataPostprocessorScalar<dim> ("VonMisesStress", update_values |
update_gradients | update_quadrature_points)
{
}
/****************************************************************************************
*
****************************************************************************************/
template <int dim>
void ComputeVonMises<dim>::
compute_derived_quantities_vector (const std::vector<Vector<double> >
&uh,
const
std::vector<std::vector<Tensor<1,dim> > > & duh,
const
std::vector<std::vector<Tensor<2,dim> > > & /*dduh*/,
const std::vector<Point<dim> >
& /*normals*/,
const std::vector<Point<dim> >
& evaluation_points,
std::vector<Vector<double> >
&computed_quantities) const
{
SymmetricTensor<2, 3> strain, stress;
SymmetricTensor<4, 3> C;
double VonMisesStress;
const unsigned int n_quadrature_points = uh.size();
Assert (duh.size() == n_quadrature_points,
ExcInternalError());
Assert (computed_quantities.size() == n_quadrature_points,
ExcInternalError());
Assert (uh[0].size() == dim,
ExcInternalError());
#if 0
fprintf(stderr, "n_quadrature_points %d\n", n_quadrature_points);
for (unsigned int q=0; q < n_quadrature_points; ++q) {
fprintf(stderr, "%7.4f %7.4f %7.4f\n", evaluation_points[q][0],
evaluation_points[q][1], evaluation_points[q][2]);
}
#endif
#if 1
C = GetStressStrainTensor3D (3846153.846154, 5769230.769231);
fprintf(stderr, "%3d - ", iCELL++);
fflush(stderr);
#else
C = GetStressStrainTensor3D(lambda[iCELL], mu[iCELL]);
fprintf(stderr, "%3d mu %f lambda %f\n", iCELL-1, mu[iCELL-1],
lambda[iCELL-1]);
fflush(stderr);
#endif
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
//#### NaN test
for (unsigned int i=0; i<dim; i++)
for (unsigned int j=0; j<dim; j++)
if (std::isnan(duh[q][i][j])) fprintf(stderr, "duh[%d][%d][%d] =
%f", q, i, j, duh[q][i][j]);
for (unsigned int i=0; i<dim; i++) { // step-18
strain[i][i] = duh[q][i][i];
for (unsigned int j=i+1; j<dim; j++) { // j=i+1 since
SymmetricTensor<2, dim>
strain[i][j] = (duh[q][i][j] + duh[q][j][i])/2;
}
}
stress = C*strain;
// Von-Mises Stress
double vonSquared =
(6*(stress[0][1]*stress[0][1]+stress[1][2]*stress[1][2]+stress[2][0]*stress[2][0])
+
(stress[0][0]-stress[1][1])*(stress[0][0]-stress[1][1]) +
(stress[1][1]-stress[2][2])*(stress[1][1]-stress[2][2]) +
(stress[2][2]-stress[0][0])*(stress[2][2]-stress[0][0]))/2;
VonMisesStress = sqrt(vonSquared);
if (VonMisesStress > MaxVonMises) {
MaxVonMises = VonMisesStress;
MaxVonMisesLoc = evaluation_points[q];
}
computed_quantities[q](0) = VonMisesStress;
//#### NaN test
if (std::isnan(VonMisesStress)) fprintf(stderr, "VonMisesStress[%d]
= %f\n", q, VonMisesStress);
}
fprintf(stderr, "%3d\n", iCELL);
}
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
