Hi Andrew,
thank-you very much for your reply. Sorry for bothering you again.
I have thought about the problem and I have read more carefully through
the mailing list and I think I have found an easier way to do it.
Instead using "QGauss" quadrature formula for "FEValues", I can use
"Quadrature" along with "support_points" because in that way the order
of points in the array matches that returned by the
|cell->get_dof_indices| function, and so, I can obtain the 2nd
derivatives on dofs.
I have two more questions.
I would do:
/Quadrature<DIMENSION>
quadrature_formula(fe_fc.get_unit_support_points());
unsigned int dofs_per_cell = fe_fc.dofs_per_cell;
Vector<double> cell_rhs;
std::vector<unsigned int> local_dof_indices;
const unsigned int n_q_points = quadrature_formula.n_quadrature_points;/
/FEValues<DIMENSION> fe_values (fe_fc, quadrature_formula,
UpdateFlags(update_values |
update_gradients |
update_q_points
|
update_second_derivatives|
update_JxW_values));/
/dofs_per_cell = fe_fc.dofs_per_cell;
local_dof_indices.resize (dofs_per_cell);
cell_rhs.reinit (dofs_per_cell);
std::vector< Tensor< 2, DIMENSION > > values;
values.resize(n_q_points);
DoFHandler<DIMENSION>::active_cell_iterator cell_p =
dof_handler_total.begin_active(),
endc = dof_handler_total.end();
solution=0.0;
/
/for (; cell_p!=endc; ++cell_p)
{
cell_rhs=0.0;
fe_values.reinit (cell_p);
fe_values.get_function_2nd_derivatives (vector_potential, values);
for (unsigned int q_point=0;q_point<n_q_points;q_point++)
cell_rhs(q_point)-=cte_dielectrica*
(values[q_point][0][0]+values[q_point][1][1]+values[q_point][2][2]);
cell_p-> get_dof_indices(local_dof_indices);
for (unsigned int i=0;i<dofs_per_cell;++i)
solution(local_dof_indices[i]) += cell_rhs(i);
}
/
Am I right?
Later, I want to integrate the vector "solution", so I do:
/cell_p = dof_handler_total.begin_active();
endc = dof_handler_total.end();
std::vector<double > values2;
values2.resize(n_q_points);
double integral= 0.0;
for (; cell_p!=endc; ++cell_p)
{
fe_values.reinit (cell_p);
fe_values.get_function_values (solution,values2);
for (unsigned int q=0; q<n_q_points; ++q)
{
integral += values2[q]* fe_values.JxW(q);
}
}/
Here, the problem I have is that "fe_values.JxW(q)" has an "Inf" value
for all "q". However, if I use "QGauss" quadrature formula for the
"FEValues" (for calculating this integral), "fe_values.JxW(q)" has the
correct value.
Cannot use "fe_values.JxW(q)" if I have used "Quadrature" quadrature
formula applied to "FEValues"?
Thanks again!
My best
Isa
Andrew McBride escribió:
Hi Isabel
The function compute_projection_from_quadrature_points_matrix provides
the matrix X that transforms your scalar data defined at the
quadrature points to the nodes. The projection is constructed based on
the criteria described in the documentation wherein product of the
integral of the projection with arbitrary (discrete) w = the product
of the integral of the function u with w (so they are equivalent in
L^2 sense).
i.e. (V_h,w) = (u,w) \forall w (*)
The problem is we don't know u we only know the values at the
quadrature points.
Thats OK though. We can integrate the LHS and RHS of (*) using a
quadrature rule. The rhs_quadrature would then, sensibly, correspond
to the quadrature rule you are using to integrate your system matrix.
The LHS quadrature rule, in your case, should be the same as the RHS.
You would only need different rules if you were projecting, say, a
pressure determined using a single point quadrature rule to the nodes
of your cell (i.e. there is a mismatch in the order of the projection).
This function gives you the projection matrix X.
The next thing to understand is that that the Laplacian projected from
the quadrature points on an element to the nodes will be globally
discontinuous (generally).
There well may be a way to do this more efficiently in deal.II but
this way would work. You could for example simply project the Hessian
of \phi (your scalar variable) to the midpoint of the cell and then
determine the trace. This way you have one Laplacian per cell and you
don't need to project back to the nodes.
Regards
Andrew
Am 17 May 2010 um 5:12 PM schrieb Isabel Gil:
Hi everybody:
In a previous email about calculating the Laplacian of a scalar field named
"phi", you have told me to use:
/
/
/FETools::compute_projection_from_quadrature_points_matrix / /(/
/const FiniteElement
<http://www.dealii.org/6.1.0/doxygen/deal.II/classFiniteElement.html><
dim > & / //fe/, /
/
/ /
/ /const Quadrature
<http://www.dealii.org/6.1.0/doxygen/deal.II/classQuadrature.html><
dim > & / //lhs_quadrature/, /
/
/ /
/ /const Quadrature
<http://www.dealii.org/6.1.0/doxygen/deal.II/classQuadrature.html><
dim > & / //rhs_quadrature/, /
/
/ /
/ /FullMatrix
<http://www.dealii.org/6.1.0/doxygen/deal.II/classFullMatrix.html><
double <http://www.dealii.org/6.1.0/doxygen/deal.II/classdouble.html>
> & / //X// / /
/
/ /)/ /
/
So I want to calculate "rho(x,y,z)=-div(grad(phi))".
And I have:
/QGauss<DIMENSION> quadrature_formula(3);
const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
FEValues<DIMENSION> fe_values (fe_fc, quadrature_formula,
UpdateFlags(update_values |
update_gradients |
update_q_points |
update_normal_vectors |
update_second_derivatives|
//update_hessians|
update_JxW_values));
std::vector< Tensor< 2, DIMENSION > > values;
values.resize(n_q_points);
std::vector<double > values2;
values2.resize(n_q_points);/
/.
.
.
/
/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 q_point=0;q_point<n_q_points;q_point++)
values2[q_point]=(values[q_point][0][0]+values[q_point][1][1]+values[q_point][2][2]);
//.
.
.
}
/
It is clear that "fe_fc" is the variable "fe" in method
FETools::compute_projection_from_quadrature_points_matrix()
But, how should I take the variables "lhs_quadrature",
"rhs_quadrature" and the matrix "X". All of them are the variables
inside the previous method.
Or perhaps there is an easier way to calculate the Laplacian of a known scalar field.
Thanks in advace.
My best.
Isa
Andrew McBride escribió:
Sorry, I meant the trace of the Hessian.
You only need the Jacobian determinant if you performing an
integration and need to transform the limits of integration. The
value of the interpolation function at the quadrature points on the
reference cell and the actual cell are the same.
You would need to divide by the volume of the cell to get the
average value on the cell.
Andrew
Am 17 May 2010 um 11:56 AM schrieb Isabel Gil:
Hi!
Thanks for your answer.
Andrew McBride escribió:
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.
What I want to calculate is the Laplacian of Phi over the nodes of
my mesh.
When I use the loop "for (unsigned int
q_point=0;q_point<n_q_points;q_point++)", I always multiply by
"fe_values.JxW(q_point)" to get the transformation from reference
to real cell.
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?
Sorry, I made a mistake :-[ . What I wanted to use (and indeed I
use) was "fe_values.get_function_values (rho,values);"
But, once I have the "integral" value, should I divide it by the
volume to get the average or does the "integral" value itself
contain the avarage?
Thanks again.
Best.
Isa
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)
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
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
