But I did define c as c = TrialFunction(W)
. Below is the entire code along with the gmsh script of the geometry.
Python script:
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
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
import meshio
msh = meshio.read("intimapdgf.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
# Normal vector to the intima
n = FacetNormal(mesh)
# Create mesh and define function space
W = FunctionSpace(mesh, "Lagrange", 1)
# Define boundary condition and neumann boundary functions
D, k, t = 0.7, 0.3, 3
f0, ts = 5, 2
gamma1 = (-1/D)*Constant(f0*t/(ts + t)) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0) # Dirichlet condition on outer circle
bc2 = DirichletBC(W, gamma2, mf, 2)
bcs=[bc2]
# Define variational problem
w = TestFunction(W)
c = TrialFunction(W)
c = Function(W)
a = inner(grad(c), grad(w))*dx
a -= -((((k/D)*c)*w)*dx) - inner(gamma1, w)*ds_custom(1)
# Compute solution
solve(a == 0, c, bcs)
# Save solution to file in PVD format for Paraview
file = File("intimapdgf.pvd");
file << c;
gmsh script:
//+
SetFactory("OpenCASCADE");
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Point(3) = {-1.6, 0.3, 0, 1.0};
//+
Point(4) = {-1.4, -0.7, 0, 1.0};
//+
Point(5) = {-1.1, -1.3, 0, 1.0};
//+
Point(6) = {-0.3, -1.7, 0, 1.0};
//+
Point(7) = {0.8, -2.1, 0, 1.0};
//+
Point(8) = {1.7, -1.3, 0, 1.0};
//+
Point(9) = {1.7, -0.4, 0, 1.0};
//+
Point(10) = {1.9, 0.5, 0, 1.0};
//+
Point(11) = {1.6, 1.2, 0, 1.0};
//+
Point(12) = {0.9, 1.5, 0, 1.0};
//+
Point(13) = {0, 1.5, 0, 1.0};
//+
Point(14) = {-0.6, 0.9, 0, 1.0};
//+
Spline(3) = {3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 3};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {1};
//+
Hide "*";
//+
Show {
Point{3}; Curve{3};
}
//+
Curve Loop(1) = {2};
//+
Curve Loop(2) = {3};
//+
Plane Surface(1) = {1, 2};
//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3};
}
//+
Curve Loop(3) = {1};
//+
Curve Loop(4) = {2};
//+
Plane Surface(2) = {3, 4};
//+
Show "*";
//+
Physical Surface("intima", 4) = {1};
//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; Surface{1}; Surface{2};
}
//+
Physical Curve("C1", 1) = {3};