Hello guys !
I’m trying to solve a homogenization problem of the geometry attached below. And i don’t know how to apply the same principle of mapping mentioned in here Periodic homogenization of linear elastic materials or how to define my vertices because i don’t understand the principle nor the meaning behind mapping periodic boundaries.
Hi @CHAIkah,
I have posted an explanation of periodic boundary mappings here: Periodic BC Problem in 3D - #6 by conpierce8
Perhaps that will clear up some of your confusion.
Note that for periodic boundary mapping to be successful, you need your mesh to be periodic, i.e. the nodes on the left/right and top/bottom boundaries should be in the same locations.
Hey @conpierce8 !
i tried to combine what you said with this code Periodic homogenization but i still don’t get a perfectly symmetrical effective tensor.
[[ 2.73e+05 1.04e+05 -1.48e-05]
[ 1.04e+05 2.72e+05 -2.63e-05]
[ 7.93e-06 1.41e-05 8.42e+04]]
This is my geometry named 4cells.geo :
Point(1) = {0, 0, 0, 1.0};
Point(2) = {2, 0, 0, 1.0};
Point(3) = {0, 2, 0, 1.0};
Point(4) = {2, 2, 0, 1.0};
Point(5) = {2, 1, 0, 1.0};
Point(6) = {0, 1, 0, 1.0};
Point(7) = {1, 0, 0, 1.0};
Point(8) = {1, 1, 0, 1.0};
Point(9) = {1, 2, 0, 1.0};
//+
Line(1) = {1, 7};
Line(2) = {7, 2};
Line(3) = {6, 8};
Line(4) = {8, 5};
Line(5) = {3, 9};
Line(6) = {9, 4};
Line(7) = {1, 6};
Line(8) = {6, 3};
Line(9) = {7, 8};
Line(10) = {8, 9};
Line(11) = {2, 5};
Line(12) = {5, 4};
//+
Curve Loop(1) = {7, 3, -9, -1};
Plane Surface(1) = {1};
Curve Loop(2) = {2, 11, -4, -9};
Plane Surface(2) = {2};
Curve Loop(3) = {3, 10, -5, -8};
Plane Surface(3) = {3};
Curve Loop(4) = {4, 12, -6, -10};
Plane Surface(4) = {4};
//+
Physical Curve(1) = {7, 8};
Physical Curve(2) = {1, 2};
Physical Curve(3) = {11, 12};
Physical Curve(4) = {6, 5};
//+
Physical Surface(0) = {1};
Physical Surface(1) = {2};
Physical Surface(2) = {3};
Physical Surface(3) = {4};
And this is the fenics code
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Define geometry
fname = "4cells"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
L = 2.
vol = 4.
tol1 =DOLFIN_EPS
class PeriodicBoundary(SubDomain):
def inside(self, y, on_boundary):
#left and front boundary are target domain
return bool((near(y[0],0) and on_boundary) \
or (near(y[1],0,tol1) and on_boundary)) \
and (not(near(y[0],L) or near(y[1],L)))
def map(self,x,y):
#map right and back boundary to target domain
if near(x[0], L) and near(x[1], L):
# Added this condition to deal with the
# case x = y = L
y[0] = x[0] - L
y[1] = x[1] - L
elif near(x[0], L):
y[0] = x[0] - L
y[1] = x[1]
elif near(x[1], L):
y[0] = x[0]
y[1] = x[1] - L
else:
y[0] = -0.5
y[1] = -0.5
#define materials
Em, num = 220e3 , 0.25
Er, nur = 210e3 , 0.3
material_parameters = [(Em, num), (Er, nur),(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)))
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))
lmbda_hom = Chom[0, 1]
mu_hom = Chom[2, 2]
print(Chom[0, 0], lmbda_hom + 2*mu_hom)
E_hom = mu_hom*(3*lmbda_hom + 2*mu_hom)/(lmbda_hom + mu_hom)
nu_hom = lmbda_hom/(lmbda_hom + mu_hom)/2
print("Apparent Young modulus:", E_hom)
print("Apparent Poisson ratio:", nu_hom)
i also noticed that with a model of an 8x8 cells, when i defined the 64 materials for each cell with a Gaussian distribution like this:
E, nu = 210e3 , 0.3
random_E = np.random.normal(E,2,64)
random_nu = np.random.normal(nu,0.03,64)
the code only works if the standard deviation is small (0.03 in this case). and it’s also not symmetrical:
[[-6.46e-03 -3.64e-01 -1.49e+00]
[-1.93e-01 -3.94e-03 1.57e+00]
[ 7.17e-01 -1.94e+00 2.45e-04]]
otherwise it gives me this error : # Solve nonlinear variational problem
Hi @CHAIkah,
I wouldn’t expect you to get a perfectly symmetric effective tensor. As the author notes in the link you provided:
Is there a reason you expect the tensor to be perfectly symmetric?
If you provide your mesh and code for this case, I can be of more assistance. I think the reason your code fails for stdev > 0.03 is because you achieve non-physical values of the Poisson’s ratio. Remember that \nu is physically restricted to -1 \leq \nu \leq 0.5, whereas any value of \nu is theoretically possible from a normal distribution. Below, I sampled 64 points from a normal distribution for various values of stdev, plotted the minimum and maximum values, and repeated 10 times to get a sense of the range of min/max values that could be expected. As you can see, for stdev as low as 0.05, I sometimes obtain values outside the theoretical limits. In addition, you might experience problems any time \nu is close to 0.5, due to the near-incompressibility of the material.
import numpy as np
from matplotlib import pyplot as plt, cm
nu = 0.3
stdev = np.arange(0.01, 0.16, 0.01)
N_domains = 64
N_trials = 10
nu_random = np.zeros((stdev.size, N_domains, N_trials))
plt.figure()
plt.plot([0.01,0.15],[0.5,0.5],'k-',linewidth=1)
for i in range(stdev.size):
s = stdev[i]
for j in range(N_trials):
nu_random[i,:,j] = np.random.normal(nu, s, N_domains)
plt.plot(s*np.ones((N_trials,)), np.amin(nu_random,axis=1)[i], 'o',
color=cm.tab20(i))
plt.plot(s*np.ones((N_trials,)), np.amax(nu_random,axis=1)[i], 'o',
color=cm.tab20(i))
plt.xlabel("Standard Deviation")
plt.ylabel("Min/Max $\\nu$")
plt.annotate("$\\nu_{min}$", [0.01,0.2])
plt.annotate("$\\nu_{max}$", [0.01,0.39])
plt.annotate("Physical Upper Limit", [0.01,0.51])
plt.show()
@conpierce8 i understand but why when i try to solve the same problem with code mentioned here Extract Stiffness tensor with the same gaussian distribution it gives us a symmetrical effective tensor. Isn’t the two models equivalent ?
The two models are equivalent. Do you get a perfectly symmetric effective tensor, or one which is only symmetric within discretization error (i.e. to about 14 digits)? When I run this code on a 2x2 grid, I obtain the following:
1.905555199155007131e+03 7.201990857886308959e+02 -3.044448448593542474e-13
7.201990857886319191e+02 1.905555199155005766e+03 3.235937159125762382e-13
3.836114830940946321e-13 -3.105551175433411946e-13 5.926780566831878332e+02
Comparing the 12 and 21 components, you can see it is not perfectly symmetric; rather, it is only symmetric to 14 digits:
7.201990857886308959e+02
7.201990857886319191e+02
| || |
|<-symmetric->||<->|
not symmetric
No @conpierce8 with the 2×2 grid everything is good
The problem is with grids that are bigger like the 8×8 one that i mentioned
Hi @CHAIkah,
Unless you post your code, I probably can’t diagnose your specific issue. I was able to create a code based on Post #3 that works for an arbitrary-sized grid and produces a symmetric tensor.
from dolfin import *
import meshio
import numpy as np
import matplotlib.pyplot as plt
import gmsh
geo = gmsh.model.geo
#%% Define geometry
Lx, Ly = 1, 1
Nx, Ny = 8, 8
vol = Lx*Ly
gmsh.initialize()
# Create points
for iy in range(Ny+1):
for ix in range(Nx+1):
geo.addPoint(ix*Lx/Nx, iy*Ly/Ny, 0, iy*(Nx+1)+ix+1)
# Create horizontal lines
line_tags = []
bnd_tags = []
for iy in range(Ny+1):
for ix in range(Nx):
tag = iy*Nx+ix+1
geo.addLine(iy*(Nx+1)+ix+1, iy*(Nx+1)+ix+2, tag)
if iy == 0 or iy == Ny:
bnd_tags.append(tag)
else:
line_tags.append(tag)
# Create vertical lines
for iy in range(Ny):
for ix in range(Nx+1):
tag = Nx*(Ny+1)+iy*(Nx+1)+ix+1
geo.addLine(iy*(Nx+1)+ix+1, (iy+1)*(Nx+1)+ix+1, tag)
if ix == 0 or ix == Nx:
bnd_tags.append(tag)
else:
line_tags.append(tag)
# Create curve loops
for iy in range(Ny):
for ix in range(Nx):
bottomLine = iy*Nx + ix + 1
rightLine = Nx*(Ny+1)+iy*(Nx+1)+ix+2
tag = iy*Nx+ix+1
geo.addCurveLoop([bottomLine, rightLine, -bottomLine-Nx, -rightLine+1],
tag)
geo.addPlaneSurface([tag], tag)
geo.addPhysicalGroup(2, [tag], tag)
# Mark boundaries
geo.addPhysicalGroup(1, bnd_tags, 1)
geo.addPhysicalGroup(1, line_tags, 2)
geo.synchronize()
# Constrain periodic boundaries
translation = [1,0,0,0, 0,1,0,Ly, 0,0,1,0, 0,0,0,1]
for ix in range(Nx):
gmsh.model.mesh.setPeriodic(1, [Nx*Ny+ix+1], [ix+1], translation)
translation = [1,0,0,Lx, 0,1,0,0, 0,0,1,0, 0,0,0,1]
for iy in range(Ny):
tag = Nx*(Ny+1)+iy*(Nx+1)+1
gmsh.model.mesh.setPeriodic(1, [tag], [tag+Nx], translation)
# Generate mesh and convert to XDMF
gmsh.model.mesh.generate(1)
gmsh.model.mesh.generate(2)
gmsh.write("grid.msh")
gmsh.finalize()
tmp_msh = meshio.read("grid.msh")
tmp_msh.prune_z_0()
msh = meshio.Mesh(tmp_msh.points, cells={"triangle": tmp_msh.get_cells_type("triangle")})
meshio.write("grid.xdmf", msh)
cellData = tmp_msh.get_cell_data("gmsh:physical", "triangle")
pg_msh = meshio.Mesh(tmp_msh.points,
cells={"triangle": tmp_msh.get_cells_type("triangle")},
cell_data={"bdry": [cellData]},
)
meshio.write("grid_{0}.xdmf".format(cell_type), pg_msh)
#%% Solve homogenization problem
mshfile = XDMFFile("grid.xdmf")
mesh = Mesh()
mshfile.read(mesh)
mshfile.close()
mshfile = XDMFFile("grid_triangle.xdmf")
subdomains = MeshFunction("size_t", mesh, 2)
mshfile.read(subdomains)
mshfile.close()
tol1 =DOLFIN_EPS
class PeriodicBoundary(SubDomain):
def inside(self, y, on_boundary):
#left and front boundary are target domain
return bool((near(y[0],0) and on_boundary) \
or (near(y[1],0,tol1) and on_boundary)) \
and (not(near(y[0],Lx) or near(y[1],Ly)))
def map(self,x,y):
#map right and back boundary to target domain
if near(x[0], Lx) and near(x[1], Ly):
# Added this condition to deal with the
# case x = y = L
y[0] = x[0] - Lx
y[1] = x[1] - Ly
elif near(x[0], Lx):
y[0] = x[0] - Lx
y[1] = x[1]
elif near(x[1], Ly):
y[0] = x[0]
y[1] = x[1] - Ly
else:
y[0] = -0.5
y[1] = -0.5
#define materials
Em, num = 220e3 , 0.25
Er, nur = 210e3 , 0.3
material_parameters = [(Er, nur) for nur in np.random.normal(nur, 0.05, Nx*Ny)]#[(Em, num), (Er, nur),(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)))
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i+1) 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+1) for i in range(nphases)]))/vol
Chom[j, :] = Sigma
print(np.array_str(Chom, precision=2))
lmbda_hom = Chom[0, 1]
mu_hom = Chom[2, 2]
print(Chom[0, 0], lmbda_hom + 2*mu_hom)
E_hom = mu_hom*(3*lmbda_hom + 2*mu_hom)/(lmbda_hom + mu_hom)
nu_hom = lmbda_hom/(lmbda_hom + mu_hom)/2
print("Apparent Young modulus:", E_hom)
print("Apparent Poisson ratio:", nu_hom)
Result:
[[2.82e+05 1.20e+05 3.15e+01]
[1.20e+05 2.82e+05 3.10e+01]
[3.15e+01 3.10e+01 8.13e+04]]
i understand now thank you !
in your code, if we want to extract the strain and stress deformation in each cell for each elementary load cases (uniaxial strain and pure shear), how can we do that ?
I’m not sure how to do this, since I’ve never done it myself. In what form do you wish to extract these quantities? (e.g. as a plot; as the average value of the strain and stress components over the cell; as a list of nodal values of the stress/strain components…)
I wanna extract them over each mesh cell
You can project the stress/strain into an appropriate function space (e.g. DG-1 since your displacements are CG-2) and use stress_fs.dofmap.cell_dofs with stress.vector to obtain the nodal values for a given cell:
stress_fs = TensorFunctionSpace(mesh, "DG", 1)
stress = project(<stress expression>, stress_fs)
stress_vector = stress.vector()
for cell in range(mesh.cells().shape[0]):
nodal_stresses = stress_vector[stress_fs.dofmap().cell_dofs(cell)]
print(nodal_stresses)
Hi I am wondering if fenics can handle aperiodic mesh for periodic boundary condition