Hello everybody!
I would like to calculate the volume of flow (m³/s) at the exit boundary
and I do it as it is pointed out below:
Once I have the velocity vector ("vector_solution"), I do:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
QGauss<DIMENSION-1> face_quadrature_formula(DIMENSION);
unsigned int dofs_per_cell = fe_total.dofs_per_cell;
//* Where FE_Q<DIMENSION> fe_velocity(2),
FE_Q<DIMENSION> fe_pressure(1);
FESystem<DIMENSION> fe_total(fe_velocity,DIMENSION, fe_pressure,1)*/
/const unsigned int n_face_q_points =
face_quadrature_formula.n_quadrature_points;
FEFaceValues<DIMENSION> fe_face_values (fe_total, face_quadrature_formula,
UpdateFlags(update_values |
update_gradients |
update_q_points |
update_normal_vectors |
update_second_derivatives|
update_JxW_values));
dofs_per_cell = fe_total.dofs_per_cell;
local_dof_indices.resize (dofs_per_cell);
std::vector<Vector<double> > values;
values.resize(n_face_q_points);
for (unsigned int i=0;i<n_face_q_points;i++)
values[i].reinit(DIMENSION+1);
DoFHandler<DIMENSION>::active_cell_iterator cell_p =
dof_handler_total.begin_active(),
endc = dof_handler_total.end();
double valueOUT=0.0;
for (; cell_p!=endc; ++cell_p)
{
for (unsigned int face=0;
face<GeometryInfo<DIMENSION>::faces_per_cell; ++face)
{
if ( (cell_p->face(face)->boundary_indicator() == out_flag)
)
{
fe_face_values.reinit (cell_p, face);
fe_face_values.get_function_values(vector_solution,
values); ///Where I have the velocity solution vector /
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int
comp_i=fe_total.system_to_component_index(i).first;
for (unsigned int q_point=0; q_point<n_face_q_points;
++q_point)
{
if ((comp_i==0))
{
valueOUT += ((
((fe_face_values.normal_vector(q_point)(0))*
(values[q_point](0)))+
(
(fe_face_values.normal_vector(q_point)(1))*
(values[q_point](1)) )+
(
(fe_face_values.normal_vector(q_point)(2))*
(values[q_point](2)) ))*
fe_face_values.JxW(q_point) *
fe_face_values.shape_value(i,
q_point)
);
}
}
}
}
}
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
So, in "valueOUT" variable, I have the outflow (m³/s), don't I?
Thanks in advance
Best regards
Isa
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
