How to output more than one value at the location of a vertex of the mesh?

[Andreas Hegendörfer]

Hello Everyone,

I am relatively new to deal.ii.

I want to output the symmetric gradient (=strain) of the solution (=displacement).

I have implemented my "own" computation of the strain at the support points as well as using the class DataPostprocessorTensor.

A comparison gives the results that can be found in the appendix for a component of the strain.

The upper picture shows the result of DataPostprocessorTensor, while the lower one is my "own" implementation. 

According to the the FAQs (https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#the-graphical-output-files-dont-make-sense-to-me----they-seem-to-have-too-many-degrees-of-freedom) the result of DataPostprocessorTensor has more than one value at the location of a vertex of the mesh.

That makes sense because there is of course more than one gradient at each vertex. The solution looks discontinous.

I am only able to output only one value at a support point, what makes my solution looks continous (especially in the middle). The objective is to make my output of the strain to look like the ouput of DataPostprocessorTensor.

The FAQs say :"The solution to the problem is then to simply output each vertex multiple times -- with different vertex numbers but at exactly the same location, once for each cell it is adjacent to."

But I dont know how to do that. How to output more than one value at the location of a vertex of the mesh?

Thank you very much in advance,

Andreas

[Wolfgang Bangerth]

On 10/25/19 5:03 AM, 'Andreas Hegendörfer' via deal.II User Group wrote:
>
> I am only able to output only one value at a support point, what makes
> my solution looks continous (especially in the middle). The objective is
> to make my output of the strain to look like the ouput of
> DataPostprocessorTensor.
>
> The FAQs say :"The solution to the problem is then to simply output each
> vertex multiple times -- with different vertex numbers but at exactly
> the same location, once for each cell it is adjacent to."
>
>
> But I dont know how to do that. How to output more than one value at the
> location of a vertex of the mesh?

You'll have to describe how exactly you compute your strain :-) If you
output it with DataOut, then I suspect that you compute the strain
somehow and put it into a finite element field? If you want to put a
discontinuous quantity into a finite element field, then you probably
also want that finite element field to use a discontinuous element
(e.g., FE_DGQ).

Best
W.

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

[Andreas Hegendörfer]

Hello Wolfgang,

thank you very much for your response. I am sorry to send you parts of my code, but I see no other way to explain what I do.

This is how I calculate the strain:

template <int dim>
void  ElasticProblem<dim>::calculate_strain()
{
hp::QCollection<dim> q_collection_strain;
Quadrature<dim> piezo_quadrature_formula_strain (piezo_fe. get_unit_support_points()) ;
Quadrature<dim> substructure_quadrature_formula_strain (substructure_fe. get_unit_support_points()) ;


q_collection_strain.push_back (substructure_quadrature_formula_strain);
q_collection_strain.push_back (piezo_quadrature_formula_strain);


hp::FEValues<dim> hp_fe_values_strain (fe_collection,  q_collection_strain, update_values   | update_gradients | update_quadrature_points | update_JxW_values);

const FEValuesExtractors::Vector     displacement (0);
std::vector<types::global_dof_index> local_dof_indices;

std::vector<Tensor<2,dim>> strain_solution;

strain_solution.resize(dof_handler.n_dofs());
    Tensor<2,dim> strain;
typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();
cell = dof_handler.begin_active();

for (; cell!=endc; ++cell)
{
        hp_fe_values_strain.reinit (cell);
        const FEValues<dim> &fe_values = hp_fe_values_strain.get_present_fe_values();
        int   n_q_points    = fe_values.n_quadrature_points;

        local_dof_indices.resize (cell->get_fe().dofs_per_cell);

          cell->get_dof_indices (local_dof_indices);

        for (unsigned int q_point=0; q_point<n_q_points;++q_point)
        {

            strain=0;
            for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
             {
            strain=fe_values[displacement].symmetric_gradient (i, q_point)*solution [local_dof_indices[i]]+strain;
              }
            strain_solution[local_dof_indices[q_point]]=strain;
        }
}



for(int i=0; i<n_x; i++) //n_x=dofs_per_component[0];  n_x=n_y=n_z
{
strain_xx_yy_zz(i)=strain_solution[i][0][0];
}

for(int i=n_x; i<n_x+n_y; i++)  //n_y=dofs_per_component[1];
{
strain_xx_yy_zz(i)=strain_solution[i][1][1];
}


for(int i=n_x+n_y; i<n_u; i++)  //n_z=dofs_per_component[2];
{
strain_xx_yy_zz(i)=strain_solution[i][2][2];
}


}

In the programm I do have only one  hp::DoFHandler<dim>   dof_handler overall. I am using the two elements: substructure_fe (FE_Q<dim>(2), dim, FE_Nothing<dim>(), 1 ) and piezo_fe (FE_Q<dim>(2), dim, FE_Q<dim>(2), 1).

I am solving for the displacement and the potential as unknowns in a piezoelectric application. In "setup_system" I do: DoFRenumbering::component_wise(dof_handler). This means, that the dofs of the vector valued problem are ordered as follows:

- displacement in x direction

- displacement in y direction 

- displacement in z direction  

-potential.

The BlockVector<double> strain_xx_yy_zz has the size of dof_handler.n_dofs().

At first I calculate the strain at each support point. Because many support points of the vector values problem are "the same" I do some calculations more than needed. 

I save the result of the last calculation at one support point in  strain_solution[local_dof_indices[q_point]]. So the ohter calculation results are lost but it should be no effort to save them, too. 

After calculation I copy the result of the xx, yy and zz components of the strain in strain_xx_yy_zz. So for example the component xx of the strain is in the vector where the component x of the displacement is in the solution vector.

This is the relevant code for my output:

template <int dim>
void ElasticProblem<dim>::output_results_response_complex () const
{

 DataOut<dim,hp::DoFHandler<dim> >data_out;
data_out.attach_dof_handler (dof_handler);
std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation(dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar);//pot

 std::vector<std::string> solution_names_strain (0, "eigen_strain");
solution_names_strain.push_back ("eigen_xx");
solution_names_strain.push_back ("eigen_yy");
solution_names_strain.push_back ("eigen_zz");
solution_names_strain.push_back ("eigen_nothing");


std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation_strain(0, DataComponentInterpretation::component_is_scalar);
data_component_interpretation_strain.push_back(DataComponentInterpretation::component_is_scalar);//eigen_yy
data_component_interpretation_strain.push_back(DataComponentInterpretation::component_is_scalar);//eigen_zz
data_component_interpretation_strain.push_back(DataComponentInterpretation::component_is_scalar);//eigen_nothing
data_out.add_data_vector (strain_xx_yy_zz, solution_names_strain,DataOut<dim,hp::DoFHandler<dim> >::type_dof_data,data_component_interpretation);

 data_out.build_patches (2);
std::ostringstream filename;
filename << "complex.vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);

}

This all is just very provisionally. The strain components are in the entries in the solution vector, that should contain the displacement. 

I use the  class DataPostprocessorTensor for the calculation of the strain just like the documentation does.

I am really happy, that you help me. Thank you in advance!

Best,

Andreas

--
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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/81e2c8f1-a2f5-f026-5a57-c59ec20b26e7%40colostate.edu.

[Wolfgang Bangerth]

On 10/25/19 11:06 AM, 'Andreas Hegendörfer' via deal.II User Group wrote:
>
>             for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
>              {
>             strain=fe_values[displacement].symmetric_gradient (i, q_point)
> *solution [local_dof_indices[i]]+strain;

You can do the first part simpler by using
fe_values[displacement].get_function_symmetric_gradients(...)

>               }
>             strain_solution[local_dof_indices[q_point]]=strain;

So here you put 'strain' into a solution vector. This solution vector
must be associated with a finite element. And if it is associated...

> In the programm I do have only one  hp::DoFHandler<dim> dof_handler
> overall. I am using the two elements: substructure_fe (FE_Q<dim>(2),
> dim, FE_Nothing<dim>(), 1 ) and piezo_fe (FE_Q<dim>(2), dim,
> FE_Q<dim>(2), 1).

...with this element, then you will necessarily have a continuous space
into which you put the strains. If you don't like that, you'll have to
create a different (discontinuous) FE, DoFHandler, and associated vector.

> I save the result of the last calculation at one support point in
> strain_solution[local_dof_indices[q_point]]. So the ohter calculation
> results are lost but it should be no effort to save them, too.

Correct.

[Andreas Hegendörfer]

Thank you very much for your great support, Wolfgang! I know now what to do.

--
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.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1e523ccf-c490-a493-0599-d5b90e5c476d%40colostate.edu.