Hi DealII's users:
(I wrote a bizarre, not clearly enough, question some days ago in
relation to this subject). Reading the previous emails, I think I can
delimit my problem.
I would like to interpolate a solution from one mesh to another, both
are unstructured meshes (generated externally) and I have seen in the
mailing list that the answer to this problem was:
*****************************************************************************************
Unfortunately, there is no simple way of doing it (this is not a lack of
deal.II, but a mathematical statement). In particular, interpolation
between general meshes is of order h only and spoils the accuracy of most
numerical methods. L^2-projection would be better.
In both cases, you'd have to evaluate a function on one mesh at the
support or quadrature points of the other. Therefore, you first have to
compute this point (using FEValues) and then find the cell containing it
in the other mesh, possibly compute its transformation to the
corresponding unit cell and compute the function values at that point. All
this will cause an incredible amount of time on fine meshes.
*****************************************************************************************
I have:
triangulation1;
dof_handler_1(triangulation1);
FE_Q<DIMENSION> fe_v_1 (2);
FE_Q<DIMENSION> fe_p_1(1);
FESystem<DIMENSION> fe_total_1((fe_v_1,DIMENSION, fe_p_1,1);
dof_handler_1.distribute_dofs(fe_total_1);
DoFRenumbering::component_wise(dof_handler_1);
vector_solution1.reinit(dof_handler_1.n_dofs());
Once I have the solution in this vector_solution1 I would like to
interpolate it into another mesh:
triangulation2;
dof_handler_2(triangulation2);
FE_Q<DIMENSION> fe_v_2 (2);
FE_Q<DIMENSION> fe_p_2(1);
FESystem<DIMENSION> fe_total_1((fe_v_2,DIMENSION, fe_p_2,1);
dof_handler_2.distribute_dofs(fe_total_2);
DoFRenumbering::component_wise(dof_handler_2);
vector_solution2.reinit(dof_handler_2.n_dofs());
So as I would like to have the best accuracy as possible, I should use
L^2-projection, shouldn't I? and you advised: you'd have to evaluate a
function on one mesh at the support or quadrature points of the other.
Therefore, you first have to compute this point (using FEValues) and
then find the cell containing it in the other mesh, possibly compute its
transformation to the corresponding unit cell and compute the function
values at that point.
Sorry, but I don't get to understand what exactly I have to do. First, I
define:
//Quadrature formule:
QGauss3<DIMENSION> quadrature_formula;
const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
//I create two FEValues objects:
FEValues<DIMENSION> fe_values_1(fe_total_1, quadrature_formula,
UpdateFlags(update_values |update_gradients |update_q_points
|update_JxW_values));
FEValues<DIMENSION> fe_values_2(fe_total_2, quadrature_formula,
UpdateFlags(update_values |update_gradients |update_q_points
|update_JxW_values));
//The number of dofs per cell will be the same, so
dofs_per_cell_1=dofs_per_cell_2
const unsigned int dofs_per_cell_1 = fe_total_1.dofs_per_cell;
const unsigned int dofs_per_cell_2 = fe_total_2.dofs_per_cell;
FullMatrix<double> cell_matrix_1(dofs_per_cell_1,dofs_per_cell_1);
Vector<double> cell_rhs_1(dofs_per_cell_1);
FullMatrix<double> cell_matrix_2(dofs_per_cell_2,dofs_per_cell_2);
Vector<double> cell_rhs_2(dofs_per_cell_2);
std::vector<unsigned int> local_dof_indices_1 (dofs_per_cell_1);
std::vector<unsigned int> local_dof_indices_2 (dofs_per_cell_2);
// I create two cell accesors because I have two dof_handler objects:
DoFHandler<DIMENSION>::active_cell_iterator
cell_1= dof_handler_1.begin_active(),
endc_1 = dof_handler_1.end();
DoFHandler<DIMENSION>::active_cell_iterator
cell_2= dof_handler_2.begin_active(),
endc_2 = dof_handler_2.end();
// I create two cell accesors because I have two triangulations objects:
Triangulation<DIMENSION>::cell_iterator celltria_1 =
triangulation1.begin (),
endctria_1 =
triangulation1.end();
Triangulation<DIMENSION>::cell_iterator celltria_2 =
triangulation2.begin (),
endctria_2 =
triangulation2.end();
std::vector<Vector<double> > values_1;
values_1.resize(n_q_points);
std::vector<Vector<double> > values_2;
values_2.resize(n_q_points);
for (unsigned int i=0;i<n_q_points;i++)
{
values_1[i].reinit(DIMENSION+1);
values_2[i].reinit(DIMENSION+1);
}
for (; cell_1!=endc_1; ++cell_1)
{
cell_matrix_1=0.0;
cell_rhs_1=0.0;
fe_values_1.reinit (cell_1);
fe_values_2.reinit (cell_2);
fe_values_1.get_function_values(vector_solution_1, values_1);
// HERE I SHOULD INITIALIZE THE PROCESS L² PROJECTION FROM ONE MESH TO
THE ANOTHER, BUT I AM REALLY MISSING.
}
Could you show me the scheme how realize this process?
Thanks in advance.
Best
Isa
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
