Output of Gauss point stress tensor

[Muhammad Mashhood]

Dear Deal.ii users,

                                Hi! I am working on the deal.ii code. I have data of stress and strain tensors at quadrature points. Currently using the results_output function, I am transferring the data at nodal points of my triangulation using " qpoint_to_dof_matrix.vmult " function and mapping "local_history_strain_values_at_qpoints[i][j]" Gauss point values on the "history_strain_field[i][j]" tensor with "dg_cell->set_dof_values" function. After this I write this nodal values data "history_strain_field[i][j]" as .pvtu and .vtu files for Paraview.

My question is that along with this, can I also export the stress and strain tensors data of quadrature points i.e. "local_history_strain_values_at_qpoints[i][j]" directly to the output file (.pvtu or .vtu etc.) i.e. without mapping or averaging etc. on the node. In this way I want to directly visualize the quadrature point values of stress and strain in Paraview.

My current simplified code is following:

FE_DGQ<dim> history_fe (1);
    DoFHandler<dim> history_dof_handler (triangulation);
    history_dof_handler.distribute_dofs (history_fe);

std::vector< std::vector< Vector<double> > >
    history_strain_field (dim, std::vector< Vector<double> >(dim)),
                         local_history_strain_values_at_qpoints (dim, std::vector< Vector<double> >(dim)),
                         local_history_strain_fe_values (dim, std::vector< Vector<double> >(dim));

history_strain_field[i][j].reinit(history_dof_handler.n_dofs());
                    local_history_strain_values_at_qpoints[i][j].reinit(quadrature_formula.size());
                    local_history_strain_fe_values[i][j].reinit(history_fe.dofs_per_cell);

for (unsigned int q=0; q<quadrature_formula.size(); ++q)
            {

for (unsigned int i=0; i<dim; i++)
                for (unsigned int j=i; j<dim; j++)
                  {

local_history_strain_values_at_qpoints[i][j](q) = local_quadrature_points_history[q].old_elastic_strain;

} }

for (unsigned int i=0; i<dim; i++)
            for (unsigned int j=i; j<dim; j++)
              {

qpoint_to_dof_matrix.vmult (local_history_strain_fe_values[i][j],
                                                            local_history_strain_values_at_qpoints[i][j]);

                                dg_cell->set_dof_values (local_history_strain_fe_values[i][j],
                                                         history_strain_field[i][j]);
}

DataOut<dim>  data_out;
              data_out.attach_dof_handler (history_dof_handler);


              //for strain
              data_out.add_data_vector (history_strain_field[0][0], "strain_xx");
              data_out.add_data_vector (history_strain_field[1][1], "strain_yy");
              data_out.add_data_vector (history_strain_field[0][1], "strain_xy");

data_out.build_patches ();


              const std::string filename_base_strain = ("strain-" + filename_base);



              const std::string filename =
                          (output_dir + filename_base_strain + "-"
                           + Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4));

              std::ofstream output_vtu((filename + ".vtu").c_str());
              data_out.write_vtu(output_vtu);


              pcout << output_dir + filename_base_strain << ".pvtu" << std::endl;

              if (this_mpi_process == 0)
                {
                  std::vector<std::string> filenames;
                                for (unsigned int i = 0; i < n_mpi_processes; ++i)
                                    filenames.push_back(filename_base_strain + "-" +
                                                      Utilities::int_to_string(i, 4) +
                                                      ".vtu");

                                std::ofstream pvtu_master_output((output_dir + filename_base_strain + ".pvtu").c_str());
                                data_out.write_pvtu_record(pvtu_master_output, filenames);

                                std::ofstream visit_master_output((output_dir + filename_base_strain + ".visit").c_str());
                                data_out.write_pvtu_record(visit_master_output, filenames);
                }


            }

[Wolfgang Bangerth]

On 7/15/20 11:05 AM, Muhammad Mashhood wrote:
>
> My question is that along with this, can I also export the stress and strain
> tensors data of quadrature points i.e.
> "local_history_strain_values_at_qpoints[i][j]" directly to the output file
> (.pvtu or .vtu etc.) i.e. without mapping or averaging etc. on the node. In
> this way I want to directly visualize the quadrature point values of stress
> and strain in Paraview.

The question is more a philosophical one than one of how you can achieve this.
Typically, when we output information in finite element contexts, we output
them as "fields", i.e., functions of space. This allows us to show them as
surfaces, with color gradients, etc. The strategy you have found of taking
values defined in (quadrature) points and converting them to fields is a way
to make that happen.

If you want to output the stress/strain information, you have two options:

* You actually do show them as values defined only at individual points,
rather than as fields. In this scheme, the stress/strain really only exists
at the quadrature points, but not anywhere in between. You can't create
surface plots, you can't create isocontour plots, etc.
* You create a field that somehow approximates these point values. There are
different ways of doing that, and the approach you are currently using is
one way for this.

