Output nodal scalar

[Jie Cheng]

Hi everyone

I'm trying to compute and output nodal stress and strain. Although there are already some helpful discussions, for example this and this, none of them addresses my problem.

I computed the stress and strain on the quadrature points in each cell and used FETools::compute_projection_from_quadrature_points_matrix to extrapolate them to the vertices as in Step-18. Eventually I end up with a vector<vector<Vector<double>>> stressGlobal, whose size is dim*dim*outputHandler.n_dofs(). Note that the outputHandler is different from the dof handler I used to solve the linear elastic problem because its number of dofs equals the number of vertices in the triangulation, rather than dim times that.

I know that I must output the global_stress component by component as scalars. The problem is, given the DataOut object that is attached to the original DoFHandler, how do I write a specific component of stressGlobal to it?  Here is what I tried:

  std::string fileName = "solution-";
  fileName += ('0' + cycle);
  fileName += ".vtu";
  std::ofstream out(fileName.c_str());
  DataOut<dim> data;
  data.attach_dof_handler(this->dofHandler);
  std::vector<std::string> solutionNames;
  switch (dim)
  {
  case 2:
    solutionNames.push_back ("x_displacement");
    solutionNames.push_back ("y_displacement");
    break;
  case 3:
    solutionNames.push_back ("x_displacement");
    solutionNames.push_back ("y_displacement");
    solutionNames.push_back ("z_displacement");
    break;
  default:
    Assert(false, ExcNotImplemented());
  }
  data.add_data_vector(this->solution, solutionNames);
  
  data.clear_data_vectors();
  data.attach_dof_handler(outputHandler);
  data.add_data_vector(stressGlobal[0][0], "Sxx");
  data.add_data_vector(stressGlobal[1][1], "Syy");
  data.add_data_vector(stressGlobal[0][1], "Sxy");
  data.build_patches();
  data.write_vtu(out);

clear_data_vectors would erase the displacement data, which is wrong. But if I do not reattach the DataOut to the outputHandler, I would get an error: the size of the vector stressGlobal does not equal n_dofs of the original DoFHandler. There must be a easy fix to this, any suggestions?

In addition, I used a DofHandler<dim> outputHandler to extrapolate stress. If I want to output the strain as well, do I have to create another DofHandler ? (Because each DofHandler has only one dof on every node.)

I hope I made my point clear. Thank you in advance!

Jie

 

[Jie Cheng]

Update:

I switched to my old trick of nodal averaging method: take the average of the values at quadrature points, then add this average to the vertex dofs. And if a vertex is shared by n cells, divide the values at the related dofs by n.

I am happy with the results because they look quite smooth.

The only concern is that: for a stress component, say Sxx, appears in my resulting vtu files as Sxx_0, Sxx_1 (in 2D) because the variable Sxx is a vector of size n_dofs as required by the add_data_vector function. Is there a way to get rid of the unnecessary duplicates? Essentially I want Sxx to be a vector of size n_vertices.

Thank you

Jie

[Jean-Paul Pelteret]

Dear Jie,

I'm glad to hear that you managed to find a solution to your problem and for posting what it was. Two quick points:

1. There's actually another add_data_vector() member function of DataOut that takes in another DoFHandler as an argument, so you could compute the stress fields using another DoFHandler, with another DataComponentInterpretation etc. This will help remedy your last concern, because you will have only as many components as you need.

2. You might want to consider solving for the smoothed stress fields via global L2-smoothing rather than manual averaging. To the best of my understanding, it should produce a more accurate representation of the stress field than your approach. This should also be a relatively inexpensive operation because solving the linear system only involves inverting a mass matrix. 

Regards,

Jean-Paul

[Jie Cheng]

Hi Jean-Paul

Thanks for your help! The specific add_data_vector() function you shared does solve my problem! For your second suggestion on L2-smoothing, I'd like to give it a try. Can you please point me to some detailed references? Examples or documentations would be great!

Jie

[Jean-Paul Pelteret]

