223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'fix-add_primitive-rotation' into 'development' 

Avoids coordinate normalization _before_ rotation

See merge request damask/DAMASK!276
This commit is contained in:
Martin Diehl 2020-11-13 14:11:54 +01:00
parent 1f59f63011 7b7b449877
commit bacfe9f462
2 changed files with 21 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
python/damask/_geom.py
View File

@ -564,14 +564,14 @@ class Geom:
If given as floats, coordinates are addressed. If given as floats, coordinates are addressed.
center : int or float numpy.ndarray of shape (3) center : int or float numpy.ndarray of shape (3)
Center of the primitive. If given as integers, grid point Center of the primitive. If given as integers, grid point
locations (cell centers) are addressed. coordinates (cell centers) are addressed.
If given as floats, coordinates are addressed. If given as floats, coordinates in space are addressed.
exponent : numpy.ndarray of shape (3) or float exponent : numpy.ndarray of shape (3) or float
Exponents for the three axes. Exponents for the three axes.
0 gives octahedron (ǀxǀ^(2^0) + ǀyǀ^(2^0) + ǀzǀ^(2^0) < 1) 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) 1 gives sphere (ǀxǀ^(2^1) + ǀyǀ^(2^1) + ǀzǀ^(2^1) < 1)
fill : int, optional 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 R : damask.Rotation, optional
Rotation of primitive. Defaults to no rotation. Rotation of primitive. Defaults to no rotation.
inverse : Boolean, optional inverse : Boolean, optional
@ -581,25 +581,27 @@ class Geom:
Repeat primitive over boundaries. Defaults to True. Repeat primitive over boundaries. Defaults to True.
""" """
# normalized 'radius' and center # radius and center
r = np.array(dimension)/self.grid/2.0 if np.array(dimension).dtype in np.sctypes['int'] else \ r = np.array(dimension)/2.0*self.size/self.grid if np.array(dimension).dtype in np.sctypes['int'] else \
np.array(dimension)/self.size/2.0 np.array(dimension)/2.0
c = (np.array(center) + .5)/self.grid if np.array(center).dtype in np.sctypes['int'] else \ c = (np.array(center) + .5)*self.size/self.grid if np.array(center).dtype in np.sctypes['int'] else \
(np.array(center) - self.origin)/self.size (np.array(center) - self.origin)
coords = grid_filters.cell_coord0(self.grid,np.ones(3)) \ coords = grid_filters.cell_coord0(self.grid,self.size,
- ((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 -(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 coords_rot = R.broadcast_to(tuple(self.grid))@coords
with np.errstate(all='ignore'): with np.errstate(all='ignore'):
mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0 mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0
if periodic: # translate back to center 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,
return Geom(material = np.where(np.logical_not(mask) if inverse else mask, self.material,fill_), np.nanmax(self.material)+1 if fill is None else fill),
size = self.size, size = self.size,
origin = self.origin, origin = self.origin,
comments = self.comments+[util.execution_stamp('Geom','add_primitive')], comments = self.comments+[util.execution_stamp('Geom','add_primitive')],

10
python/tests/test_Geom.py
View File

@ -265,12 +265,12 @@ class TestGeom:
@pytest.mark.parametrize('inverse',[True,False]) @pytest.mark.parametrize('inverse',[True,False])
@pytest.mark.parametrize('periodic',[True,False]) @pytest.mark.parametrize('periodic',[True,False])
def test_add_primitive_rotation(self,center,inverse,periodic): def test_add_primitive_rotation(self,center,inverse,periodic):
"""Rotation should not change result for sphere (except for discretization errors).""" """Rotation should not change result for sphere."""
g = np.array([32,32,32]) g = np.random.randint(8,32,(3))
s = np.random.random(3)+.5
fill = np.random.randint(10)+2 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'),s).add_primitive(.3,center,1,fill,inverse=inverse,periodic=periodic)
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'),s).add_primitive(.3,center,1,fill,Rotation.from_random(),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)
assert geom_equal(G_1,G_2) assert geom_equal(G_1,G_2)