So the question is mostly: What's your goal with this? Do you want to show
these quantities as fields, or as data to be visualized at individual points?

Best
W.


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

[Muhammad Mashhood]

Dear Prof. Wolfgang,

                                      Thank you for the comprehensive and informative response. As my goal is to access the original Gauss point values and I do not want to involve any interpolation or approximation scheme so I would like to opt for the first option among the two you mentioned:

"If you want to output the stress/strain information, you have two options:

* You actually do show them as values defined only at individual points,
rather than as fields. In this scheme, the stress/strain really only exists
at the quadrature points, but not anywhere in between. You can't create
surface plots, you can't create isocontour plots, etc.
* You create a field that somehow approximates these point values. There are
different ways of doing that, and the approach you are currently using is
one way for this."

 

So far it would be enough if I have the Gauss point values only at the points rather than having complete field. I like to access these point values of stress and strain tensors in .vtk file or .pvtu file through ParaView (same as I am accessing the results currently).

I would be grateful if you or any other user has experience in this regard or know about relevant deal.ii tutorial dealing with same feature. Thank you!

[Wolfgang Bangerth]

On 7/16/20 9:15 AM, Muhammad Mashhood wrote:
> So far it would be enough if I have the Gauss point values only at the points
> rather than having complete field. I like to access these point values of
> stress and strain tensors in .vtk file or .pvtu file through ParaView (same as
> I am accessing the results currently).
> I would be grateful if you or any other user has experience in this regard or
> know about relevant deal.ii tutorial dealing with same feature.

deal.II has functionality to output data on individual points if that data is
associated with particles. I don't know how much work you want to invest in
getting this to work for you, but here are two options:

* Easy: Create a ParticleHandler object, loop over all of your quadrature
points, and for each quadrature point you create a particle. You can then
associate "properties" with each particle, and use Particles::DataOut to
output these properties in the same way as you would use DataOut for field
data. You can probably figure out how this works by looking at the draft
step-19 tutorial program:
https://github.com/dealii/dealii/pull/10301

* More work, but more elegant: You could write a class, let's say
QuadraturePointData::DataOut or something similar, to which you could describe
the information of location + values without the detour of particles. I'm sure
we'd be happy to walk you through the process of doing so if you wanted to go
for that, and it would be a very nice addition to deal.II if you wanted to get
it merged.

[Muhammad Mashhood]

Dear Prof. Wolfgang,

                                      Thank you very much for the suggestions. To have a quick solution, I am opting the first suggestion i.e. to make particles on quadrature points and assigning them the tensor quantities. And if I gain some reasonable results specific to my project then I might also want to add it as a permanent feature in my code by spending more time to implement the second suggested solution of elegant way. So I have updated my version from 9.1.1 to 9.2.0 and started writing the code for particles.

My original code is similar as step-42 but first I am using step-42 to fulfill at least the task of creating particles on quadrature point locations.

I have declared the concerned variables:

    Particles::ParticleHandler<dim> particle_handler;
    types::particle_index           next_unused_particle_id;

    MappingQ<dim>             mapping;

and then added the following syntax during construction of  class (PlasticityContactProblem<dim>::PlasticityContactProblem) (after line 833):

    , mapping(1)
    , particle_handler(triangulation, mapping, /*n_properties=*/dim)
    , next_unused_particle_id(0)

When I run the step-42 then it returns following error:

.~/working_dir/step-42$ /step-42 p1_chinese.prm
Segmentation fault (core dumped)

I tried to comment out the line "  , particle_handler(triangulation, mapping, /*n_properties=*/dim) " and by doing so the code starts running normally.

I have tried to run it in the debug mode but over there the program exits with following error:

(gdb) r
Starting program: ~/working_dir/step-42/step-42
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
*** Call this program as <./step-42 input.prm>
[Inferior 1 (process 5082) exited with code 01]

Can there be a possible problem in using the " particle_handler(triangulation, mapping, /*n_properties=*/dim) " ? Thank you!

[Wolfgang Bangerth]

On 7/28/20 11:05 AM, Muhammad Mashhood wrote:
>
> /.~/working_dir/step-42$ /step-42 p1_chinese.prm
> Segmentation fault (core dumped)/

>
> I tried to comment out the line "  , particle_handler(triangulation, mapping,
> /*n_properties=*/dim) " and by doing so the code starts running normally.
>
> I have tried to run it in the debug mode but over there the program exits with
> following error:
>

> /(gdb) r

> Starting program: ~/working_dir/step-42/step-42
> [Thread debugging using libthread_db enabled]
> Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
> *** Call this program as <./step-42 input.prm>

