Hi everybody
I am a new deal.II user and I am trying to make my way through the
distribution.I am trying to set the normal component of the Traction
vector for Stokes flow however I am doing something (or plenty) wrong
Which i cannot figure out here are some relevant snippets of the code
template <int dim>
class Traction:public Function<dim>
{
public:
Traction():Function<dim>(dim+1){}
virtual void vector_value(const Point<dim>& p,Vector<double>&
values)const;
virtual void vector_value_list(const std::vector<Point<dim> >&
points,std::vector<Vector<double> >& value_list)const;
};
template<int dim>
void Traction<dim>::vector_value(const Point<dim>& p,Vector<double>&
values)const {
Assert(values.size()==dim,ExcDimensionMismatch(values.size(),dim));
values(0)=0.02;
values(1)=0.0;
values(2)=0.0;
}
template<int dim>
void Traction<dim>::vector_value_list( const std::vector<Point<dim>
>&points,std::vector<Vector<double> >& value_list)const { unsigned int
n_points=points.size(); for(unsigned int p=0;p<n_points;++p) {
Traction<dim>::vector_value(points[p],value_list[p]);
}
}
Then I assemble the right hand side as
for(unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
{
if(cell->face(face)->at_boundary() &&
(cell->face(face)->boundary_indicator()==1))
{
fe_face_values.reinit(cell,face);
//applied_traction.vector_value_list(fe_face_values.get_quadrature_point
s(),traction_values);
//perform the integration on the contour integral by using quadrature
for(unsigned int q_point=0; q_point<n_face_q_points; ++q_point) {
//const double
traction=applied_traction.vector_value_list(fe_face_values.quadrature_po
int(q_point));//*fe_face_values.normal_vector(q_point);
for(unsigned int i=0; i<dofs_per_cell; ++i) { const unsigned int
component_i =
fe.system_to_component_index(i).first;
//if(component_i==0)
local_rhs(i) += fe_face_values.shape_value(i,q_point) *
traction_values[q_point] (component_i)*
fe_face_values.JxW(q_point);
//local_rhs(i)
+=(traction*fe_face_values.shape_value(i,q_point))*fe_face_values.JxW(q_
point);
}
}
}
}
The rest boundary conditions are u1=u2=0 on the upper and lower parts of
the channel and at the most left and right I set u2=0 and leave the u1
component unspecified by doing the following
std::vector<bool> component_mask_1 (dim+1,true);
//exclude u1 velocity and pressure
component_mask_1[0]=false;
component_mask_1[dim]=false;
Apart from these changes the rest is the same as in step-22 tutorial.The
results that I got are somehow strange When I set the normal component
at the left boundary the pressure distribution does not look right
however when I move to the right boundary the pressure decays linearly
while the velocity u1 is, qualitatively, correct in absolute values.i am
sure that I am missing something here but I cannot figure out what so
any comments are mostly welcome.
Thanks in advance Alex
This message has been checked for viruses but the contents of an attachment
may still contain software viruses, which could damage your computer system:
you are advised to perform your own checks. Email communications with the
University of Nottingham may be monitored as permitted by UK legislation.
_______________________________________________
