[yt-users] Computing the surface area of an iso-quantity surface

trobolo dinni trobolo.trobolo.dinni5 at gmail.com
Mon Feb 3 16:49:31 PST 2014


Hi Matt,

it is not a YT deficiency, it is me doing something wrong!

For testing purposes, I am using a spherically symmetric distribution
(centered in the box center) for which I know precisely the surface and the
radius at the density cutting point, so that I can compare the areas
obtained in different ways. I am trying to do everything in code units,
hence I think is not a units problem, but something wrong in my computation.
The values I get are the following:

1. by summing up the triangles areas: A1=0.0411385641439

2. by using the geometric formula for the sphere area: A2=0.130534407181

3. by using the geometric formula for the sphere area and as radius the
distance of one of the triangles vertexes from the center of the box:
A3=0.129874695188

So for sure yt is cutting correctly where I want, since the last two values
are similar. Then the error is in my procedure.
The difference in surface is not constant, I tried spheres with various
radii and A1 is always smaller than A2 and A3 of a variable factor.

Finally, I tried the flux with ones
(surface.calculate_flux("Ones","Ones","Ones") ) and the value I get
is: -7.95604693874e-05, according to the yt docs this should be in code
units as well.
I checked the implementation of the function march_cubes_grid_flux in
marching_cubes.pyx (attached below), and it does exactly what I am trying
to do, that is using Heron's formula to get the areas starting from the
coordinates of the vertices:

*                        # Now we need the area of the triangle.  This will
take*
*                        # a lot of time to calculate compared to the rest.*
*                        # We use Heron's formula.*
*                        for n in range(3):*
*                            fv[n] = 0.0*
*                        for n in range(3):*
*                            fv[0] += (current.p[0][n] -
current.p[2][n])**2.0*
*                            fv[1] += (current.p[1][n] -
current.p[0][n])**2.0*
*                            fv[2] += (current.p[2][n] -
current.p[1][n])**2.0*
*                        s = 0.0*
*                        for n in range(3):*
*                            fv[n] = fv[n]**0.5*
*                            s += 0.5 * fv[n]*
*                        area = (s*(s-fv[0])*(s-fv[1])*(s-fv[2]))*
*                        area = area**0.5*
*                        flux += temp*area*

so the result confuses me. I expected a result similar to A2 and A3.
Can I ask if you have any other suggestion about this?

Thanks for the help,
                                 Roberto










On 4 February 2014 02:19, Matthew Turk <matthewturk at gmail.com> wrote:

> Hi Roberto,
>
> On Sun, Feb 2, 2014 at 10:49 PM, trobolo dinni
> <trobolo.trobolo.dinni5 at gmail.com> wrote:
> > Dear yt users,
> >
> > I am trying to use the pf.h.surface routine, for example on density:
> >
> > surface=pf.h.surface(primary,"Density",external_density)
> >
> > and then using the sum of the triangles area (via surface.triangles) as
> an
> > estimate of the surface area. But this method is giving me a wrong value
> in
> > a test dataset, where I already know the right value for the surface
> area.
> >
> > I would like hence to ask if there is a better way to estimate the
> surface
> > area of an iso-quantity surface, or if the process I am following is
> wrong.
> > Below are the commands I am using to estimate the surface:
> >
> > #computes the surface of the equidensity surface at the desired density
> > #obtains the vertices of the triangles composing the surface tassellation
> > approximation
> > surface=pf.h.surface(primary,"Density",density)
> > surface_density = surface['Density']
> > surface_trg_vertices = surface.triangles
> > #rearranges the arrays with the vertex of the triangles (coordinates
> names:
> > a, b, c) to simplify the vectorial computation
> > a_x = surface_trg_vertices[:,0][:,0]
> > a_y = surface_trg_vertices[:,0][:,1]
> > a_z = surface_trg_vertices[:,0][:,2]
> > b_x = surface_trg_vertices[:,1][:,0]
> > b_y = surface_trg_vertices[:,1][:,1]
> > b_z = surface_trg_vertices[:,1][:,2]
> > c_x = surface_trg_vertices[:,2][:,0]
> > c_y = surface_trg_vertices[:,2][:,1]
> > c_z = surface_trg_vertices[:,2][:,2]
> > #finds the sides lengths of the triangles (an array for each set of side)
> > side_ab = ((a_x - b_x)**2.0 + (a_y - b_y)**2.0 + (a_z - b_z)**2.0)**0.5
> > side_bc = ((b_x - c_x)**2.0 + (b_y - c_y)**2.0 + (b_z - c_z)**2.0)**0.5
> > side_ca = ((c_x - a_x)**2.0 + (c_y - a_y)**2.0 + (c_z - a_z)**2.0)**0.5
> > #computes the altitude of the triangles wrt the ab side
> > s = 0.5 * (side_ab + side_bc + side_ca) #semiperimeter
> > altitude_side_ab = (2.0 * (s * (s - side_ab) * (s - side_bc) * (s -
> > side_ca))**0.5) / side_ab
> > #computes the area of the triangles using ab sides and altitudes
> > trg_surface_area = (side_ab * altitude_side_ab) / 2.0
> > #computes the total surface area
> > total_surface_area = np.sum(trg_surface_area)
>
> Hmm, I'm not sure.  Inside the flux computation we do the computation,
> and it's part of our testing suite, so I am not aware of any
> deficiencies on the yt side; I am not sure I am following your code
> completely to get the answer.  But one immediate thing you might want
> to be careful of is how the units are converted, as the units you get
> back will be in code_length**2.  What is it off by?  A constant, or a
> factor?
>
> You can also try taking the flux of the fields "ones" over your
> surface to see what yt's internal area computation will return.
>
> -Matt
>
> >
> > Thanks,
> >                   Roberto
> >
> >
> > _______________________________________________
> > yt-users mailing list
> > yt-users at lists.spacepope.org
> > http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
> >
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20140204/0471745b/attachment.htm>


More information about the yt-users mailing list