From b3e8b85d25dd4d7676e6104ee8c8d80a69b47d5c Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 12 Nov 2020 18:34:29 -0500 Subject: [PATCH 1/2] normalizing coordinates _before_ rotation causes wrong primitive shape for non-cubic VEs --- python/damask/_geom.py | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index d6e0c1b9e..3d73861ad 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -564,14 +564,14 @@ class Geom: If given as floats, coordinates are addressed. center : int or float numpy.ndarray of shape (3) Center of the primitive. If given as integers, grid point - locations (cell centers) are addressed. - If given as floats, coordinates are addressed. + coordinates (cell centers) are addressed. + If given as floats, coordinates in space are addressed. exponent : numpy.ndarray of shape (3) or float Exponents for the three axes. 0 gives octahedron (ǀxǀ^(2^0) + ǀyǀ^(2^0) + ǀzǀ^(2^0) < 1) 1 gives sphere (ǀxǀ^(2^1) + ǀyǀ^(2^1) + ǀzǀ^(2^1) < 1) fill : int, optional - Fill value for primitive. Defaults to material.max() + 1. + Fill value for primitive. Defaults to material.max()+1. R : damask.Rotation, optional Rotation of primitive. Defaults to no rotation. inverse : Boolean, optional @@ -581,25 +581,27 @@ class Geom: Repeat primitive over boundaries. Defaults to True. """ - # normalized 'radius' and center - r = np.array(dimension)/self.grid/2.0 if np.array(dimension).dtype in np.sctypes['int'] else \ - np.array(dimension)/self.size/2.0 - c = (np.array(center) + .5)/self.grid if np.array(center).dtype in np.sctypes['int'] else \ - (np.array(center) - self.origin)/self.size + # radius and center + r = np.array(dimension)/2.0*self.size/self.grid if np.array(dimension).dtype in np.sctypes['int'] else \ + np.array(dimension)/2.0 + c = (np.array(center) + .5)*self.size/self.grid if np.array(center).dtype in np.sctypes['int'] else \ + (np.array(center) - self.origin) - coords = grid_filters.cell_coord0(self.grid,np.ones(3)) \ - - ((np.ones(3)-(1./self.grid if np.array(center).dtype in np.sctypes['int'] else 0))*0.5 if periodic else c) # periodic center is always at CoG + coords = grid_filters.cell_coord0(self.grid,self.size, + -(0.5*(self.size + (self.size/self.grid + if np.array(center).dtype in np.sctypes['int'] else + 0)) if periodic else c)) coords_rot = R.broadcast_to(tuple(self.grid))@coords with np.errstate(all='ignore'): mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0 if periodic: # translate back to center - mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2)) + mask = np.roll(mask,((c/self.size-0.5)*self.grid).round().astype(int),(0,1,2)) - fill_ = np.full_like(self.material,np.nanmax(self.material)+1 if fill is None else fill) - - return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_), + return Geom(material = np.where(np.logical_not(mask) if inverse else mask, + self.material, + np.nanmax(self.material)+1 if fill is None else fill), size = self.size, origin = self.origin, comments = self.comments+[util.execution_stamp('Geom','add_primitive')], From 7b7b449877c349a220ed91195f585bcefb2142d6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Nov 2020 06:50:57 +0100 Subject: [PATCH 2/2] testing fixed behavior --- python/tests/test_Geom.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 5c2777da7..a9b9c91b5 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -265,12 +265,12 @@ class TestGeom: @pytest.mark.parametrize('inverse',[True,False]) @pytest.mark.parametrize('periodic',[True,False]) def test_add_primitive_rotation(self,center,inverse,periodic): - """Rotation should not change result for sphere (except for discretization errors).""" - g = np.array([32,32,32]) + """Rotation should not change result for sphere.""" + g = np.random.randint(8,32,(3)) + s = np.random.random(3)+.5 fill = np.random.randint(10)+2 - eu=np.array([np.random.randint(4),np.random.randint(2),np.random.randint(4)])*.5*np.pi - G_1 = Geom(np.ones(g,'i'),[1.,1.,1.]).add_primitive(.3,center,1,fill,inverse=inverse,periodic=periodic) - G_2 = Geom(np.ones(g,'i'),[1.,1.,1.]).add_primitive(.3,center,1,fill,Rotation.from_Eulers(eu),inverse,periodic=periodic) + G_1 = Geom(np.ones(g,'i'),s).add_primitive(.3,center,1,fill,inverse=inverse,periodic=periodic) + G_2 = Geom(np.ones(g,'i'),s).add_primitive(.3,center,1,fill,Rotation.from_random(),inverse,periodic=periodic) assert geom_equal(G_1,G_2)