Hi Wolfang, thanks for your time on this.
Ok, re the first step-20 bc issue, I'll have another think, but I still am
not sure why then it isn't giving me the exact figure, whilst my suggestion
is (I understand your point here and I would have said I agreed with you,
but my implementation does work). I've verified this part of the code using
exact solutions so I really have no idea what the issue is, but I guess
i'll have to have another think.
On the second, step-22 issue, i don't have u.n=0 conditions. I have partial
boundary conditions in the form of that stated in the tutorial (1st in the
fourth set of Bc explanations) of:
u_t = 0
n.(n.(pI-2*epsilon)) = nonzero value.
So this is imposed weakly, and no strong conditions are set on that
boundary at all, and that nonzero value equals the topstress_value I had
before in the weak form.
On Saturday, March 16, 2019 at 3:53:10 PM UTC, Wolfgang Bangerth wrote:
>
> On 3/16/19 7:16 AM, jane...@jandj-ltd.com <javascript:> wrote:
> > So that's what I meant by having extra contributions. I thought you
> needed a
> > local_rhs = 0 after/within:
> > for (unsigned int face_n = 0;
> > face_n < GeometryInfo<dim>::faces_per_cell;
> > ++face_n)
> > if (cell->at_boundary(face_n))
> > {
> > fe_face_values.reinit(cell, face_n);
> >
> > In the same way that if you enforce them strongly, it basically ignores
> the
> > contributions to the cell vertices that are on the boundary. Because
> this is a
> > Dirichlet condition, should we not reset the local_rhs contribution that
> have
> > accumulated just prior to the boundary condition? putting local_rhs = 0
> again
> > within the boundary contribution terms seems to actually do the job.
>
> You need to distinguish between Dirichlet/Neumann boundary conditions, and
> strong/weak boundary conditions. Weak boundary conditions are the ones
> that go
> into the weak formulation; for the typical Laplace equation, Neumann
> boundary
> conditions are weak and Dirichlet boundary conditions are strong. But for
> the
> mixed form you are using in step-20, it is the other way around.
>
> As for whether you should reset local_rhs=0: Think about this from a data
> flow
> perspective. In the block you quote above, you are computing something for
> each face. You are *writing* it into local_rhs. But what are you *doing*
> with
> what you have computed? If you set local_rhs=0 between each face you
> visit,
> then you are simply throwing away what you have computed on previous
> faces,
> before you use it for something (namely, putting it into the global
> matrix).
>
> So what the code is doing is to first accumulate the contributions of all
> (boundary) faces of a cell before it puts it into the global matrix. Then
> it
> moves to the next cell, sets local_rhs=0, accumulates contributions of all
> (boundary) faces of that cell, puts it into the global matrix, and so on.
>
> If you reset local_rhs=0 between each two faces, then you are computing
> the
> contribution of each (boundary) face alright, but because you overwrite it
> with zeros when you move to the next face, the only thing that ends up in
> the
> global matrix are only the contributions of the last (boundary) face you
> visit
> on each cell.
>
> If this makes a difference, you may ask why that is the case and when it
> actually matters. Specifically, if you put local_rhs=0 right before or
> after
> the call to fe_face_values.reinit(...) above, then you will be zeroing out
> any
> computation you may have done on the previous boundary faces of the
> current
> cell -- but this is only going to make a difference on cells that actually
> have two or more boundary faces, i.e., the ones in the corners of your
> domain.
> I don't know why that would matter in your case, but it's an observation
> that
> may help you think about what could be wrong in your program.
>
>
> > For Stokes - normal component of the normal stress:
> > Ok, that's interesting. I seems to be getting a non-zero contribution
> when I do:
> >
> > local_rhs(i) = topstress_values[q] * fe_face_values.normal_vector(q)*
> >
> fe_face_values[velocities].value(i,q_point)*
> > fe_face_values_rock.JxW(q);
> >
> > where topstress_values is the value of the normal component of the
> normal stress.
>
> Yes, this may be nonzero, but when you call
> constraints.distribute_local_to_global(), it zeros out rows and columns
> that
> correspond to degrees of freedom with strong boundary conditions. So, for
> example, if topstress_values*normal_vector has a nonzero normal component,
> then you will get something nonzero here, but you have prescribed strong
> boundary conditions of the form u.n=0, then the nonzeroness here will be
> zeroed out by the constraints.
>
> 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.