Hi Jie,

Sure, we used it in this paper

@Article{vogel2014a-preprint,
  author    = {Vogel, F. and Pelteret, J-P. V. and Kaessmair, S. and Steinmann, P.},
  title     = {Magnetic force and torque on particles subject to a magnetic field},
  journal   = {European Journal of Mechanics A/Solids},
  year      = {2014},
  volume    = {48},
  pages     = {23--31},
  month     = {November--December}
}

but, to be honest, almost any finite element textbook should have a blurb on it. E.g. page 227 of

@Book{Hughes2000a,
  Title                    = {The Finite Element Method: Linear Static and Dynamic Finite Element Analysis},
  Author                   = {Hughes, T. J.},
  Publisher                = {Dover Publications Inc.},
  Year                     = {2000},
  Address                  = {New York, USA},
  Note                     = {ISBN: 978-0486411811},
}

Regards,

Jean-Paul

[Jie Cheng]

Jean-Paul

Any examples programming-wise? Like a deal.II manual/tutorial...

Jie

[Jean-Paul Pelteret]

I can't think of any off of the top of my head, but I think that MatrixCreator::create_mass_matrix() plus VectorTools::create_right_hand_side() collectively assemble what you're wanting. Its really not difficult to do by hand - you should just give it a go!

[Jean-Paul Pelteret]

There's also the MeshWorker framework which has some integrators that can help. But I don't know how that works so I can't help you further other than to point it out to you.

[Jie Cheng]

Thank you Jean-Paul. I did not hear about L2-projection method before, but it seems that the L2-projection method also evaluates the stress tensor at the quadrature points and then finds the cell-averaged stress. The only difference compared with the heuristic method is that it solves AX = Y where A is mass matrix, Y is integration of cell-averaged stress component and X is the nodal stress component that we want. Is that right?

[Jean-Paul Pelteret]

Dear Jie,

 

it seems that the L2-projection method also evaluates the stress tensor at the quadrature points and then finds the cell-averaged stress. The only difference compared with the heuristic method is that it solves AX = Y where A is mass matrix, Y is integration of cell-averaged stress component and X is the nodal stress component that we want. Is that right?

No, this is not correct. I can think of only two simplified situations where the least-squares projection is reduced to either of the two calculation approaches you've mentioned in this thread:

1. If you use a FE_DGQ of degree 0, then the projection is reduced to the cell-average value (a constant, and thus obviously  piece-wise discontinuous).

2. If your grid is uniform and Cartesian, and you're using linear FE_Q elements then you can compute the equivalent average (continuous) nodal value using the geometric arguments you put forward earlier.

However, if you add any complexity to your problem then immediately you run into issues:

- If your grid is non-uniform or unstructured, then your have more geometric considerations to make: If you have a large cell adjacent to a small cell, then is the weighting of the values that come from the large cell equal to the small cell?

- What if you have higher order continuous elements, and the support points no longer all coincide with the cell vertices?

- What if your finite element has no support points?

In any of these situations your geometric averaging idea will either fail or it will not produce a smoothed field that is representative of the underlying quadrature data. Like a simple linear regression, the least-squares method will minimise the error of your function (as given by the quadrature data) projected onto the given finite element space (be it continuous or discontinuous, low or high order, etc.).

Does that make sense?

Regards,

Jean-Paul

[Jie Cheng]

Dear Jean-Paul

Now I understand why the heuristic method is less preferable and have a better sense of the least square method in this context.

Like a simple linear regression, the least-squares method will minimise the error of your function (as given by the quadrature data) projected onto the given finite element space (be it continuous or discontinuous, low or high order, etc.).

