Hi everyone, I am currently stuck on an example problem where I have found the pressure and velocity field around an obstacle immersed in a fluid domain by solving the Navier Stokes equation in the considered domain. Further, I compute the stress matrix and use that to compute the traction on each element. The example proceeds to utilize the traction on each element to compute the lift and the drag on the immersed body, however, I diverge at this point and want to find out the traction on each element of the mesh in the form of a NumPy array. I am not able to do this as a component of the traction is returning an indexed element instead of a scalar field which I am not able to proceed with. Let u and p be the velocity and the pressure fields in the given domain which have already been correctly computed.
A minimal example of the problem in the code is as follow:
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (on_boundary and (-x_shape < x[0] < x_shape) and (-y_shape < x[1] < y_shape))
domain = Rectangle(Point(-5.0, -5.0), Point(5.0,5.0)) - Circle(Point(0, 0),2)
mesh = generate_mesh(domain, 40)
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace (mesh, 'P', 1)
n = FacetNormal(mesh)
obstacle = Obstacle()
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
obstacle.mark(boundaries, 5)
ds = Measure('ds', subdomain_data=boundaries)
gamma_shape = ds(5)
eps = 0.5*(nabla_grad(u) + nabla_grad(u).T)
sigma = 2.0*eps - p*Identity(len(u))
traction = dot(sigma, n)
fx = traction[0]
fy = traction[1]
These (fx, fy) computed are of the Indexed type and hence are not scalar fields, I want to store their values in a NumPy array for further post-processing. Is there any way to do this? I would be grateful for any help in this regard.