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

Matthew Turk matthewturk at gmail.com
Mon Feb 3 07:19:40 PST 2014


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
>



More information about the yt-users mailing list