One final remark for today, if I remove the two asymmetrical terms in my
equation so it becomes :
*Then all three codes gives back the same solution...*
So I thought it would come from my assembly of this part of the system but
I have checked and checked again... It should not come from my weak
formulation since the fenics and dealii have the same one.
Something to do with the matrix no longer being symmetrical ? But I used
UMFPack as a solver so it should not be an issue...
Here is my assembly process for the stokes problem in case you see
something wrong
> void Elfe::stokes_assemble_system ()
> { //Fonction qui assemble le problème de Stokes
>
> std::cout << BLUE << "Assemblage Stokes " << RESET
> << std::endl;
>
> //On réinitialise la matrice et le rhs
> stokes_system_matrix = 0;
> stokes_system_rhs = 0;
>
> //On charge de quoi extraire les valeurs
> QGauss<2> quadrature_formula(4);
> const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell;
> const unsigned int n_q_points = quadrature_formula.size();
>
> std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell
> );
>
> FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
> Vector<double> local_rhs (dofs_per_cell);
>
> //Stokes values
> FEValues<2> stokes_fe_values(stokes_fe, quadrature_formula,
> update_values | update_gradients |
> update_JxW_values);
>
> double u_i, v_i;
> double u_i_x, u_i_y, v_i_x, v_i_y;
> double u_j, v_j;
> double u_j_x, u_j_y, v_j_x, v_j_y;
> double P_i, P_j;
>
> const FEValuesExtractors::Vector Velocities(0);
> const FEValuesExtractors::Scalar Pressure(2);
>
> //Phi values
> FEValues<2> phi_fe_values(phi_fe, quadrature_formula,
> update_values | update_gradients |
> update_JxW_values);
>
> double phi, phi_x, phi_y;
> double DX;
>
> double philin = 2*pi*n;
> double phiscale = 2*pi*s;
>
> std::vector<Tensor<1, 2> > phi_solution_gradients(n_q_points);
> std::vector<double> phi_solution_values(n_q_points);
>
> //Il nous faut 2 cell iterator, un pour phi, un pour stokes
> typename DoFHandler<2>::active_cell_iterator
> stokes_cell = stokes_dof_handler.begin_active(),
> stokes_endc = stokes_dof_handler.end();
>
> typename DoFHandler<2>::active_cell_iterator
> phi_cell = phi_dof_handler.begin_active();
>
> for (; stokes_cell!=stokes_endc; ++stokes_cell, ++phi_cell)
> { //Pour chaque cell
> stokes_fe_values.reinit(stokes_cell);
> phi_fe_values.reinit (phi_cell);
>
> local_matrix = 0;
> local_rhs = 0;
>
> //On récupère les valeurs pour les angles phis
> phi_fe_values.get_function_values(phi_solution,
> phi_solution_values);
> phi_fe_values.get_function_gradients(phi_solution,
> phi_solution_gradients);
>
> for (unsigned int q=0; q<n_q_points; ++q)
> { //Pour chaque point de quadrature
> DX = stokes_fe_values.JxW(q);
>
> phi = phi_solution_values[q];
> phi_x = phi_solution_gradients[q][0];
> phi_y = phi_solution_gradients[q][1];
>
> for (unsigned int i=0; i<dofs_per_cell; ++i)
> {
> u_i = stokes_fe_values[Velocities].value (i, q)[0];
> v_i = stokes_fe_values[Velocities].value (i, q)[1];
> u_i_x = stokes_fe_values[Velocities].gradient(i, q)[0][0];
> u_i_y = stokes_fe_values[Velocities].gradient(i, q)[0][1];
> v_i_x = stokes_fe_values[Velocities].gradient(i, q)[1][0];
> v_i_y = stokes_fe_values[Velocities].gradient(i, q)[1][1];
> P_i = stokes_fe_values[Pressure] .value (i, q);
>
> local_rhs(i) += (
>
> u_i*(
> 8*philin*phiscale*X
> *(cos(2*phiscale*phi)*phi_x + sin(2*phiscale*phi)*
> phi_y)
>
> + 2*phiscale*phiscale*phi_x*(-4*philin/phiscale)
> )
>
> + v_i*(
> 8*philin*phiscale*X
> *(sin(2*phiscale*phi)*phi_x - cos(2*phiscale*phi)*
> phi_y)
>
> + 2*phiscale*phiscale*phi_y*(-4*philin/phiscale)
> )
>
> )*DX;
>
> for (unsigned int j=0; j<dofs_per_cell; ++j)
> {
> u_j = stokes_fe_values[Velocities].value (j, q)[0
> ];
> v_j = stokes_fe_values[Velocities].value (j, q)[1
> ];
> u_j_x = stokes_fe_values[Velocities].gradient(j, q)[0
> ][0];
> u_j_y = stokes_fe_values[Velocities].gradient(j, q)[0
> ][1];
> v_j_x = stokes_fe_values[Velocities].gradient(j, q)[1
> ][0];
> v_j_y = stokes_fe_values[Velocities].gradient(j, q)[1
> ][1];
> P_j = stokes_fe_values[Pressure] .value (j, q);
>
> local_matrix(i, j) += (
> //Classical Stokes terms
> - u_i_x*u_j_x - u_i_y*u_j_y - v_i_x*v_j_x - v_i_y*
> v_j_y
> + (u_i_x+v_i_y)*P_j
> + P_i*(u_j_x+v_j_y)
>
> //Leslie Asymetrical terms
> + phiscale/a
> *(v_i_x-u_i_y)
> *(u_j*phi_x + v_j*phi_y + 1/2/phiscale*(u_j_y-
> v_j_x))
>
> //Terme provenant du laplacien phi
> -2*phiscale*phiscale/a
> *(u_i*phi_x + v_i*phi_y)
> *(u_j*phi_x + v_j*phi_y + 1/2/phiscale*(u_j_y-
> v_j_x))
>
> )*DX;
> }
> }
> }
>
> stokes_cell->get_dof_indices (local_dof_indices);
>
> //On applique les boundary conditions
> stokes_constraints.distribute_local_to_global(local_matrix,
> local_rhs,
> local_dof_indices,
> stokes_system_matrix,
> stokes_system_rhs);
> }
>
> }
>
>
>
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/a984ed95-fa81-4758-bc46-1357ec5ce2b5%40googlegroups.com.