You need to run the executable with
r p1_chinese.prm
to make sure it loads the appropriate input file, or you will not see the
segmentation fault in the debugger.


> [Inferior 1 (process 5082) exited with code 01]/

>
> Can there be a possible problem in using the "

> /particle_handler(triangulation, mapping, /*n_properties=*/dim) /" ? Thank you!

I don't know. Can you get a backtrace in the debugger?

It's something that ought to work. If you can't figure out what is happening,
send me the exact .cc file you're trying to run and I'll take a look!

[Muhammad Mashhood]

Dear Prof. Wolfganag,

                                       Thank you for the suggestion. I ran debugger with the ".prm" file and was successful to backtrace the error and resolved it. After this, using the step-19 particles tutorial refered by you, so far I am also successful in creating and exporting out the particles in ".vtu" file to visualize their ids and location (which is quadrature points location in my case).

As a next step, After creating particle on quadrature point, I am now trying to assign it the quadrature point value which is "SymmetricTensor<2, dim>" in my case. I am using following simplified function similar to the one in particles tutorial file of step-19:

template <int dim>
    void  Problem<dim>::create_particles()
    {
      next_unused_particle_id = 1;

        FEValues<dim> fe_values (fe, quadrature_formula,
                                     update_values | update_gradients |
                                     update_quadrature_points);

        for (const auto &cell : dof_handler.active_cell_iterators())
              {
            PointHistory<dim> *local_quadrature_points_history
                               = reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
                                      Assert (local_quadrature_points_history >=
                                              &quadrature_point_history.front(),
                                                  ExcInternalError());
                                   Assert (local_quadrature_points_history <
                                           &quadrature_point_history.back(),
                                                  ExcInternalError());

                   fe_values.reinit(cell);

                for (const unsigned int q_point :
                        fe_values.quadrature_point_indices())
                  {
                        const Point<dim> location =
                                fe_values.quadrature_point(q_point);

                        Particles::Particle<dim> new_particle;
                        new_particle.set_location(location);
                        new_particle.set_reference_location(
                          mapping.transform_real_to_unit_cell(cell, location));
                        new_particle.set_id(next_unused_particle_id);

                        SymmetricTensor<2, dim> total_strain = local_quadrature_points_history[q_point].old_strain;

                        new_particle.set_properties(make_array_view(total_strain));

                        particle_handler.insert_particle(new_particle, cell);

                        ++next_unused_particle_id;
                  }
              }
        particle_handler.update_cached_numbers();
    }

But it gives the following error on running:

An error occurred in line <298> of file </home/muhammad/software/dealii-9.2.0/source/particles/particle.cc> in function
    void dealii::Particles::Particle<dim, spacedim>::set_properties(const dealii::ArrayView<const double>&) [with int dim = 3; int spacedim = 3]
The violated condition was:
    property_pool != nullptr
Additional information:
    This exception -- which is used in many places in the library -- usually indicates that some condition which the author of the code thought must be satisfied at a certain point in an algorithm, is not fulfilled. An example would be that the first part of an algorithm sorts elements of an array in ascending order, and a second part of the algorithm later encounters an element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an exception since it indicates an error in deal.II, not in your own program. Try to come up with the smallest possible program that still demonstrates the error and contact the deal.II mailing lists with it to obtain help.

I tried to fix it by initializing and setting the property pool to the new_particle in my current function :

Particles::PropertyPool propertypool(6); // six entries present in the symmetric tensor of second order

new_particle.set_property_pool(propertypool);

But still it stops with the following error:

An error occurred in line <531> of file </home/muhammad/software/dealii-9.2.0/include/deal.II/base/array_view.h> in function
    dealii::ArrayView<ElementType, MemorySpaceType>::value_type& dealii::ArrayView<ElementType, MemorySpaceType>::operator[](std::size_t) const [with ElementType = double; MemorySpaceType = dealii::MemorySpace::Host; dealii::ArrayView<ElementType, MemorySpaceType>::value_type = double; std::size_t = long unsigned int]
The violated condition was:
    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(i), decltype(n_elements)>::type)>:: type>(i) < static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(i), decltype(n_elements)>::type)>:: type>(n_elements)
Additional information:
    Index 3 is not in the half-open range [0,3).

I explored the mailing list and found someone from the community already had similar problem of assigning the properties to the particles where it was suggested to use following way of assigning the properties:

for( auto iter=particleHandler.begin(); iter!=particleHandler.end();
++iter, ++part )

iter->set_properties
(make_array_view(interpolatedParticleQuantities.begin()+
part*ncomponents,
interpolatedParticleQuantities.begin()+
part*ncomponents + ncomponents));

If I have related it correctly to my case, I also then tried the similar way e.g.

  

new_particle.set_properties(make_array_view(total_strain.begin_raw(),
total_strain.end_raw()));

