dolfin team mailing list archive

[Harish Narayanan <harish.mlists@xxxxxxxxx>]

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
}

Harish