Computing the surface area of an iso-quantity surface
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)* Thanks, Roberto
Hi Roberto, On Sun, Feb 2, 2014 at 10:49 PM, trobolo dinni <trobolo.trobolo.dinni5@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
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@gmail.com> wrote:
Hi Roberto,
On Sun, Feb 2, 2014 at 10:49 PM, trobolo dinni <trobolo.trobolo.dinni5@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
_______________________________________________ yt-users mailing list yt-users@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
participants (2)
-
Matthew Turk
-
trobolo dinni