Hi Isabel
Creating a quadrature rule based on the the cell support points is a neat way
to have data extrapolated there. The problem is that by creating the
quadrature_formula in this way the quadrature weighting is unspecified
(actually, its set to inf). That is why you can't use the same quadrature rule
for integration.
So you have the option of creating your own quadrature rule with weights and
points defined at the support points of the cell. Quadrature of this form is
used to create diagonal mass matrices (see, e.g. The Finite Element Method,
Hughes).
Otherwise some of the methods that have been suggested before in this post will
definitely work.
My pragmatic suggestion would be to project to the midpoint of the cell using a
QMidPoint quadrature rule and take this as the average value (i.e. no need to
integrate to get the average). Then project the scalar from the cell support
points to the cell support points to obtain the Hessian using a different
quadrature rule based on the cell support points and take this is the value at
the nodes.
As we have discussed though, there are many ways to do this it appears. Perhaps
someone else has a better idea but all these approaches will work.
Cheers
Andrew
Am 18 May 2010 um 6:32 PM schrieb Isabel Gil:
> 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< dim > & fe,
>>>
>>>
>>> const Quadrature< dim > & lhs_quadrature,
>>>
>>>
>>> const Quadrature< dim > & rhs_quadrature,
>>>
>>>
>>> FullMatrix< double > & 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
