I highly appreciate all of your answers. It becomes clearly for me to look
through the library.
Let me show all of my function and my ideas below.
I set the point for evaluating its solution to be on the vertex of a cell,
and using the following function to get the solution
template<int dim>
void OutputTools<dim>::point_solution(const Mapping<dim> &mapping,
//input mapping
const DoFHandler<dim>
&dof_handler, //input dof_handler
const Vector<double> &solution,
//input solution vector
const Point<dim> &point,
//input point in real coordinates to get the solution
Vector<double> &point_solution)
const { //output
const std::pair<typename DoFHandler<dim>::active_cell_iterator, Point<dim>>
cell_point =
GridTools::find_active_cell_around_point(mapping, dof_handler, point);
//find active cell around input point
const double delta = 1e-20;
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < delta,
ExcInternalError());
const Quadrature<dim>
quadrature(GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
//get quadrature by project the point to unit cell
FEValues<dim> fe_values(mapping, dof_handler.get_fe(), quadrature,
update_values | update_gradients |
update_quadrature_points);
fe_values.reinit(cell_point.first);
const FEValuesExtractors::Vector vector_re(0);
const FEValuesExtractors::Vector vector_im(dim);
std::vector<Tensor<1, dim>> temp_solution_re(quadrature.size());
std::vector<Tensor<1, dim>> temp_solution_im(quadrature.size());
//Get e_field at point
fe_values[vector_re].get_function_values(solution, temp_solution_re);
fe_values[vector_im].get_function_values(solution, temp_solution_im);
for (unsigned int j = 0; j < dim; ++j) {
point_solution(j) = temp_solution_re[0][j];
point_solution(j + dim) = temp_solution_im[0][j];
}
//using point_value() for comparing the result
Vector<double> values (6);
VectorTools::point_value(dof_handler,solution,point,values);
//print out all the results for testing the code
std::cout<<"\n";
std::cout<<"cell: \t"<<cell_point.first<<"\n";
std::cout<<"quadrature size: \t"<<quadrature.size()<<"\n";
std::cout<<"measuring point: \t"<<fe_values.get_quadrature_points()[0]<<"\n";
std::cout<<"size of temp_solution_re: "<<temp_solution_re.size() <<"\n";
std::cout<<"result with get_function_values real: \t"<<
temp_solution_re[0] << ", "<<temp_solution_re[1] << ",
"<<temp_solution_re[2] <<"\n";
std::cout<<"result with get_function_values imag: \t"<<
temp_solution_im[0] << ", "<<temp_solution_im[1] << ",
"<<temp_solution_im[2] <<"\n";
std::cout<<"result with point_value(): \t"<<values<<"\n";
}
The results printed by this function like this:
cell: 0.1141
quadrature size: 1
measuring point: -2000 0 -500
size of temp_solution_re: 1
result with get_function_values real: 3.03537e-09 5.77993e-11
-6.4355e-09, 1.11165e-321 4.68544e-310 0.1, 0.1 0.443649 0.443649
result with get_function_values imag: -1.51054e-09 -1.09549e-11
8.0214e-10, 1.63042e-322 3.03537e-09 5.77993e-11, -6.4355e-09
1.11165e-321 4.68544e-310
result with point_value(): 3.035e-09 5.780e-11 -6.436e-09 -1.511e-09
-1.095e-11 8.021e-10
1. When evaluating point is on the vertex of a cell, this means that there
should be several cells in 3D have this vertex.
Thus, my first question is the function :
const std::pair<typename DoFHandler<dim>::active_cell_iterator, Point<dim>>
cell_point =
GridTools::find_active_cell_around_point(mapping, dof_handler, point);
I personally think that as I set it to be const, so the return cell
(cell_point.first) is only one cell, rather than a group of cells.
Could you please tell me about this?
2. My second idea is with the codes:
const Quadrature<dim>
quadrature(GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
FEValues<dim> fe_values(mapping, dof_handler.get_fe(), quadrature,
update_values | update_gradients |
update_quadrature_points);
fe_values.reinit(cell_point.first);
I have only one quadrature defined on the cell, and the printing
result show me the same thing.
Can I use this way to get the solution, or it would be more efficient
when using codes like:
QGaussLobatto<dim> quadrature_formula(2);
FEValues<dim> fe_values(fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
fe_values.reinit(cell_point.first);
And then finding the quadrature point coincide with the evaluating
point to get the solution?
3. My third thought is in my function I have only 1 quadrature point
so that temp_solution_re[0] and temp_solution_im[0] are the variables
containing the solution.
The printing results showed that the solution from point_value() is
identical with those from temp_solution_re[0] and temp_solution_im[0].
Do temp_solution_re[1], temp_solution_re[2],
temp_solution_im[1],temp_solution_im[2] make any sense for my
solution? Because I saw the code can run with these variables.
Thank you very much!
Vào Th 6, 15 thg 3, 2019 vào lúc 00:29 Wolfgang Bangerth <
bange...@colostate.edu> đã viết:
> On 3/14/19 1:28 AM, Phạm Ngọc Kiên wrote:
> >
> > So for the first function
> > VectorTools::point_value(dof_handler,solution,point,values);
> > I can specify the observed point in real coordiante, solution vector and
> > dof_handler to get the values, right?
>
> Correct.
>
>
> > For the second one, assuming that I am using QGauss for solving the
> > system, does it mean that if the point is on the edge of the cell then I
> > cannot get the solution at that point?
>
> Yes. It can be any quadrature formula defined on the reference cell. It
> doesn't matter where the points are located.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
