Computing wall normal stress and shear stress on mesh

[sarshia]

Hello,

I am trying to compute and visualize the shear stress and normal stress/traction but having issues with the dimension. I wanted to ask if you might know the issue or a better way to approach this.

The computed velocity and pressure profiles are stored in u (3D in V = VectorFunctionSpace(mesh, "CG", 2)) and p (in FunctionSpace(mesh, "CG", 1)), respectively. I am trying to do the following:

n = FacetNormal(mesh)
I = as_tensor([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

sigma_u = ((sym(grad(u))-p*I))
P = TensorFunctionSpace(mesh, "DG", 1)
sigma_u_func = Function(P)
sigma_u_func.interpolate(sigma_u)

sigma_n = sigma_u_func*n
T = VectorFunctionSpace(mesh, "DG", 1)
sigma_n_func = Function(T)
sigma_n_func.interpolate(sigma_n)

I want to be able to visualize the vectors in Paraview so they need to be stored as functions, but I get the following errors for the lines above:

Traceback (most recent call last):
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/code/CylinderStokes.py", line 163, in
sigma_n_func.interpolate(sigma_n)
File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/function.py", line 410, in interpolate
return interpolation.interpolate(
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/interpolation.py", line 143, in interpolate
return Interpolator(
^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/pyadjoint/pyadjoint/tape.py", line 109, in wrapper
return function(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/interpolation.py", line 638, in init
self.callable, arguments = make_interpolator(expr, V, subset, access, bcs=bcs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/interpolation.py", line 818, in make_interpolator
loops.extend(_interpolator(V, tensor, expr, subset, arguments, access, bcs=bcs))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "", line 2, in _interpolator
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/utils.py", line 85, in wrapper
return f(*args, **kwargs)
^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/interpolation.py", line 894, in _interpolator
kernel = compile_expression(cell_set.comm, expr, to_element, V.ufl_element(),
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/PyOP2/pyop2/caching.py", line 300, in wrapper
v = func(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/firedrake/firedrake/interpolation.py", line 1028, in compile_expression
return compile_expression_dual_evaluation(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/driver.py", line 268, in compile_expression_dual_evaluation
evaluation, basis_indices = to_element.dual_evaluation(fn)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/FInAT/finat/tensorfiniteelement.py", line 175, in dual_evaluation
expr = fn(x)
^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/driver.py", line 340, in call
gem_expr, = fem.compile_ufl(self.expression, translation_context, point_sum=False)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/fem.py", line 721, in compile_ufl
result = map_expr_dags(context.translator, expressions)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/ufl/ufl/corealg/map_dag.py", line 98, in map_expr_dags
r = handlersv.ufl_typecode
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/ufl_utils.py", line 162, in _modified_terminal
return self.modified_terminal(o)
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/fem.py", line 365, in modified_terminal
return translate(mt.terminal, mt, self.context)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/local/Cellar/python@3.11/3.11.6_1/Frameworks/Python.framework/Versions/3.11/lib/python3.11/functools.py", line 909, in wrapper
return dispatch(args[0].class)(*args, **kw)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/fem.py", line 441, in translate_reference_normal
return ctx.entity_selector(callback, mt.restriction)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/fem.py", line 98, in entity_selector
return callback(self.entity_ids[0])
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/tsfc/tsfc/fem.py", line 439, in callback
n = ctx.fiat_cell.compute_reference_normal(ctx.integration_dim, facet_i)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/arshiasinghal/Desktop/Rice/Research/firedrake/src/fiat/FIAT/reference_element.py", line 445, in compute_reference_normal
assert facet_dim == self.get_spatial_dimension() - 1
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError

Is this an issue with the polynomial degree or perhaps something else? Please let me know if you would like any additional information, thank you.

[dham]

The proximate cause of your error is that you are attempting to interpolate an expression containing FacetNormal. This won't work because interpolation is implemented cell-wise, so the interpolation code can't possibly know what the facet normal direction is. Indeed, DG1 degrees of freedom won't even have a facet association (all DG DoFs are associated with the cell), so the expression doesn't really make mathematical sense.

I'm slightly struggling to see what quantity it is that you are attempting to represent in the vector space. Do you have the maths of what you are trying to do written up somewhere?

[sarshia]

Hello, I wanted to follow up and ask if you might have an idea based on the clarification below. Thank you again for your time.

[sarshia]

Hello, thank you for your response. As a simpler case, I am trying to calculate grad(u)n where u is the computed (vector) solution and n is the normal vector. This would ultimately correspond with the shear stress. However, in order to write this to the .pvd file to visualize it, it needs to be a function. I thought this might involve interpolation but I see this causes an issue with the cell-wise interpolation. Is this something that might be possible? Please let me know if there is anything else I can clarify.

[rckirby]

I think the application-level problem is well-defined here. For each piece of the boundary, you want to take some quantity that depends on the normal on that bit of the boundary and look at it. But, the code is a bit fussy here since we have symbolic normals and don't have a brute-force way to "iterate over boundary facets, do X with the normal vector"

Is it possible to stick the quantity you want to visualize in some other space like Raviart-Thomas or BDM, then the boundary dofs are normal components?