Johan Hoffman wrote:
Patrick Riesen wrote:
Johan Hoffman wrote:
Patrick Riesen wrote:
Patrick Riesen wrote:
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.
Ok. You can check in the Unicorn solver, this is what we use:
http://www.fenics.org/wiki/Unicorn
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?
do you think i can keep the weak statement and just replace the
facet normal 'n' (the facet normal from ufl) by a linear normal
function of the vertex normals?
Yes, I would expect this to work.
ok, i think it did not work (or not much better).
I might try the strong way.
But how to specify strong Dirichlet conditions for a variable 'v_m' as
a function of an other variable 'v' in the mixed problem (v_m, v) in
fenics?
We apply this type of boundary conditions to the discrete matrix
equation after having assembled the system, see the implementation in
Unicorn. This is not very neat, and I think you instead probably would
like to apply strong bc through the definition of the function space.
so far i came up with something similar. i use the dirichletBC class
together with a special expression which sets the vertex values on the
free surface to v_m = v in the space of v_m, i.e the boundary is fully
lagrangian. This is more restrictive than specifying the condition on
the normal components, but was easier to implement for initial testing.
i here plotted the element-wise area differences across timestepping:
http://people.ee.ethz.ch/~priesen/public/tmsoim_ale.gif
i still get a total cumulative area loss of about 10e-2 from the O(1)
square. I'm not sure if this is acceptable. although, the deformation is
rather significant. should i be happy or not? ;-)
an other approach on the formulation level which might be possible is to
use the Lifting functionality in ufl for the dirichlet moving boundary
part. i will look into that.
thanks for the help so far,
regards,
patrick
Best,
Johan
regards,
patrick
Best,
Johan
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@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
