Hi dokken! Finally I created the below MWE:
Code:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
class boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class disk(SubDomain):
def inside(self, x, on_boundary):
return True
mesh = generate_mesh(Circle(Point(0, 0), 3), 50)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)
ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
n = FacetNormal(mesh)
W1 = FunctionSpace(mesh, "Lagrange", 1)
n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D=Constant(0.0)
bc = DirichletBC(W1, u_D, boundary_markers, 2)
# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx
# Compute solution
u = Function(W1)
solve(a == L, u, bc)
file = File("fsp1.pvd");
file << u;
# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Uy, u))
N = as_matrix([[0.0,0.5],
[0.0,0.0]])
Sn = mu/2 * dot(grad(U) + grad(U).T, (N.T)*-n)
print(assemble(Sn[0]*ds), assemble(Sn[1]*ds))
S=Sn[1]
u = TrialFunction(W1)
v = TestFunction(W1)
a = inner(u,v)*ds
l = inner(S, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(W1)
solve(A, nh.vector(), L)
File("fsp2.pvd") << nh
#print(mesh.geometry().dim())
class expr(UserExpression):
def _init_(self, mesh, marker, f1, **kwargs):
super(expr,self)._init_(kwargs)
self.mesh = mesh
self.marker = marker
self.f1 = f1
def eval_cell(self, value, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
if self.marker[cell] == 1:
value[0] = 0
def value_shape(self):
return()
new_nh = expr(mesh, surface_markers, nh, degree = 1)
W0 = FunctionSpace(mesh, "Lagrange", 1)
new_nh1 = interpolate(new_nh, W0)
gamma1 = Expression(("(-1/D)*(f0 + f1*(1-exp(-beta*(pow(t,2)))))"), degree = 1, D = 0.7, t = new_nh1, f0 = 5, f1 = 2, beta = 0.3)
gamma2 = Constant(0.0) # Dirichlet condition on outer circle
bc2 = DirichletBC(W1, gamma2, boundary_markers, 2)
bcs=[bc2]
k, D = 0.3, 0.7
# Define variational problem
w = TestFunction(W1)
c = TrialFunction(W1)
c = Function(W1)
a = inner(grad(c), grad(w))*dx + (((k/D)*c)*w)*dx
a -= - inner(gamma1, w)*ds
# Compute solution
solve(a == 0, c, bcs)
# Save solution to file in PVD format for Paraview
file = File("fsp3.pvd");
file << c;
The full error message:
Traceback (most recent call last):
File "fspex.py", line 113, in <module>
new_nh1 = interpolate(new_nh, W0)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/interpolation.py", line 71, in interpolate
Pv.interpolate(v._cpp_object)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py", line 383, in interpolate
self._cpp_object.interpolate(u)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/expression.py", line 56, in wrapped_eval_cell
self.user_expression.eval_cell(values, x, cell)
File "fspex.py", line 104, in eval_cell
cell = Cell(self.mesh, ufc_cell.index)
AttributeError: 'expr' object has no attribute 'mesh'
What is the reason for this error and how to fix it?