instead of new_particle.set_properties(make_array_view(total_strain));

It still stops with the same (2nd) error as mentioned above. Is there another way of doing it or am I not using the correct approach as per my limited knowledge so far? Sorry for the lengthy thread, I wanted to clear out the approaches which I tried so far as much as possible.

Thank you very much for your ongoing and productive support.

Regards,

Muhammad 

[Wolfgang Bangerth]

On 7/30/20 5:26 PM, Muhammad Mashhood wrote:
>
> _But it gives the following error on running:_
>
> /An error occurred in line <298> of file

> </home/muhammad/software/dealii-9.2.0/source/particles/particle.cc> in function
>     void dealii::Particles::Particle<dim, spacedim>::set_properties(const
> dealii::ArrayView<const double>&) [with int dim = 3; int spacedim = 3]
> The violated condition was:
>     property_pool != nullptr
> Additional information:

> [...]

> I explored the mailing list and found someone from the community already had
> similar problem of assigning the properties to the particles where it was
> suggested to use following way of assigning the properties:

I believe that this is a bug that I fixed a while back:
https://github.com/dealii/dealii/pull/10589
https://github.com/dealii/dealii/issues/10590
Unfortunately, this happened after the 9.2 release. Are you in a position to
work with the current development sources?

If this doesn't solve the problem, can you come up with a simplified piece of
code that illustrates the issue?

[Muhammad Mashhood]

Dear Prof. Wolfgang,

                                     Thanks for the guidance. I tried replacing the " source/particles/particle_handler.cc" file with the one present in "https://github.com/dealii/dealii/pull/10589/files". Also include/deal.II/particles/particle_accessor.h and source/particles/particle_accessor.cc from "https://github.com/dealii/dealii/pull/10319/files". Then compiled again my deal.ii 9.2.0 but it gives error during compilation.

Therefore as a second option, I have created a simplified code which represents my problem. Kindly receive the files attached. Looking forward for your guiding response!

Thank you very much for ongoing cooperation.

Best regards,

Muhammad Mashhood

[Wolfgang Bangerth]

>                                      Thanks for the guidance. I tried
> replacing the " source/particles/particle_handler.cc

> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fpull%2F10589%2Ffiles%23diff-df7869b8ff6741c5bba620988f7bd995&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cd00ea92f27a5466d8a5e08d8395d0f41%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637322419138637595&sdata=CQSVVYYrHp5wMfy3INtkJBdS%2BNjFzBqLAdPwVHI1r%2FU%3D&reserved=0>"

> file with the one present in
> "https://github.com/dealii/dealii/pull/10589/files". Also
> include/deal.II/particles/particle_accessor.h

> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fpull%2F10319%2Ffiles%23diff-ba98c140727c3e55b479f34172de381b&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cd00ea92f27a5466d8a5e08d8395d0f41%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637322419138637595&sdata=Y9F5eICnrFT8jNpnmg03tO3gb%2BWO%2B8xEqt9%2F9fucTBQ%3D&reserved=0>
> and source/particles/particle_accessor.cc
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fpull%2F10319%2Ffiles%23diff-0b2cb32ab20618cd2a77e2ddfda9eddd&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cd00ea92f27a5466d8a5e08d8395d0f41%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637322419138647592&sdata=7ZHhT3M4M%2Bu2FaPwpR2qwSkJevjIip4d7%2BvaQrqmsG0%3D&reserved=0>

> from "https://github.com/dealii/dealii/pull/10319/files". Then compiled again
> my deal.ii 9.2.0 but it gives error during compilation.

Yes -- we put about 10 changes every day into deal.II. You can't just replace
individual files in 9.2 with the current development sources :-)

> Therefore as a second option, I have created a simplified code which
> represents my problem. Kindly receive the files attached. Looking forward for
> your guiding response!

The code in question looks like this:

Particles::Particle<dim> new_particle;
new_particle.set_location(location);
new_particle.set_reference_location(
mapping.transform_real_to_unit_cell(cell, location));
new_particle.set_id(next_unused_particle_id);

SymmetricTensor<2, dim> strain; strain = 0;
new_particle.set_properties(make_array_view(strain));

particle_handler.insert_particle(new_particle, cell);

It is correct that set_properties() throws an exception here, because there
really is no property pool associated with this particle. If you write it like
this:

auto inserted_particle
= particle_handler.insert_particle(new_particle, cell);
inserted_particle->set_properties(make_array_view(strain));

then things will work if you also change the number of properties stored by
the ParticleHandler object to dim*(dim+1)/2=6 (it is currently 'dim').

[Muhammad Mashhood]

Thanks a lot Prof. Wolfgang. Your guidance was very fruitful and right on point. Solved my problem! :)