Hi!

I have the value of the scalar phi(x,y,z) and I would like to calculate rho(x,y,z)=- const*div(grad(phi)), so I do:
.
.
.
for (; cell_p!=endc; ++cell_p)
   {
   cell_rhs=0.0;
   fe_values.reinit (cell_p);
fe_values.get_function_2nd_derivatives (phi, values);
   for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int q_point=0;q_point<n_q_points;q_point++)
                 cell_rhs(i)-=const*fe_values.shape_value(i,q_point)*
(values[q_point][0][0]+values[q_point][1][1]+values[q_point][2][2])*fe_values.JxW(q_point);

   cell_p-> get_dof_indices(local_dof_indices);
   for (unsigned int i=0;i<dofs_per_cell;++i)
rho(local_dof_indices[i]) += cell_rhs(i); }

Secondly, I would like to calculate the average rho, so I do:

double integral=0.0;

for (; cell_p!=endc; ++cell_p)
   {
   fe_values.reinit (cell_p);
   fe_face_values.get_function_values (rho,values);

         for (unsigned int q_point=0;q_point<n_q_points;q_point++)
                integral += values[q] * fe_face_values.JxW(q);
   }

Are both pieces of code all right?

Does the value of integral correspond to {(1/Volume)*integral (rho)} or only to {integral(rho)}?

Thanks in advance.
Best
Isa

_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii

Reply via email to