Can I  express this method as following: given the scalar data lower case p_j at the quadrature points, find the nodal values upper case P_i such that sum(p'_j - p_j)^2 is minimum where p'_j is the values at the quadrature points obtained by interpolating P_i ?

I did not derive, instead learned from this link, that the least square method yields the following equation:

MP = R

where M is mass matrix, P is nodal values, and R is the volume integration of the scalar variable p.

This involves solving additional linear equation but the equation itself is easy to assemble. Is my understanding right?

Thank you so much

Jie

[Jean-Paul Pelteret]

Yes, thats correct :-) The book link that you mention outlines the problem exactly as I understand it, and I interpret it in the same way as you've mentioned.

[Wolfgang Bangerth]

> I'm trying to compute and output nodal stress and strain. Although there
> are already some helpful discussions, for example this

> <https://groups.google.com/forum/#!searchin/dealii/project$20stress$20and$20strain|sort:relevance/dealii/lvwwrT6zOFQ/C4r3JVAQAQAJ>
> and this
> <https://groups.google.com/forum/#!searchin/dealii/output|sort:relevance/dealii/uME5kwaJIEU/s5EN7hutAAAJ>,

> none of them addresses my problem.
>
> I computed the stress and strain on the quadrature points in each cell
> and used FETools::compute_projection_from_quadrature_points_matrix

> [...]

This is one way of doing it, but an easier way -- at least for
visualization purposes -- is to use the DataPostprocessor framework to
directly output the gradient (or the stress, if you want) to the viz files.

Rather than describing this here, I thought I'd take the time to
actually write this up properly since your question is not the first one
to go in this direction. Take a look here:
https://github.com/dealii/dealii/pull/4991

This is not quite what you want: I only show how to do it for a scalar
solution where the gradient is a vector. In your case, the solution is a
vector and the gradient/stress a tensor. We don't currently support
tensor-valued output, so you'll have to output dim^2 components via the
DataPostprocessor class, but maybe you get the idea.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Jie Cheng]

Hi Wolfgang

Thank you very much for adding this documentation! I'm trying to use the same technique on my linear elasticity problem. To compute the displacement gradient, the idea is to extract ux, uy, uz from the solution and then compute their individual gradients as you suggested. However I have some troubles:

My gradient processor is exactly the same as the one your wrote in include/deal.II/numerics/data_postprocessor.h and the way I compute the gradient of ux is like this:

  Vector<double> ux(this->solution.size()/dim);

  for (unsigned int i = 0; i < ux.size(); ++i)

  {

    ux[i] = this->solution[i];

  }

  GradientPostProcessor<dim> gradientProcessor;

  data.add_data_vector(scalarDofHandler, ux, gradientProcessor);

 

where I assume the sequence of the solution is ux, uy, uz.. I guess this is not the right way to do it? And the scalarDofHandler is constructed just for the purpose of outputing nodal scalars: 

FE_Q<dim> scalarFe(1); 

DoFHandler<dim> scalarDofHandler(this->tria); 

scalarDofHandler.distribute_dofs(scalarFe); 

When I run this code, I got:

 

--------------------------------------------------------

An error occurred in line <104> of file </home/jie/dealii/source/base/subscriptor.cc> in function

    void dealii::Subscriptor::check_no_subscribers() const

The violated condition was: 

    counter == 0

Additional information: 

    (none)

Stacktrace:

-----------

#0  /home/jie/dealii-build/lib/libdeal_II.g.so.9.0.0-pre: dealii::Subscriptor::check_no_subscribers() const

#1  /home/jie/dealii-build/lib/libdeal_II.g.so.9.0.0-pre: dealii::Subscriptor::~Subscriptor()

#2  /home/jie/dealii-build/lib/libdeal_II.g.so.9.0.0-pre: dealii::DataPostprocessor<2>::~DataPostprocessor()

#3  ./main: dealii::DataPostprocessorVector<2>::~DataPostprocessorVector()

#4  ./main: IFEM::LinearElasticSolver<2>::output(unsigned int) const

#5  ./main: main

--------------------------------------------------------

 I pulled the latest master from dealii yesterday, and pass Test #2976 data_out_postprocessor_vector_03. My code crashed at add_data_vector. Can you see why this happened?

Thank you

Jie

[Wolfgang Bangerth]

