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