dolfin team mailing list archive

[Patrick Riesen <priesen@xxxxxxxxxxxxxxx>]

Johan Hoffman wrote:
> Patrick Riesen wrote:
> > Patrick Riesen wrote:
> > > Garth N. Wells wrote:
> > > > On 17/05/10 08:50, Patrick Riesen wrote:
> > > > > Anders Logg wrote:
> > > > > > On Fri, May 07, 2010 at 02:39:41PM +0200, Patrick Riesen wrote:
> > > > > > > hello,
> > > > > > >
> > > > > > > i want to zero the normal velocity on a boundary part of my
> > > > > > > 2D-domain, so i used something like
> > > > > > >
> > > > > > > n = cell.n
> > > > > > > a = ... + inner(w,dot(v,n)*n)*ds(2) + ...
> > > > > > >
> > > > > > > in my form-file (w=testfunction, v=trialfunction).
> > > > > > >
> > > > > > > i marked the boundary part using the "exterior facet domains"
> > > > > > > meshfunction, which will be picked up in assemble(..):
> > > > > > >
> > > > > > > MeshFunction< uint >* _ext_facet_domains =
> > > > > > > mesh.data().create_mesh_function("exterior facet domains", 1)
> > > > > > >
> > > > > > > *_ext_facet_domains = 5;
> > > > > > > boundary.mark(*_ext_facet_domains, 2)
> > > > > > >
> > > > > > > can i do this? or what is wrong with the above way? because i don't
> > > > > > > get zero velocities...
> > > > > > >
> > > > > > > thank your for help,
> > > > > > >
> > > > > > > patrick
> > > > > >
> > > > > > What if you add a (big) parameter in front of that term? Does it have
> > > > > > any effect at all on the solution?
> > > > > >
> > > > > > I would also suggest writing the term as
> > > > > >
> > > > > > dot(w, n)*dot(v, n)*ds(2)
> > > > >
> > > > > thank you for the help anders,
> > > > >
> > > > > so i noticed the following:
> > > > > without parameter in front, i could see no effect at all. if i add a
> > > > > parameter C in front and increase it stepwise until C is about 10^6 i
> > > > > can observe how the velocities cancel out on that boundary.
> > > >
> > > > You can probably improve this significantly by enforcing zero normal 
> > > > velocity using Nitsches's method. This involves just adding two 
> > > > extra terms to your form. See demo/pde/dg/poisson for an example.
> > > >
> > > > Garth
> >
> > hello,
> > so based on the literature from peter hansbo
> > http://www.math.chalmers.se/~hansbo/nitscherev_preprint.pdf
> >
> > i used the nitsche's method so impose a ALE boundary condition of the 
> > form
> >
> > (v - v_m)*n = 0
> >
> > on my (moving) boundary (garth: it ends up to 5 additional terms).
> > (v=fluid velocity, v_m = mesh velocity, n = outward normal).
> > I think it works correctly but i observe a loss of volume because the 
> > above bcs is only fullfilled to differences to about 1.0e-4 to 1.0e-5, 
> > so in every timestep i loose some mass/volume due to that, i guess.
> >
> > in hansbo's script there is not much about the stability parameter 
> > h^-1 needed and i just used the same as in the dg/poisson example, but 
> > does somebody know more about how to choose or estimate it?
> >
> > regards,
> > patrick
> If your boundary is curved (that is: not a line or plane), applying your 
> normal bc (slip bc) weakly in this way based on the default 
> discontinuous face/edge normals will result in artificial friction, or 
> if your penalty is large enough zero velocity.

yes, i noticed when my penalty coefficient is large, it doesn't help the 
solution, but i just get closer to v = v_m.

> What you can do is:
>
> (1) apply your bc strongly at each surface node instead, based on a  
> vertex normal, e.g. from a weighted average from surrounding face/edge 
> normals (this is also mass conservative), see:

this was where i first started before the weak imposition, but i had 
difficulties in creating those vertex normals.

> http://www.csc.kth.se/~jhoffman/archive/papers/scm.pdf
>
> (2) form a continuous normal field to apply your bc weakly, for example 
> a p.w. linear normal field over the surface. This may also work 
> depending on what method you use to weakly apply your bc.

so i would solve for a normal field? can you elaborate some more on this?

> Good luck,
> Johan
thank you (for the help) ;-)

regards,
patrick


> > > sounds interesting, i'll look into it.
> > > thank you,
> > > patrick
> > >
> > > > > i guess this approach works, but there is something else:
> > > > > looking at the (time-dependent) solution (*pvd file set or single vtu
> > > > > files) in paraview, i ran into several
> > > > >
> > > > > ERROR: In
> > > > > /home/kitware/Dashboard/MyTests/ParaView-3-8/ParaView-3.8/ParaView/VTK/IO/vtkXMLDataReader.cxx, 
> > > > >
> > > > > line 509
> > > > > vtkXMLUnstructuredGridReader (0x1398ce0): Cannot read point data array
> > > > > "U" from PointData in piece 0. The data array in the element may be 
> > > > > too
> > > > > short.
> > > > >
> > > > > I noticed that the reader crashes on some absurd numbers appearing in
> > > > > PointData "U" as
> > > > >
> > > > > [...]
> > > > > -3.949617e-105 9.769857e-104 0.000000e+00
> > > > > [...]
> > > > >
> > > > > If i set those to zero, there is no problem. but the appearance of 
> > > > > such
> > > > > 'almost-zero' numbers is clearly connected to using the parameter C.
> > > > >
> > > > > do you suggest me another way how to do this?
> > > > >
> > > > > or do you think that the VTK output writing should check for numbers <
> > > > > DOLFIN_EPS to avoid such corrupted output, and the approach is ok?
> > > > >
> > > > > best regards and thanks for the support,
> > > > > patrick
> > > > >
> > > > > > --
> > > > > > Anders
> > > > >
> > > > >
> > > > > _______________________________________________
> > > > > Mailing list: https://launchpad.net/~dolfin
> > > > > Post to : dolfin@xxxxxxxxxxxxxxxxxxx
> > > > > Unsubscribe : https://launchpad.net/~dolfin
> > > > > More help : https://help.launchpad.net/ListHelp
> > >
> > >
> > > _______________________________________________
> > > Mailing list: https://launchpad.net/~dolfin
> > > Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> > > Unsubscribe : https://launchpad.net/~dolfin
> > > More help   : https://help.launchpad.net/ListHelp
> >
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help   : https://help.launchpad.net/ListHelp