Ok, so I am providing the .geo script and the full python script below.
gmsh code:
//+
SetFactory("OpenCASCADE");
//+
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("C1", 1) = {3};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {1};
//+
Hide "*";
//+
Show {
Point{3}; Curve{3};
}
//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3};
}
//+
Show "*";
//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; Surface{1}; Surface{2};
}
//+
Curve Loop(5) = {3};
//+
Plane Surface(3) = {5};
//+
Physical Surface("lumen", 6) = {3};//+
Hide "*";
//+
Show {
Point{1}; Point{2}; Point{3}; Curve{1}; Curve{2}; Curve{3}; Surface{3};
}
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("lumen.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)
#print(assemble(1*ds_custom(1)) ,assemble(1*ds_custom(2)))
# Normal vector to the intima
n = FacetNormal(mesh)
# Create mesh and define function space
V = FunctionSpace(mesh, "Lagrange", 1)
# Define boundary condition
G , mu = 10, 0.000004
u_D=Constant(0.0)
bc = DirichletBC(V, u_D, mf, 1)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Ux, Uy, u))
StressSpace = TensorFunctionSpace(mesh, "DG", 0)
S_func = project(S, StressSpace)
S = (mu/2) * (dot(grad(U) + grad(U).T), n)*ds_custom(1)
#print('shear stress =', S)
plot(S)
plot(mesh)
plt.show()
# Save solution to file in PVD format for Paraview
file = File("lumenshear.pvd");
file << u;
The lumen.msh is the msh file of the gmsh script provided.