I just tried your script with some data I have here. I've attached the images. I can't really tell if this is right or not. I think another way to test this would be to do a projection with a periodic sphere as the source object. If you throw a sphere that overlaps the domain boundary to a projection, it should be quite clear if it's doing what it's supposed to. Does that sound right?<br>
<br>Britton<br><br><div class="gmail_quote">On Fri, Apr 17, 2009 at 8:17 PM, Matthew Turk <span dir="ltr"><<a href="mailto:matthewturk@gmail.com">matthewturk@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Oops, just as a note, there's an error in the paste -- the second plot<br>
of "NewRadius" should be of "Radius".<br>
<br>
Thanks guys...<br>
<br>
-Matt<br>
<div><div></div><div class="h5"><br>
On Fri, Apr 17, 2009 at 7:16 PM, Matthew Turk <<a href="mailto:matthewturk@gmail.com">matthewturk@gmail.com</a>> wrote:<br>
> Hi guys,<br>
><br>
> Can somebody review this? I chatted with Jeff, and I believe it is<br>
> correct for periodic boxes, but I need somebody to look at it before I<br>
> commit it, as it is a rather crucial aspect of yt. This changes the<br>
> way radius is calculated and should now work with periodic boundary<br>
> conditions.<br>
><br>
> Here's the new field definition:<br>
><br>
> def _NewRadius(field, data):<br>
> center = data.get_field_parameter("center")<br>
> DW = <a href="http://data.pf" target="_blank">data.pf</a>["DomainRightEdge"] - <a href="http://data.pf" target="_blank">data.pf</a>["DomainLeftEdge"]<br>
> radius = na.zeros(data["x"].shape, dtype='float64')<br>
> for i, ax in enumerate('xyz'):<br>
> r = na.abs(data[ax] - center[i])<br>
> radius += na.minimum(r, na.abs(DW[i]-r))**2.0<br>
> na.sqrt(radius, radius)<br>
> return radius<br>
><br>
> the final argument in na.sqrt is the 'out' argument, to avoid temporary arrays.<br>
><br>
> I've committed a sample script here:<br>
><br>
> <a href="http://paste.enzotools.org/show/101/" target="_blank">http://paste.enzotools.org/show/101/</a><br>
><br>
> which loads up a data dump, calculates the radius from [0.1, 0.1, 0.1]<br>
> at all points in the root grid, then plots the output. It gives<br>
> correct min & max values for the radius.<br>
><br>
> If I can get somebody else to sign off on this, I will commit the<br>
> necessary changes to _Radius and _ParticleRadius, as well as grid &<br>
> point selection for spheres. We'll then have a fully-periodic<br>
> profiler, which will close #165.<br>
><br>
> -Matt<br>
><br>
> On Thu, Apr 16, 2009 at 11:39 AM, Matthew Turk <<a href="mailto:matthewturk@gmail.com">matthewturk@gmail.com</a>> wrote:<br>
>> Hi guys,<br>
>><br>
>> I've been exploring the idea of changing radius to be correct for<br>
>> periodic boxes. Right now it is incorrect; each component (x,y,z)<br>
>> should not have a distance greater than 0.5 * domain_size. The<br>
>> easiest way to do this would be:<br>
>><br>
>> rx = min( abs( x - center_x ), abs( x - center_x - domain_x) )<br>
>><br>
>> I think the best way to do this is to set up a NumPy ufunc<br>
>><br>
>> <a href="http://docs.scipy.org/doc/numpy/user/c-info.beyond-basics.html#creating-a-new-universal-function" target="_blank">http://docs.scipy.org/doc/numpy/user/c-info.beyond-basics.html#creating-a-new-universal-function</a><br>
>><br>
>> that accepts three arrays, along with the center and the domain width,<br>
>> and then returns the radius. What do you all think? Alternatively,<br>
>> I'm thinking maybe just a regular function that gets the arrays and<br>
>> returns one would be better; the ufunc machinery is a bit complicated<br>
>> and I might get confused. Once I come up with it, will somebody be<br>
>> able to look over my work?<br>
>><br>
>> -Matt<br>
>><br>
><br>
_______________________________________________<br>
Yt-dev mailing list<br>
<a href="mailto:Yt-dev@lists.spacepope.org">Yt-dev@lists.spacepope.org</a><br>
<a href="http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org</a><br>
</div></div></blockquote></div><br>