dolfin team mailing list archive

[Patrick Riesen <priesen@xxxxxxxxxxxxxxx>]

Harish Narayanan wrote:
> On 5/17/10 9:50 AM, 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.
> >
> > 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?
>
> I have seen this for a while and use the following bash script to get
> around this error:
>
> 1. cleanup.sh:
>
> #!/usr/bin/env bash
>
> for file in *.vtu
> do
>     gawk -f truncate.gawk $file  > tmp_$file
>     mv tmp_$file $file
>
>
> 2. truncate.gawk:
>
> function truncate(x)
> {
>     EPS=1.e-16
>     absx = x >= 0 ? x : -x
>     if (absx < EPS) {
> 	return "0.0"
>     }
>     else {
> 	return x
>     }
> }
>
> {
>     if (NF == 1)
> 	$i = truncate($i);
>     print
> }

cool! thank you, harish!
regards,
patrick

> Harish