Hi Isabel
I'm confused by a couple of things.
You want to determine the trace of the Laplacian of Phi at the quadrature
points and then extrapolate the calculated values back to the nodes? Why are
you multiplying the value at the quadrature point by fe_values.JxW(q_point)?
This implies that you are doing integration based on nodal values as opposed to
extrapolation. What you want to do is just extrapolate the quadrature point
data back to the nodes without integrating. The function you need to do this
FETools::compute_projection_from_quadrature_points_matrix
The projection is on an element so the result will be globally discontinuous.
In calculating the integral why are you using the
fe_face_values.get_function_values (rho,values)? I would extrapolate the values
at the nodes to the quadrature points and then integrate them using a sum over
the quadrature points and multiplying each quadrature point contribution by
fe_values.JxW(q_point). There may even be a deal function that does this?
Cheers
Andrew
Am 17 May 2010 um 10:58 AM schrieb Isabel Gil:
> 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
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
