Dear all,
I’m practicing fenics base on the problem “Periodic homogenization of linear elastic materials”. In this example, the considered 2D plane strain problem deals with a skewed unit cell of dimensions 1×3–√/21×3/2 consisting of circular inclusions (numbered 11) of radius RR with elastic properties (Er,νr)(Er,νr) and embedded in a matrix material (numbered 00) of properties (Em,νm)(Em,νm) following an hexagonal pattern.
I’m trying to redo this example but consider a rectangle unit cell of dimension of 1xsqrt(3) as following:
The mesh of this problem is created using Gmsh with the geometry as following:
// Gmsh project created on Mon Sep 14 09:33:57 2015
a = 1.;
R = 0.2;
t = Pi/3;
b = 2*Sin(t)*a;
d = 1.;
nx = 80;
nR = R*nx;
ns = nx-2*nR;
Point(1) = {0, 0, 0, d};
Point(2) = {a, 0, 0, d};
Point(3) = {a, b, 0, d};
Point(4) = {c, b, 0, d};
Point(5) = {R,0,0,d};
Point(6) = {a-R,0,0,d};
Point(7) = {a,R,0,d};
Point(8) = {a,b-R,0,d};
Point(9) = {a-R,b,0,d};
Point(10) = {R,b,0,d};
Point(11) = {0,b-R,0,d};
Point(12) = {0,R,0,d};
Point(13) = {a/2,b/2,0,d};
Line(1) = {1, 5};
Line(2) = {5, 6};
Line(3) = {6, 2};
Line(4) = {2, 7};
Line(5) = {7, 8};
Line(6) = {8, 3};
Line(7) = {3, 9};
Line(8) = {9, 10};
Line(9) = {10, 4};
Line(10) = {4, 11};
Line(11) = {11, 12};
Line(12) = {12, 1};
Circle(13) = {5, 1, 12};
Circle(14) = {7, 2, 6};
Circle(15) = {9, 3, 8};
Circle(16) = {11, 4, 10};
Line Loop(19) = {1, 13, 12};
Plane Surface(20) = {19};
Line Loop(21) = {3, 4, 14};
Plane Surface(22) = {21};
Line Loop(23) = {7, 15, 6};
Plane Surface(24) = {23};
Line Loop(25) = {10, 16, 9};
Plane Surface(26) = {25};
Transfinite Line{1,3,4,6,7,9,10,12} = nR;
Transfinite Line{2,5,8,11} = ns;
Transfinite Line{14,16} = nR;
Transfinite Line{13,15} = nR;
Transfinite Line{18,19,20,17} = nR;
Periodic Line{1} = {9};
Periodic Line{2} = {8};
Periodic Line{3} = {7};
Periodic Line{4} = {12};
Periodic Line{5} = {11};
Periodic Line{6} = {10};
Physical Line(1) = {1,2,3};
Physical Line(2) = {4,5,6};
Physical Line(3) = {7,8,9};
Physical Line(4) = {10,11,12};
Physical Surface(0) = {18};
Physical Surface(1) = {20,22,24,26,27};
//+
Point(14) = {a/2+R, b/2, 0, 1.0};
//+
Point(15) = {a/2-R, b/2, 0, 1.0};
//+
Point(16) = {a/2, b/2-R, 0, 1.0};
//+
Point(17) = {a/2, b/2+R, 0, 1.0};
//+
Circle(17) = {14, 13, 17};
//+
Circle(18) = {17, 13, 15};
//+
Circle(19) = {15, 13, 16};
//+
Circle(20) = {16, 13, 14};
//+
Curve Loop(26) = {18, 19, 20, 17};
Plane Surface(27) = {26};//+
Curve Loop(17) = {11, -13, 2, -14, 5, -15, 8, -16};
//+
Plane Surface(18) = {26, 17};
//+
Transfinite Curve {17, 18, 19, 20} = nR Using Progression 1;
and then was mesh and implemented to fenics to solve the problem of Periodic homogenization of linear elastic materials with the code following:
from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
a = 1. # unit cell width
b = 1. # unit cell height
R = 0.2 # inclusion radius
vol = a*b # unit cell volume
# we define the unit cell vertices coordinates for later use
vertices = np.array([[0, 0.],
[a, 0.],
[a, b],
[0, b],
[a/2, b/2]])
mesh = Mesh()
with XDMFFile("homo_rec.xdmf") as in_file:
in_file.read(mesh)
mvc_cells = MeshValueCollection("size_t", mesh, mesh.topology().dim())
in_file.read(mvc_cells, "name_to_read")
subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)
plot(subdomains)
plt.show()
# class used to define the periodic boundary map
class PeriodicBoundary(SubDomain):
def __init__(self, vertices, tolerance=DOLFIN_EPS):
""" vertices stores the coordinates of the 4 unit cell corners"""
SubDomain.__init__(self, tolerance)
self.tol = tolerance
self.vv = vertices
self.a1 = self.vv[1,:]-self.vv[0,:] # first vector generating periodicity
self.a2 = self.vv[3,:]-self.vv[0,:] # second vector generating periodicity
# check if UC vertices form indeed a parallelogram
assert np.linalg.norm(self.vv[2, :]-self.vv[3, :] - self.a1) <= self.tol
assert np.linalg.norm(self.vv[2, :]-self.vv[1, :] - self.a2) <= self.tol
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the
# bottom-right or top-left vertices
return bool((near(x[0], self.vv[0,0] + x[1]*self.a2[0]/self.vv[3,1], self.tol) or
near(x[1], self.vv[0,1] + x[0]*self.a1[1]/self.vv[1,0], self.tol)) and
(not ((near(x[0], self.vv[1,0], self.tol) and near(x[1], self.vv[1,1], self.tol)) or
(near(x[0], self.vv[3,0], self.tol) and near(x[1], self.vv[3,1], self.tol)))) and on_boundary)
def map(self, x, y):
if near(x[0], self.vv[2,0], self.tol) and near(x[1], self.vv[2,1], self.tol): # if on top-right corner
y[0] = x[0] - (self.a1[0]+self.a2[0])
y[1] = x[1] - (self.a1[1]+self.a2[1])
elif near(x[0], self.vv[1,0] + x[1]*self.a2[0]/self.vv[2,1], self.tol): # if on right boundary
y[0] = x[0] - self.a1[0]
y[1] = x[1] - self.a1[1]
else: # should be on top boundary
y[0] = x[0] - self.a2[0]
y[1] = x[1] - self.a2[1]
Em = 50e3
num = 0.2
Er = 210e3
nur = 0.3
material_parameters = [(Em, num), (Er, nur)]
nphases = len(material_parameters)
def eps(v):
return sym(grad(v))
def sigma(v, i, Eps):
E, nu = material_parameters[i]
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
return lmbda*tr(eps(v) + Eps)*Identity(2) + 2*mu*(eps(v)+Eps)
Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary())
V = FunctionSpace(mesh, Ve)
v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
dx = Measure('dx')(subdomain_data=subdomains)
Eps = Constant(((0, 0), (0, 0)))
#print(Eps)
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
def macro_strain(i):
"""returns the macroscopic strain for the 3 elementary load cases"""
Eps_Voigt = np.zeros((3,))
Eps_Voigt[i] = 1
return np.array([[Eps_Voigt[0], Eps_Voigt[2]/2.],
[Eps_Voigt[2]/2., Eps_Voigt[1]]])
def stress2Voigt(s):
return as_vector([s[0,0], s[1,1], s[0,1]])
Chom = np.zeros((3, 3))
for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
print("Solving {} case...".format(case))
Eps.assign(Constant(macro_strain(j)))
solve(a == L, w, [], solver_parameters={"linear_solver": "cg"})
(v, lamb) = split(w)
Sigma = np.zeros((3,))
for k in range(3):
Sigma[k] = assemble(sum([stress2Voigt(sigma(v, i, Eps))[k]*dx(i) for i in range(nphases)]))/vol
Chom[j, :] = Sigma
print(np.array_str(Chom, precision=2))
However the results are not correct in compare to those of original example
could you please explain and help me to solve this problem
Best regards
