Commit 376a9274 authored by MARTIN Antoine's avatar MARTIN Antoine

add elliptical shape

parent c303da6d
......@@ -230,11 +230,14 @@ N0 = 128
d=2
N=N0*np.ones(d,dtype=np.int64)
radius = 8
shape = {2:"disc",3:"sphere"}
a=14
b=3
#shape = {2:"disc",3:"sphere"}
shape = {2:"ellipse",3:"ellipsoïd"}
V = pi*radius**2 if d==2 else 4*pi/3*radius**3
k1 = 1.
k2 = 0.001
for model in ["B"]:
for model in ["A"]:
if model == "A":
phaseProps = np.array([k1,k2],dtype=np.float64)
elif model == "B":
......@@ -244,11 +247,16 @@ for model in ["B"]:
kCA_list=[]
kSA_list=[]
kref_list=[]
for target_phi in [98,96,94,92,90,88,86,84,78,77,76,66,65]:
ncenter = int(-np.log(1-target_phi/100.)*N0**d/V)
for target_phi in range(2,99,5):
#ncenter = int(-np.log(1-target_phi/100.)*N0**d/V)
ncenter = int(-np.log(1-target_phi/100.)*N0**d/np.pi/a/b)
centers = np.random.rand(ncenter,d)*N0
radii = np.ones(ncenter)*radius
microstructure = booleanSpheres(centers,radii,N)
#radii = np.ones(ncenter)*radius
ai=np.ones(ncenter)*a
bi=np.ones(ncenter)*b
thetai = np.random.rand(ncenter,)*np.pi
#microstructure = booleanSpheres(centers,radii,N)
microstructure = booleanEllipses(centers,ai,bi,thetai,N)
#microstructure = booleanSpheres(np.array([[0,0],[N0/3,N0/3],[3*N0/4,N0/2]]),[N0/3,N0/3,N0/3],N)
mapStrides = np.array([N0**(d-i-1) for i in range(d)],dtype=np.int64)
phasemap = microstructure.phasemap.flatten() #TODO : flatten to allow scalar index (see numba issue with tuple index)
......@@ -333,24 +341,24 @@ plt.show()
plt.close()
x=np.linspace(0,1,101)
plt.plot(x,Kmoritanaka(x,k1,k2,d),label="Hashin-Shtrikman")
plt.plot(x,Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
plt.plot(x,[Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
# x=np.linspace(0,1,101)
# plt.plot(x,Kmoritanaka(x,k1,k2,d),label="Hashin-Shtrikman")
# plt.plot(x,Kmoritanaka(1-x,k2,k1,d),label="Hashin-Shtrikman")
# plt.plot(x,[Kselfconsistent(f2,k1,k2,d) for f2 in x],label="Self-Consistent")
resB = np.loadtxt("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resB[:,0],resB[:,1],'k.',label="hierachical num. self-consistent B")
resBSA = np.loadtxt("boolean_%s_contrast%dmodel%sSA.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resBSA[:,0],resBSA[:,1],'k.',color='red',label="hierachical num. self-consistent B SA")
resBCA = np.loadtxt("boolean_%s_contrast%dmodel%sCA.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resBCA[:,0],resBCA[:,1],'k.',color='green',label="hierachical num. self-consistent B CA")
resrefB = np.loadtxt("boolean_%s_contrast%dmodel%sref.txt" % (shape[d],int(k1/k2),"B"))
plt.plot(1-resrefB[:,0],resrefB[:,1],'k+',color='red',label="B ref")
# resB = np.loadtxt("boolean_%s_contrast%dmodel%s.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resB[:,0],resB[:,1],'k.',label="hierachical num. self-consistent B")
# resBSA = np.loadtxt("boolean_%s_contrast%dmodel%sSA.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resBSA[:,0],resBSA[:,1],'k.',color='red',label="hierachical num. self-consistent B SA")
# resBCA = np.loadtxt("boolean_%s_contrast%dmodel%sCA.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resBCA[:,0],resBCA[:,1],'k.',color='green',label="hierachical num. self-consistent B CA")
# resrefB = np.loadtxt("boolean_%s_contrast%dmodel%sref.txt" % (shape[d],int(k1/k2),"B"))
# plt.plot(1-resrefB[:,0],resrefB[:,1],'k+',color='red',label="B ref")
plt.legend(loc="upper right")
plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
plt.xlabel("volume fraction of insulating phase")
plt.ylabel(r"$K_{hom} / K_{conductive \, phase}$")
plt.savefig("boolean_%s_contrast%d_model%s.pdf" % (shape[d],int(k1/k2),"B"))
plt.show()
plt.close()
# plt.legend(loc="upper right")
# plt.title(r"boolean %s contrast = %d" % (shape[d],int(k1/k2)))
# plt.xlabel("volume fraction of insulating phase")
# plt.ylabel(r"$K_{hom} / K_{conductive \, phase}$")
# plt.savefig("boolean_%s_contrast%d_model%s.pdf" % (shape[d],int(k1/k2),"B"))
# plt.show()
# plt.close()
......@@ -43,3 +43,33 @@ class booleanSpheres(object):
return self.phasemap[n]
class booleanEllipses(object): #only in 2D
def __init__(self,centers,ai,bi,thetai,N):
self.centers=centers
self.ai=ai
self.bi=bi
self.thetai=thetai
self.N = N
self.phasemap = np.zeros(tuple(N),dtype=np.int8)
self.initializePhaseMap()
def initializePhaseMap(self):
for i,center in enumerate(self.centers):
a=self.ai[i]
b=self.bi[i]
m=max(a,b)
theta=self.thetai[i]
c = np.rint(center).astype(int)
for s in np.ndindex(tuple(2*[2*int(m)+1])):
ds = np.array(s)-int(m)
n= c-ds
x=n+0.5
A=np.array([[np.cos(theta),np.sin(theta)],[-np.sin(theta),np.cos(theta)]])
x=np.dot(A,x-center)
if (x[0])**2/a**2+(x[1])**2/b**2 <= 1:
self.phasemap[tuple(n % self.N)]=1
def get_phase(self,n):
return self.phasemap[n]
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment