Note
Click here to download the full example code or to run this example in your browser via Binder
Four phase material composite#
Calculating the effective conductivity/resistivity.
import logging
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize_scalar
from scipy.special import ellipk, ellipkm1, lpmv
import gyptis as gy
def _ellipk(m):
return ellipkm1(m) if np.abs(m - 1) < 1e-6 else ellipk(m)
logging.basicConfig(format="%(levelname)s:%(message)s", level=logging.WARNING)
np.set_printoptions(precision=4, suppress=True)
Results are compared with [Craster2001].
class FourPhaseComposite:
def __init__(self, l, h, rho):
self.h = h
self.l = l
self.rho = np.array(rho)
self.rho_dict = dict(
W1=rho[0],
W2=rho[1],
W3=rho[2],
W4=rho[3],
)
sigma1 = np.sum(rho)
sigma2 = rho[0] * rho[2] - rho[1] * rho[3]
sigma3 = (
rho[0] * rho[1] * rho[2]
+ rho[0] * rho[1] * rho[3]
+ rho[0] * rho[2] * rho[3]
+ rho[1] * rho[2] * rho[3]
)
self.sigma = np.array([sigma1, sigma2, sigma3])
self.Delta2 = sigma2**2 / (sigma1 * sigma3 + sigma2**2)
self.lambda_ = np.arccos(1 - 2 * self.Delta2) / np.pi
def build_geometry(self, lmin=0.05):
logging.info("Building geometry and mesh")
ax, ay = l, h # lattice constant
hx, hy = 0.5 * ax, 0.5 * ay
lattice = gy.Lattice(dim=2, vectors=((ax, 0), (0, ay)))
cell = lattice.cell
sq1 = lattice.add_rectangle(hx, hy, 0, ax - hx, ay - hy)
sq1, cell = lattice.fragment(sq1, cell)
sq2 = lattice.add_rectangle(hx, 0, 0, ax - hx, hy)
sq2, cell = lattice.fragment(sq2, cell)
sq3 = lattice.add_rectangle(0, 0, 0, hx, hy)
sq3, cell = lattice.fragment(sq3, cell)
sq4 = lattice.add_rectangle(0, hy, 0, hx, ay - hy)
sq4 = lattice.fragment(sq4, cell)
lattice.remove_all_duplicates()
lattice.add_physical(sq1, "W1")
lattice.add_physical(sq2, "W2")
lattice.add_physical(sq3, "W3")
lattice.add_physical(sq4, "W4")
lattice.set_size("W1", lmin)
lattice.set_size("W2", lmin)
lattice.set_size("W3", lmin)
lattice.set_size("W4", lmin)
lattice.build()
self.lattice = lattice
def _homogenize_numerical(self):
logging.info("Computing homogenization problem with FEM")
hom = gy.Homogenization2D(
self.lattice,
self.rho_dict,
degree=2,
)
eps_eff = hom.get_effective_permittivity()[:2, :2].real
rho_y = eps_eff[0, 0]
rho_x = eps_eff[1, 1]
self.model = hom
return rho_x, rho_y
def _homogenize_analytical(self):
logging.info("Computing homogenization problem analytically")
m = self.compute_m()
self.m = m
ratio_sigma = (
self.l
/ self.h
* lpmv(0, 0.5 * (self.lambda_ - 1), 2 * m - 1)
/ lpmv(0, 0.5 * (self.lambda_ - 1), 1 - 2 * m)
)
Q = (
((rho[1] + rho[2]) * (rho[3] + rho[0]))
/ ((rho[0] + rho[1]) * (rho[2] + rho[3]))
) ** 0.5
p = (self.sigma[2] / self.sigma[0]) ** 0.5
rho_x_ana = 1 / ratio_sigma * Q * p
rho_y_ana = ratio_sigma / Q * p
return rho_x_ana, rho_y_ana
def check_product(self, rho, rtol=1e-3):
rho_x, rho_y = rho
assert np.allclose(rho_x * rho_y, self.sigma[2] / self.sigma[0], rtol=rtol)
def compute_m(self):
logging.info("Computing parameter m")
def objective(m):
# m = m[0]
obj = np.abs(_ellipk(m) / _ellipk(1 - m) - self.l / self.h)
logging.info(f"m = {m}, objective = {obj}")
return obj
logging.info(" Root finding for m")
opt = minimize_scalar(objective, bounds=(0, 1), method="bounded")
logging.info(opt)
m = opt.x
return m
def homogenize(self, method):
if method in ["numerical", "analytical"]:
return (
self._homogenize_numerical()
if method == "numerical"
else self._homogenize_analytical()
)
else:
raise ValueError("Wrong method: choose between analytical or numerical.")
Loop over lattice size along x:
L = np.linspace(0.3, 1, 21)
h = 1
rho = [1, 4, 5, 10]
num = []
ana = []
for l in L:
prob = FourPhaseComposite(l, h, rho)
prob.build_geometry(lmin=0.05)
rho_x, rho_y = prob.homogenize("numerical")
rho_x_ana, rho_y_ana = prob.homogenize("analytical")
prob.check_product([rho_x, rho_y], rtol=1e-2)
print(rho_x, rho_y)
print(rho_x_ana, rho_y_ana)
assert np.allclose(rho_x, rho_x_ana, rtol=1e-2)
assert np.allclose(rho_y, rho_y_ana, rtol=1e-2)
num.append([rho_x, rho_y])
ana.append([rho_x_ana, rho_y_ana])
num = np.array(num).T
ana = np.array(ana).T
fig, ax = plt.subplots(1, 2)
ax[0].plot(L, ana[0], label="analytical")
ax[0].plot(L, num[0], "o", label="numerical")
ax[1].plot(L, ana[1])
ax[1].plot(L, num[1], "o")
ax[0].set_xlabel(r"$l/h$")
ax[1].set_xlabel(r"$l/h$")
ax[0].set_ylabel(r"$\rho_x$")
ax[1].set_ylabel(r"$\rho_y$")
ax[0].legend()
plt.tight_layout()
Out:
Level 25:Calling FFC just-in-time (JIT) compiler, this may take some time.
INFO:Compiling form ffc_form_3e233341a92a70dc476c1f7d4de9ac1497018eb1
INFO:Compiler stage 1: Analyzing form(s)
INFO:-----------------------------------
DEBUG: Preprocessing form using 'uflacs' representation family.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:
INFO: Geometric dimension: 2
Number of cell subdomains: 5
Rank: 2
Arguments: '(v_0, v_1)'
Number of coefficients: 2
Coefficients: '[f_9081, f_9082]'
Unique elements: 'Mixed<CG2(?,?), CG2(?,?)>, CG2(?,?), Vector<2 x CG1(?,
?)>'
Unique sub elements: 'Mixed<CG2(?,?), CG2(?,?)>, CG2(?,?), Vector<2 x CG1(?,
?)>, CG1(?,?)'
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO:
INFO:Compiler stage 1 finished in 0.0408006 seconds.
INFO:Compiler stage 2: Computing intermediate representation
INFO:-------------------------------------------------------
INFO: Computing representation of 0 elements
INFO: Computing representation of 0 dofmaps
INFO: Computing representation of 0 coordinate mappings
INFO: Computing representation of integrals
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
16 full
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
16 full
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
16 full
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
16 full
INFO: Computing representation of forms
INFO:
INFO:Compiler stage 2 finished in 0.0891988 seconds.
INFO:Compiler stage 3: Optimizing intermediate representation
INFO:--------------------------------------------------------
INFO: Optimizing uflacs representation
INFO: Optimizing uflacs representation
INFO: Optimizing uflacs representation
INFO: Optimizing uflacs representation
INFO:
INFO:Compiler stage 3 finished in 0.000209808 seconds.
INFO:Compiler stage 4: Generating code
INFO:---------------------------------
INFO: Generating code for 0 finite_element(s)
INFO: Generating code for 0 dofmap(s)
INFO: Generating code for 0 coordinate_mapping(s)
INFO: Generating code for integrals
INFO: Generating code from ffc.uflacs representation
INFO: Generating code from ffc.uflacs representation
INFO: Generating code from ffc.uflacs representation
INFO: Generating code from ffc.uflacs representation
INFO: Generating code for forms
INFO:
INFO:Compiler stage 4 finished in 0.0204771 seconds.
INFO:Compiler stage 4.1 finished in 2.6226e-06 seconds.
INFO:Compiler stage 5: Formatting code
INFO:---------------------------------
INFO:
INFO:Compiler stage 5 finished in 0.000459433 seconds.
INFO:FFC finished in 0.151397 seconds.
Level 25:Calling FFC just-in-time (JIT) compiler, this may take some time.
INFO:Compiling form ffc_form_bf1c1f953a7c1643fd6b2bc2c7ddb4572bdc535f
INFO:Compiler stage 1: Analyzing form(s)
INFO:-----------------------------------
DEBUG: Preprocessing form using 'uflacs' representation family.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:Adjusting missing element cell to triangle.
INFO:
INFO: Geometric dimension: 2
Number of cell subdomains: 5
Rank: 1
Arguments: '(v_0)'
Number of coefficients: 4
Coefficients: '[f_9096, f_9097, f_9106, f_9107]'
Unique elements: 'Mixed<CG2(?,?), CG2(?,?)>, CG2(?,?), Vector<2 x R0(?,?
)>, Vector<2 x CG1(?,?)>'
Unique sub elements: 'Mixed<CG2(?,?), CG2(?,?)>, CG2(?,?), Vector<2 x R0(?,?
)>, Vector<2 x CG1(?,?)>, R0(?,?), CG1(?,?)'
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO: representation: auto --> uflacs
INFO: quadrature_rule: auto --> default
INFO: quadrature_degree: 4
INFO:
INFO:Compiler stage 1 finished in 0.0294549 seconds.
INFO:Compiler stage 2: Computing intermediate representation
INFO:-------------------------------------------------------
INFO: Computing representation of 0 elements
INFO: Computing representation of 0 dofmaps
INFO: Computing representation of 0 coordinate mappings
INFO: Computing representation of integrals
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
4 full
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
4 full
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
4 full
INFO: Computing uflacs representation
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG: Reusing element from cache
DEBUG:Blocks of each mode:
4 full
INFO: Computing representation of forms
INFO:
INFO:Compiler stage 2 finished in 0.0835199 seconds.
INFO:Compiler stage 3: Optimizing intermediate representation
INFO:--------------------------------------------------------
INFO: Optimizing uflacs representation
INFO: Optimizing uflacs representation
INFO: Optimizing uflacs representation
INFO: Optimizing uflacs representation
INFO:
INFO:Compiler stage 3 finished in 0.00018239 seconds.
INFO:Compiler stage 4: Generating code
INFO:---------------------------------
INFO: Generating code for 0 finite_element(s)
INFO: Generating code for 0 dofmap(s)
INFO: Generating code for 0 coordinate_mapping(s)
INFO: Generating code for integrals
INFO: Generating code from ffc.uflacs representation
INFO: Generating code from ffc.uflacs representation
INFO: Generating code from ffc.uflacs representation
INFO: Generating code from ffc.uflacs representation
INFO: Generating code for forms
INFO:
INFO:Compiler stage 4 finished in 0.0124834 seconds.
INFO:Compiler stage 4.1 finished in 1.19209e-06 seconds.
INFO:Compiler stage 5: Formatting code
INFO:---------------------------------
INFO:
INFO:Compiler stage 5 finished in 0.000355959 seconds.
INFO:FFC finished in 0.126241 seconds.
4.797566778490196 3.2232183969024915
4.803310207742164 3.2269412820801144
4.782168890257173 3.235248574757013
4.787025356627207 3.2379189256939367
4.765704205116796 3.2462646926231615
4.770470263456939 3.2491555641241683
4.750051011282792 3.257813901516091
4.753905875500621 3.2604768386096277
4.73391005490445 3.268769138869265
4.737986112386069 3.2714321300942224
4.718005164758221 3.279703656152755
4.722016661971629 3.282495829552651
4.702852773156859 3.2909985004391396
4.706383806145529 3.293399059328803
4.686858147462029 3.3013161192163474
4.6910166490971195 3.304187803934418
4.672184549055963 3.312279507321037
4.675884045273595 3.314881175393447
4.657785070813897 3.323109733755722
4.661070752416577 3.325416159358636
4.643286161415294 3.3333923835901142
4.6466126541198385 3.335763308408587
4.629140264343214 3.3434857270745586
4.632523865565228 3.3459082888305423
4.615586265594989 3.3535818802208
4.618787663178187 3.3558589678345276
4.602410392574912 3.363397632618094
4.605441179984758 3.3655841849338968
4.589465569784161 3.3728608642005975
4.592482654531937 3.3750807931096594
4.577224011249757 3.3823439986048482
4.579941735525802 3.3843225296446957
4.565022081735377 3.391268475843669
4.567804845926401 3.393314846588292
4.553249334266513 3.3999516325676935
4.556067497891292 3.402056709470161
4.542040849649582 3.408539887329329
4.544741812492325 3.410534776121823
4.53124532089714 3.4168370467971907
4.533808640740242 3.4187591996536697
4.520671194565713 3.424750527896771
4.523274940517169 3.4267207286382213
Total running time of the script: ( 0 minutes 7.863 seconds)
Estimated memory usage: 14 MB