Hi Wolfgang, thanks so much for getting back to me
1. I was wondering re dimensions because I could find a component mask
function or something similar to fe_values[velocities], eg, when using
fe_face_values which you need to apply the neumann conditions.
2. This is my fault. I meant working in 2d with a solution that only acts
vertically, hence choosing a solution that is only dependent on y and that
is 0 in the xdirection.
3. indeed i was showing the full system right hand side component
4. for the top and bottom boundaries, you are right, i do mean the stress
as mentioned in previous messages. Either way, I am now implementing things
in the neumann BC way as in step-22. with the solution as before, i add the
following to my top boundary:
for (unsigned int face_number=0;
face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)
if (cell->face(face_number)->at_boundary()
&&
(cell->face(face_number)->boundary_id() == 1))
{
fe_face_values.reinit (cell, face_number);
topstress.vector_value_list(fe_face_values.get_quadrature_points(),
topstress_values);
for (unsigned int q_point=0; q_point<n_face_q_points;
++q_point)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
local_rhs(i) +=
(topstress_values[q_point](component_i)*
fe_face_values.
shape_value(i,q_point) *
fe_face_values.JxW(q_point));
}
}
}
where the top stress is g_N as in step-22 (where i calculate the normal as
plus or minus (0,1) depending on whether it's the top or bottom boundary).
I am running this off exactly using the equatoins in step-22 with the right
hand sides modified and using a direct solver instead.
The thing i have worked out is that as my program runs currently, but with
DIRICHLET conditions on top and bottom, it works. but when i try to use the
nuemann conditions as i have shown in the code snippet above does not (with
everything else kept the same). So i am wondering whether the snippet i
have shown here contains an error I have not figured out.
Thanks for your help
On Monday, February 19, 2018 at 5:44:39 PM UTC-5, Wolfgang Bangerth wrote:
>
>
> Jane,
>
> > Firstly, I decided not to use the normal vector for now. Since the
> > normal vector is 2D, i wasn't sure how to implement the rest so that it
> > is a 'double' since my g (neumann condition vector) has 3 components
> > when the normal vector will only have 2?
>
> I'm not sure I understand this -- in 2 space dimensions, the normal
> vector has 2 components; in 3 space dimensions, it has 3 components.
>
>
> > I went back to the beginning - step-22 tutorial and started from
> > absolute basics. Testing with 1D, i have started with exact solution
> > u=(0,-y^2)^T and p = y^3, ie purely vertical.
>
> I also don't understand this, I'm afraid: in d space dimensions, the
> velocity vector has d components. So if you're in 1d, how come your
> velocity vector has 2 components?
>
> Similarly...
>
> > Doing a manufactured
> > solution type, i get the right hand side vector f = (0, 4+3z^2, 2z)^T.
>
> the right hand side has d components. I suspect you are here showing the
> right hand side of the whole system, which should then have d+1
> components (which in 1d should be 2, not 3).
>
> Or do you want to suggest in all of this that you really have a
> two-dimensional domain but that the flow is one-dimensional?
>
>
> > with my solution, on a rectangular domain, i have no tangential stresses
> > anywhere, no flux on sides (unit normal being (/pm 1,0)), and plus and
> > minus (0, z^3 + 4z) at the top.
>
> Your domain is [0,1]^2, so I understand how there is no flux through the
> bottom because there u=[0,0]. But how do you arrive at the flux you show
> here? The velocity field is quadratic in z (y?), and so should the flux,
> but you show something that is cubic with the distance.
>
> If instead of 'flux' you really mean 'stress', then it should actually
> only be linear in z since it involves the *derivative* of the velocity.
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> <javascript:>
> 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.
