On 3/16/19 7:16 AM, jane....@jandj-ltd.com 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: bange...@colostate.edu
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.
