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
Am Fr., 25. Okt. 2019 um 18:18 Uhr schrieb Wolfgang Bangerth <
bange...@colostate.edu>:
> 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: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> 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+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/81e2c8f1-a2f5-f026-5a57-c59ec20b26e7%40colostate.edu
> .
>
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/CAF_nQ7AAbjr4Ar_EOHA_OOz3PLaj5W%3DvRZxYVTu5YUOtV5qkxA%40mail.gmail.com.
