Dear all,
I’m trying to make a fenics code to solve the thermal problem of 2D two-phases material with periodicboundary condition. I want to make a boundary condition that the temperature T = 1 is applied on the top-side of REV, the others sides are zero temperature and I want to calculate the flux at macro. The code is given bellow:
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 = sqrt(3) # unit cell height
R = 0.2 # inclusion radius
vol = a*b # unit cell volume
Ro1 = pi*R*R/vol
# we define the unit cell vertices coordinates for later use
vertices = np.array([[0, 0.],
[a, 0.],
[a, b],
[0, b]])
mesh = Mesh()
with XDMFFile("Mesh/t3.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 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], 0, self.tol) or
near(x[1], 0, self.tol)) and
(not ((near(x[0], a, self.tol) and near(x[1], 0, self.tol)) or
(near(x[0], 0, self.tol) and near(x[1], b, self.tol)))) and on_boundary)
def map(self, x, y):
if near(x[0], a, self.tol) and near(x[1], b, self.tol): # if on top-right corner
y[0] = x[0] - a
y[1] = x[1] - b
elif near(x[0], a, self.tol): # if on right boundary
y[0] = x[0] - a
y[1] = x[1]
else: # should be on top boundary
y[0] = x[0]
y[1] = x[1] - b
km = 1
kr = 0.5
material_parameters = [km, kr]
nphases = len(material_parameters)
def Flux(u, i):
k = material_parameters[i]
return k*grad(u)
km = 1
kr = 0.5
material_parameters = [km, kr]
nphases = len(material_parameters)
def Flux(u, i):
k = material_parameters[i]
return k*grad(u)
# Define the boundary conditions
V = FunctionSpace(mesh, "P", 1, constrained_domain=PeriodicBoundary(vertices))
T = 1
bc_top = DirichletBC(V, T, PeriodicBoundary(vertices), "near(x[1], b, tol)")
bc_other = DirichletBC(V, 0, PeriodicBoundary(vertices), "on_boundary && !near(x[1], b, tol)")
# Define the variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
A = sum([inner(Flux(u, i), grad(v))*dx(i) for i in range(nphases)])
L = f*v*dx
# Apply the boundary conditions
bcs = [bc_top, bc_other]
# Solve the variational problem
u = Function(V)
solve(A == L, u, bcs)
flux = sum([Flux(u, i)*dx(i) for i in range(nphases)])
print(flux)
when I run this code, there is an error like that:
RuntimeError Traceback (most recent call last)
<ipython-input-6-c7b32d48174a> in <module>
2 V = FunctionSpace(mesh, "P", 1, constrained_domain=PeriodicBoundary(vertices))
3 T = 1
----> 4 bc_top = DirichletBC(V, T, PeriodicBoundary(vertices), "near(x[1], b, tol)")
5 bc_other = DirichletBC(V, 0, PeriodicBoundary(vertices), "on_boundary && !near(x[1], b, tol)")
6
/usr/local/lib/python3.6/dist-packages/dolfin/fem/dirichletbc.py in __init__(self, *args, **kwargs)
129 raise RuntimeError("Invalid keyword arguments", kwargs)
130
--> 131 super().__init__(*args)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to create Dirichlet boundary condition.
*** Reason: unknown method ("near(x[1], b, tol)").
*** Where: This error was encountered inside DirichletBC.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
Could you help me to solve this problem please. Thank you very much!