<div dir="ltr">Hi Matt,<div><br></div><div>it is not a YT deficiency, it is me doing something wrong! </div><div><br></div><div>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.</div>
<div>The values I get are the following:</div><div><br></div><div>1. by summing up the triangles areas: A1=0.0411385641439</div><div><br></div><div>2. by using the geometric formula for the sphere area: A2=0.130534407181</div>
<div><br></div><div>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</div><div><br></div>
<div>So for sure yt is cutting correctly where I want, since the last two values are similar. Then the error is in my procedure. </div><div>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.</div>
<div><br></div><div>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.</div>
<div><span style="text-align:justify">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:</span></div>
<div style="text-align:justify"><br></div><div style="text-align:justify"><div><i> # Now we need the area of the triangle. This will take</i></div><div><i> # a lot of time to calculate compared to the rest.</i></div>
<div><i> # We use Heron's formula.</i></div><div><i> for n in range(3):</i></div><div><i> fv[n] = 0.0</i></div><div><i> for n in range(3):</i></div>
<div><i> fv[0] += (current.p[0][n] - current.p[2][n])**2.0</i></div><div><i> fv[1] += (current.p[1][n] - current.p[0][n])**2.0</i></div><div><i> fv[2] += (current.p[2][n] - current.p[1][n])**2.0</i></div>
<div><i> s = 0.0</i></div><div><i> for n in range(3):</i></div><div><i> fv[n] = fv[n]**0.5</i></div><div><i> s += 0.5 * fv[n]</i></div>
<div><i> area = (s*(s-fv[0])*(s-fv[1])*(s-fv[2]))</i></div><div><i> area = area**0.5</i></div><div><i> flux += temp*area</i></div></div>
<div><br></div><div>so the result confuses me. I expected a result similar to A2 and A3.</div><div>Can I ask if you have any other suggestion about this?</div><div><br></div><div>Thanks for the help,</div><div> Roberto</div>
<div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On 4 February 2014 02:19, Matthew Turk <span dir="ltr"><<a href="mailto:matthewturk@gmail.com" target="_blank">matthewturk@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Roberto,<br>
<div><div class="h5"><br>
On Sun, Feb 2, 2014 at 10:49 PM, trobolo dinni<br>
<<a href="mailto:trobolo.trobolo.dinni5@gmail.com">trobolo.trobolo.dinni5@gmail.com</a>> wrote:<br>
> Dear yt users,<br>
><br>
> I am trying to use the pf.h.surface routine, for example on density:<br>
><br>
> surface=pf.h.surface(primary,"Density",external_density)<br>
><br>
> and then using the sum of the triangles area (via surface.triangles) as an<br>
> estimate of the surface area. But this method is giving me a wrong value in<br>
> a test dataset, where I already know the right value for the surface area.<br>
><br>
> I would like hence to ask if there is a better way to estimate the surface<br>
> area of an iso-quantity surface, or if the process I am following is wrong.<br>
> Below are the commands I am using to estimate the surface:<br>
><br>
> #computes the surface of the equidensity surface at the desired density<br>
> #obtains the vertices of the triangles composing the surface tassellation<br>
> approximation<br>
> surface=pf.h.surface(primary,"Density",density)<br>
> surface_density = surface['Density']<br>
> surface_trg_vertices = surface.triangles<br>
> #rearranges the arrays with the vertex of the triangles (coordinates names:<br>
> a, b, c) to simplify the vectorial computation<br>
> a_x = surface_trg_vertices[:,0][:,0]<br>
> a_y = surface_trg_vertices[:,0][:,1]<br>
> a_z = surface_trg_vertices[:,0][:,2]<br>
> b_x = surface_trg_vertices[:,1][:,0]<br>
> b_y = surface_trg_vertices[:,1][:,1]<br>
> b_z = surface_trg_vertices[:,1][:,2]<br>
> c_x = surface_trg_vertices[:,2][:,0]<br>
> c_y = surface_trg_vertices[:,2][:,1]<br>
> c_z = surface_trg_vertices[:,2][:,2]<br>
> #finds the sides lengths of the triangles (an array for each set of side)<br>
> side_ab = ((a_x - b_x)**2.0 + (a_y - b_y)**2.0 + (a_z - b_z)**2.0)**0.5<br>
> side_bc = ((b_x - c_x)**2.0 + (b_y - c_y)**2.0 + (b_z - c_z)**2.0)**0.5<br>
> side_ca = ((c_x - a_x)**2.0 + (c_y - a_y)**2.0 + (c_z - a_z)**2.0)**0.5<br>
> #computes the altitude of the triangles wrt the ab side<br>
> s = 0.5 * (side_ab + side_bc + side_ca) #semiperimeter<br>
> altitude_side_ab = (2.0 * (s * (s - side_ab) * (s - side_bc) * (s -<br>
> side_ca))**0.5) / side_ab<br>
> #computes the area of the triangles using ab sides and altitudes<br>
> trg_surface_area = (side_ab * altitude_side_ab) / 2.0<br>
> #computes the total surface area<br>
> total_surface_area = np.sum(trg_surface_area)<br>
<br>
</div></div>Hmm, I'm not sure. Inside the flux computation we do the computation,<br>
and it's part of our testing suite, so I am not aware of any<br>
deficiencies on the yt side; I am not sure I am following your code<br>
completely to get the answer. But one immediate thing you might want<br>
to be careful of is how the units are converted, as the units you get<br>
back will be in code_length**2. What is it off by? A constant, or a<br>
factor?<br>
<br>
You can also try taking the flux of the fields "ones" over your<br>
surface to see what yt's internal area computation will return.<br>
<br>
-Matt<br>
<br>
><br>
> Thanks,<br>
> Roberto<br>
><br>
><br>
> _______________________________________________<br>
> yt-users mailing list<br>
> <a href="mailto:yt-users@lists.spacepope.org">yt-users@lists.spacepope.org</a><br>
> <a href="http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org</a><br>
><br>
_______________________________________________<br>
yt-users mailing list<br>
<a href="mailto:yt-users@lists.spacepope.org">yt-users@lists.spacepope.org</a><br>
<a href="http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org</a><br>
</blockquote></div><br></div>