output results at q_points

[Frederik S.]

Hello!

I have the following problem: To output my results I use the function I attached. The first part gives me a vector with three components (the displacements in x-,y- and z-direction) for every quadrature point in the geometry. In the second part I calculate the corresponding von mises strain at every quadrature point using the function get_mises_eps, that gives me the scalar valued von mises strain corresponding to the displacement field. If I use a vector with one entry for every cell (mises_e) and calculate the mean value over the cell (dividing the sum of weighted values at the q_point through the sum of weights) this works without any problem and I get correct results. But if I now try to write out the solution at every q_point, using a vector with one entry for every q_point in the geometry (n_cells times n_qpoints_per_cell) (mises_e_2), the program runs without errors, but paraview shows me three output vectors (E_Mises_2_0, E_Mises_2_1, E_Mises_2_2) with one component per q_point although I only attach one scalar-valued vector named E_Mises_2. The values in these three vectors don't seem to make much sense and differ from E_Mises about a factor from 0,1 to 10.

Does anyone know where I did make a mistake?

Thanks in advance,

greetings,

Frederik

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

template <int dim> void OutputResults<dim>::output_volume (const unsigned int timestep){

std::vector<std::string> solution_names (dim, "displacement");

std::vector<DataComponentInterpretation::DataComponentInterpretation>  data_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector);

DataOut<dim> data_out;

data_out.attach_dof_handler (dof_handler);

data_out.add_data_vector (solution, solution_names, DataOut<dim>::type_dof_data,data_component_interpretation);

        const unsigned int n_cells = triangulation.n_active_cells();

Vector<double> mises_e (n_cells);

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

const unsigned int   n_q_points    = quadrature_formula.size();

Vector<double>  mises_e_2(n_cells*n_q_points);

int index=0;

int index2=0;

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();

for (; cell!=endc; ++cell,++index, index2=index2+n_q_points){    

fe_values.reinit (cell);

std::vector<std::vector<Tensor<1,dim> > > gradients (n_q_points, std::vector<Tensor<1,dim> >(dim));

std::vector<Vector<double> > values (n_q_points, Vector<double> (dim));

fe_values.get_function_gradients (solution, gradients);

fe_values.get_function_values    (solution, values);

double vol=0;

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

vol+=fe_values.JxW (q_point);

mises_e(index)+=get_mises_eps(gradients[q_point])*fe_values.JxW (q_point);

mises_e_2(index2+q_point)=get_mises_eps(gradients[q_point]);

}

mises_e(index)/=vol;

}

data_out.add_data_vector (mises_e_2, "E_Mises_2");

data_out.add_data_vector (mises_e, "E_Mises");

        data_out.build_patches ();

std::string filename = "solution/solution." + Utilities::int_to_string(timestep,4)+ ".vtu";

std::ofstream output (filename.c_str());

        data_out.write_vtu (output);

}

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

[Jean-Paul Pelteret]

Hi Frederik,

The problem is likely that you have a vector-valued problem (from line 1 your solution has dim components) and you are attaching a solution with only 1 component (namely mises_e_2). The output (solution) vector either has to either correspond directly with that associated with the DoFHandler, or can have the size of the number of cells (i.e. a constant value per cell). 

To fix this you could in general create another DoFHandler and pass this new DoFHandler at the same time as you attach the scalar-valued solution to the DataOut object (see this function). However since you've got data at the quadrature points that you're needing to output, you're going to have to be a little more careful how you do this (since the position of the data does not necessarily align with the support points of the FE). You can use a standard FE_DGQ and extrapolate the data to the support points on an element-wise basis. Alternatively, you might be able to use FE_DGQ elements with support points defined at the QPs (see this specialised constructor), but maybe someone else had best comment on whether or not this is the best way to approach the problem. 

I hope that this helps,

J-P

[Frederik S.]

Hi Jean Paul,

I am very, verry sorry that it took me so long to answer! 

While I was trying to implement what you suggested, I found an other mistake I made, that was more serious than problems with output, so I had to solve that first.

I now think you're right, that it might be troublesome to put out qpoint-data. While trying out your first suggestion I stumbled across the DataPostprocessor structure and tried to implement that (following step-29; the ComputeOutput-class I posted further down). It gives me some good results, if I try to give out the gradients in x-/y-/z-direction of my 3D-displacement vector solution, but the results are only good perpendicular to the gradients main direction. So for example the x-component of vector computed_quantities (coming from duh[i][0][0]) has smooth transitions between cells in the y-z-plane, but it has jumps between cells in x-direction (like I calculated the mean value over the cells x-diection). I attached an image to illustrate what I tried to describe.

Do you (or someone else here) know an answer to my problem?

Greetings and thanks in advance,

Frederik

#################################################################

a part of my output_results - function:

. . .

DataOut<dim> data_out;

data_out.attach_dof_handler (dof_handler);

data_out.add_data_vector (solution, solution_names,DataOut<dim>::type_dof_data,data_component_interpretation);

ComputeOutput<dim> grads;

data_out.add_data_vector (solution, grads);

data_out.build_patches ();

. . .

#################################################################

#################################################################

the ComputeOutput - class:

template <int dim>

class ComputeOutput : public DataPostprocessorVector<dim>

{

public:

  ComputeOutput ();

  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> ComputeOutput<dim>::ComputeOutput ()  

    :  DataPostprocessorVector<dim> ("EMises",  update_values|update_gradients) {}

template <int dim>

void ComputeOutput<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

{

  Assert(computed_quantities.size() == uh.size(),  ExcDimensionMismatch (computed_quantities.size(), uh.size()));

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

    {

      Assert(computed_quantities[i].size() == 1, ExcDimensionMismatch (computed_quantities[i].size(), 1));

      Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));

      

      computed_quantities[i](0) = duh[i][0][0];

      computed_quantities[i](1) = duh[i][1][1];

      computed_quantities[i](2) = duh[i][2][2];

    }

}

#################################################################