Hi there!
Here is the problem:
Consider an arbitrary closed domain (not circular). In the MWE, I considered a circular domain just to check if the code works. Consider a point of reference inside the domain. The shear stress on the boundary is nh(x,y). In order to get nh(\theta) from nh(x,y), I need r and \theta. Now since the domain is not circular, r is a function of \theta. In the below MWE, I modified a code suggested by @kamensky to get the changing r. But it returns an error (which I understand but don’t know how to get around it). Here is the MWE followed by the error message:
MWE:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
from sympy import symbols
from sympy import atan2,Abs,cos,sin
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 as np
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
# Allow extrapolation of nh in case a point falls slightly
# outside of the mesh.
nh.set_allow_extrapolation(True)
# Iterate over a bunch of points around the circular
# boundary and evaluate nh at each one.
import math
from numpy import array, zeros
N_theta = 100
dtheta = 2*math.pi/N_theta
nh_theta = zeros(N_theta)
theta = zeros(N_theta)
r = zeros(N_theta)
for i in range(0,N_theta):
theta[i] = i*dtheta
r[i] = ds/theta[i]
x = np.array([r[i]*math.cos(theta[i]), r[i]*math.sin(theta[i])])
nh_theta[i] = nh(x)
plt.plot(theta,nh_theta)
plt.show()
Error message:
r[i] = ds/theta[i]
TypeError: unsupported operand type(s) for /: 'Measure' and 'float'
Any suggestions would be very helpful.