Hi David,<br><br>Hmm, this smells similar to this bug:<div><br></div><div><a href="https://bitbucket.org/yt_analysis/yt/issue/854/different-dimensions-for-coordinates-and">https://bitbucket.org/yt_analysis/yt/issue/854/different-dimensions-for-coordinates-and</a><br></div><div><br></div><div>I think that there is a rounding error in the logic for constructing a "smoothed" covering grid, but I have so far been unable to find it. What the smoothing does is attempt to construct level-by-level cascading interpolation in the ghost regions. To do that, it tries to expand by a dimension each level, to make sure there's enough room for a trilinear (cell-centered) interpolation. I believe the issue is likely that at some point, it's rounding incorrectly. Would you mind bumping that bug report, and then when I have a bit of time (unfortunately the bus factor on the smoothed covering grid is rather low...) I will take a look, hopefully today?</div><div><br></div><div>-Matt</div><br><div class="gmail_quote">On Tue Nov 11 2014 at 7:38:01 AM David Ketcheson <<a href="mailto:dketch@gmail.com">dketch@gmail.com</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">I'm trying to plot isosurfaces from AMRClaw data. I've been able to successfully plot slices and projections. Here is the script I'm using for isosurfaces:<div><br></div><div><a href="https://gist.github.com/ketch/30d8bdf41422dc20bb95" target="_blank">https://gist.github.com/ketch/30d8bdf41422dc20bb95</a><br></div><div><br></div><div>(based on <a href="http://yt-project.org/doc/cookbook/complex_plots.html#plotting-isocontours" target="_blank">http://yt-project.org/doc/cookbook/complex_plots.html#plotting-isocontours</a>)</div><div><br></div><div>I get the following error:</div><div><br></div><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">---------------------------------------------------------------------------</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">ValueError Traceback (most recent call last)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">/Users/ketch/Downloads/plot_yt_amrclaw.py in <module>()</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 35 ad = ds.all_data()</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 36 surface = ds.h.surface(ad,"q",0.2)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">---> 37 p3dc = Poly3DCollection(surface.triangles, linewidth=0)#, linewidth=0.1,alpha=0.9)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 38 colors = yt.apply_colormap((surface["q"]), cmap_name="hot")</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 39 p3dc.set_facecolors(colors[0,:,:]/255.)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><br></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">/Users/ketch/anaconda/lib/python2.7/site-packages/yt/data_objects/construction_data_containers.pyc in triangles(self)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 1021 def triangles(self):</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 1022 if self.vertices is None:</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">-> 1023 self.get_data()</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 1024 vv = np.empty((self.vertices.shape[1]/3, 3, 3), dtype="float64")</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 1025 for i in range(3):</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><br></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">/Users/ketch/anaconda/lib/python2.7/site-packages/yt/data_objects/construction_data_containers.pyc in get_data(self, fields, sample_type)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 913 my_verts = self._extract_isocontours_from_grid(</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 914 block, self.surface_field, self.field_value,</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">--> 915 mask, fields, sample_type)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 916 if fields is not None:</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 917 my_verts, svals = my_verts</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><br></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">/Users/ketch/anaconda/lib/python2.7/site-packages/yt/data_objects/construction_data_containers.pyc in _extract_isocontours_from_grid(self, grid, field, value, mask, sample_values, sample_type)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 933 mask, sample_values = None,</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 934 sample_type = "face"):</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">--> 935 vals = grid.get_vertex_centered_data(field, no_ghost = False)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 936 if sample_values is not None:</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 937 svals = grid.get_vertex_centered_data(sample_values)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><br></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">/Users/ketch/anaconda/lib/python2.7/site-packages/yt/data_objects/grid_patch.py in get_vertex_centered_data(self, field, smoothed, no_ghost)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 289 else:</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 290 cg = self.retrieve_ghost_zones(1, field, smoothed=smoothed)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">--> 291 np.add(new_field, cg[field][1: ,1: ,1: ], new_field)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 292 np.add(new_field, cg[field][:-1,1: ,1: ], new_field)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"> 293 np.add(new_field, cg[field][1: ,:-1,1: ], new_field)</blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px"><br></blockquote><blockquote style="margin:0px 0px 0px 40px;border:none;padding:0px">ValueError: operands could not be broadcast together with shapes (47,47,47) (48,48,48) (47,47,47)</blockquote></blockquote><div><div><br></div><div>I get the same error if I just ask for surface.triangles.</div><div><br></div><div>After a bit of hunting, I found that the whole thing works fine if I set smoothed=False in line 290. Any ideas as to why the array size comes out wrong?</div><div><br></div></div></div>
______________________________<u></u>_________________<br>
yt-users mailing list<br>
<a href="mailto:yt-users@lists.spacepope.org" target="_blank">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/<u></u>listinfo.cgi/yt-users-<u></u>spacepope.org</a><br>
</blockquote></div>