Project kappa into a suitable function space, and visualize it with the plot command or with paraview
Thank you very much, one last question. So, when I want to use kappa in the variational formulation, should I use the kappa projection or with the kappa I already calculated in my previous code?
You should use the UserExpression, as it is evaluated at quadrature points rather than at degrees of freedom.
Hi
I am trying to solve a problem with different materials.
The code i am using is:
from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
#使用网格函数定义子域
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5 - tol
materials = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
#这会将属于 Ω0 的每个单元上的网格函数材质值设置为 0，将属于 Ω1 的所有单元上的网格函数材质值设置为 1
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
plot(materials)
#File(’materials.xml.gz’) << materials
k_0 = 1
k_1 = 0.01
#使用网格函数定义可变系数K
class K(UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
def values_shape(self):
return (1,)
Kappa = K(materials, k_0, k_1, degree=0)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = Kappa*dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Plot solution and mesh
plot(u)
plot(mesh)
However, I always get this Error:
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
Solving linear variational problem.
RuntimeError Traceback (most recent call last)
Cell In[24], line 68
66 # Compute solution
67 u = Function(V)
---> 68 solve(a == L, u, bc)
70 # Plot solution and mesh
71 plot(u)
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/solving.py:220, in solve(*args, **kwargs)
217 # Call variational problem solver if we get an equation (but not a
218 # tolerance)
219 elif isinstance(args[0], ufl.classes.Equation):
--> 220 _solve_varproblem(*args, **kwargs)
222 # Default case, just call the wrapped C++ solve function
223 else:
224 if kwargs:
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/solving.py:247, in _solve_varproblem(*args, **kwargs)
245 solver = LinearVariationalSolver(problem)
246 solver.parameters.update(solver_parameters)
--> 247 solver.solve()
249 # Solve nonlinear variational problem
250 else:
251
252 # Create Jacobian if missing
253 if J is None:
Here is a corrected version of your code that runs, please study it carefully:
from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
#使用网格函数定义子域
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5 - tol
materials = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
#这会将属于 Ω0 的每个单元上的网格函数材质值设置为 0，将属于 Ω1 的所有单元上的网格函数材质值设置为 1
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
plot(materials)
#File(’materials.xml.gz’) << materials
k_0 = 1
k_1 = 0.01
#使用网格函数定义可变系数K
class K(UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
Kappa = K(materials, k_0, k_1, degree=0)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = Kappa*dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Plot solution and mesh
plot(u)
plot(mesh)
It works now thank you very much.