> My gradient processor is exactly the same as the one your wrote in
> include/deal.II/numerics/data_postprocessor.h

> <https://github.com/dealii/dealii/pull/4991/files#diff-0e837098279c080004cb4372efaf37e6> and

> the way I compute the gradient of ux is like this:
>
>
>   Vector<double> ux(this->solution.size()/dim);
>
>   for (unsigned int i = 0; i < ux.size(); ++i)
>
>   {
>
>     ux[i] = this->solution[i];

Yes, this does not work. You cannot just take one block out of the block
vector because you do not have a DoFHandler that describes only this one
block (your DoFHandler describes the collection of all elements of all
blocks).

You will want to output *all* gradients at once. It turns out that I
just wrote a class for that that you may want to copy, along with the
example given there:

https://github.com/dealii/dealii/pull/4997

[Jie Cheng]

Hi Wolfgang

I tried DataPostprocessorTensor as soon as the commits merged in. It is easy to use, I really like it! Thank you! Unfortunately there is still a bug that I couldn't figure out, which is basically the same as my previous post: after successfully writing the output, the program aborts with the following message:

--------------------------------------------------------

An error occurred in line <104> of file </home/jie/dealii/source/base/subscriptor.cc> in function

    void dealii::Subscriptor::check_no_subscribers() const

The violated condition was: 

    counter == 0

Additional information: 

    (none)

Stacktrace:

-----------

#0  /home/jie/dealii-latest/lib/libdeal_II.g.so.9.0.0-pre: dealii::Subscriptor::check_no_subscribers() const

#1  /home/jie/dealii-latest/lib/libdeal_II.g.so.9.0.0-pre: dealii::Subscriptor::~Subscriptor()

#2  /home/jie/dealii-latest/lib/libdeal_II.g.so.9.0.0-pre: dealii::DataPostprocessor<2>::~DataPostprocessor()

#3  ./main: IFEM::LinearElasticSolver<2>::output(unsigned int) const

#4  ./main: IFEM::LinearElasticSolver<2>::runStatics(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&)

#5  ./main: main

--------------------------------------------------------

Aborted (core dumped)

This looks like the early destruction bug in step 6, but I couldn't think of any similar mistake in my program. I basically copied your StrainPostprocessor class in data_out_postprocessor_tensor_01.cc, and used it like this:

    std::vector<DataComponentInterpretation::DataComponentInterpretation>
      interpretation(dim, DataComponentInterpretation::component_is_part_of_vector);
    data.add_data_vector(this->dofHandler, this->solution,
      std::vector<std::string>(dim, "displacement"), interpretation);
    StrainPostprocessor<dim> strain;
    data.add_data_vector(this->dofHandler, this->solution, strain);
    
    data.build_patches();
    data.write_vtu(out);

 

Thanks

Jie

[Jean-Paul Pelteret]

Dear Jie,

Your query relates to this FAQ post. You’ll likely find that DataOut is holding a pointer to the DataPostprocessorTensor, and the post processor is going out of scope (so being removed from the stack) before DataOut. So I think that if you switch the order of construction of these two objects, ie.

StrainPostprocessor<dim> strain;

DataOut<dim> data;

…

then this would probably fix the problem.

Best,

J-P

--
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 dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jie Cheng]

Hi Jean-Paul

Thanks for explaining it. Now I got both StrainPostprocessor and StressPostprocessor working! But it seems a little wasteful to have two separate classes because we are not reusing the gradient of displacement computed in StrainPostprocessor. Although it is tolerable, but it'll be even better if there is a way to output both in one class derived from DataPostprocessorTensor. (:

Jie

[Wolfgang Bangerth]

It can't be derived from DataPostprocessorTensor because that class is only
meant to output a single tensor, whereas you want to output two tensors. But
it should not be very hard to write a class derived from DataPostprocessor
that outputs both the strain and the stress.

In practice, I suspect that you will not be able to measure much of a
difference between using two classes and merging it all into one.