# Shear stress Calculation, converting a FEniCS code to FEniCSx

**URL:** https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803
**Category:** General
**Created:** [November 21, 2022, 6:56am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803 "2022-11-21T06:56:29Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)
#### Post date: [November 21, 2022, 6:56am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/1 "2022-11-21T06:56:29Z")

</div>

My intension is to calculate the shear stresses and normal stresses on a pipe under Poiseuille flow ([Test problem 1: Channel flow](https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code1.html) )

I found a FEniCS code that have done the same thing, but now I’m having some trouble to convert it into the FEniCSx syntax.

The FEniCS code is as follows: (Source [here](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwjWrqrKzb77AhUcElkFHXz3DA4QFnoECAwQAQ&url=https%3A%2F%2Fwww.simula.no%2Ffile%2Fsimulaacdc40pdf%2Fdownload&usg=AOvVaw03V4zPY4SjiwoGI_006xSi&cshid=1669011044955189))

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/e1329dac7cc299890af68dc099327d4fd1e75f5a.png)

And I was able to convert upto the following:

Required imports (including the items necessary for the test problem):

```auto
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista

from dolfinx.fem import Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
from dolfinx.plot import create_vtk_mesh
from ufl import (FacetNormal,FacetArea,as_vector, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym,inner)

```

My attempt to convert the FEniCS code into the FEniCSx syntax:

```auto
def epsilon(u):
    return sym(nabla_grad(u))

def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())
#The above two values will anyway be calculated in the test problem. 

n = FacetNormal(mesh)
T = -sigma(u, p)*n

Tn = inner(T, n) # Inner product: https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals_code.html?highlight=inner%20product#expressing-inner-products
Tt = T - Tn*n

#I'm changing the piecewise constant to, order 1 and order 2
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
scalar = FunctionSpace(mesh, p_cg1)
vector = FunctionSpace(mesh, v_cg2)

v = TestFunction(scalar)
w = TestFunction(vector)

normal_stress = Function(scalar)
shear_stress = Function(vector)

Ln = (1/FacetArea(mesh))* v*Tn * ds
Lt = (1/FacetArea(mesh))* inner(w,Tt) * ds

```

But, then I’m confused on how to properly assemble it and obtain the shear stresses on the boundary (We have a unit square as the domain in the [Test problem 1](https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code1.html))

Can you please help me to proceed after this step.  
Thank you in advance

**A minimum reproducible example for the NS problem is:**  
(The exact code in the [Test problem 1: Channel flow](https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code1.html) )

> ```auto
> from mpi4py import MPI
> from petsc4py import PETSc
> import numpy as np
> import pyvista
> 
> from dolfinx.fem import Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
> from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
> from dolfinx.io import XDMFFile
> from dolfinx.mesh import create_unit_square
> from dolfinx.plot import create_vtk_mesh
> from ufl import (FacetNormal,FacetArea,as_vector, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement,
> div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym,inner)
> 
> mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
> t = 0
> T = 10
> num_steps = 500
> dt = T/num_steps
> 
> v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
> s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
> V = FunctionSpace(mesh, v_cg2)
> Q = FunctionSpace(mesh, s_cg1)
> 
> u = TrialFunction(V)
> v = TestFunction(V)
> p = TrialFunction(Q)
> q = TestFunction(Q)
> 
> def walls(x):
> return np.logical_or(np.isclose(x[1],0), np.isclose(x[1],1))
> wall_dofs = locate_dofs_geometrical(V, walls)
> u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
> bc_noslip = dirichletbc(u_noslip, wall_dofs, V)
> 
> def inflow(x):
> return np.isclose(x[0], 0)
> inflow_dofs = locate_dofs_geometrical(Q, inflow)
> bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)
> 
> def outflow(x):
> return np.isclose(x[0], 1)
> outflow_dofs = locate_dofs_geometrical(Q, outflow)
> bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
> bcu = [bc_noslip]
> bcp = [bc_inflow, bc_outflow]
> 
> u_n = Function(V)
> u_n.name = "u_n"
> U = 0.5 * (u_n + u)
> n = FacetNormal(mesh)
> f = Constant(mesh, PETSc.ScalarType((0,0)))
> k = Constant(mesh, PETSc.ScalarType(dt))
> mu = Constant(mesh, PETSc.ScalarType(1))
> rho = Constant(mesh, PETSc.ScalarType(1))
> 
> # Define strain-rate tensor
> def epsilon(u):
> return sym(nabla_grad(u))
> 
> # Define stress tensor
> def sigma(u, p):
> return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())
> # https://jorgensd.github.io/dolfinx-tutorial/chapter2/navierstokes.html
> 
> # Define the variational problem for the first step
> p_n = Function(Q)
> p_n.name = "p_n"
> F1 = rho*dot((u - u_n) / k, v)*dx
> F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
> F1 += inner(sigma(U, p_n), epsilon(v))*dx
> F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds
> F1 -= dot(f, v)*dx
> a1 = form(lhs(F1))
> L1 = form(rhs(F1))
> 
> A1 = assemble_matrix(a1, bcs=bcu)
> A1.assemble()
> b1 = create_vector(L1)
> 
> # Define variational problem for step 2
> u_ = Function(V)
> a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
> L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
> A2 = assemble_matrix(a2, bcs=bcp)
> A2.assemble()
> b2 = create_vector(L2)
> 
> # Define variational problem for step 3
> p_ = Function(Q)
> a3 = form(dot(u, v)*dx)
> L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
> A3 = assemble_matrix(a3)
> A3.assemble()
> b3 = create_vector(L3)
> 
> # Solver for step 1
> solver1 = PETSc.KSP().create(mesh.comm)
> solver1.setOperators(A1)
> solver1.setType(PETSc.KSP.Type.BCGS)
> pc1 = solver1.getPC()
> pc1.setType(PETSc.PC.Type.HYPRE)
> pc1.setHYPREType("boomeramg")
> 
> # Solver for step 2
> solver2 = PETSc.KSP().create(mesh.comm)
> solver2.setOperators(A2)
> solver2.setType(PETSc.KSP.Type.BCGS)
> pc2 = solver2.getPC()
> pc2.setType(PETSc.PC.Type.HYPRE)
> pc2.setHYPREType("boomeramg")
> 
> # Solver for step 3
> solver3 = PETSc.KSP().create(mesh.comm)
> solver3.setOperators(A3)
> solver3.setType(PETSc.KSP.Type.CG)
> pc3 = solver3.getPC()
> pc3.setType(PETSc.PC.Type.SOR)
> 
> xdmf = XDMFFile(mesh.comm, "poiseuille.xdmf", "w")
> xdmf.write_mesh(mesh)
> xdmf.write_function(u_n, t)
> xdmf.write_function(p_n, t)
> 
> import numpy as np
> def u_exact(x):
> values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
> values[0] = 4*x[1]*(1.0 - x[1])
> return values
> u_ex = Function(V)
> u_ex.interpolate(u_exact)
> 
> L2_error = form(dot(u_ - u_ex, u_ - u_ex)*dx)
> 
> for i in range(num_steps):
> # Update current time step
> t += dt
> 
> # Step 1: Tentative veolcity step
> with b1.localForm() as loc_1:
> loc_1.set(0)
> assemble_vector(b1, L1)
> apply_lifting(b1, [a1], [bcu])
> b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
> set_bc(b1, bcu)
> solver1.solve(b1, u_.vector)
> u_.x.scatter_forward()
>     
> # Step 2: Pressure corrrection step
> with b2.localForm() as loc_2:
> loc_2.set(0)
> assemble_vector(b2, L2)
> apply_lifting(b2, [a2], [bcp])
> b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
> set_bc(b2, bcp)
> solver2.solve(b2, p_.vector)
> p_.x.scatter_forward()
> 
> # Step 3: Velocity correction step
> with b3.localForm() as loc_3:
> loc_3.set(0)
> assemble_vector(b3, L3)
> b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
> solver3.solve(b3, u_.vector)
> u_.x.scatter_forward()
> # Update variable with solution form this time step
> u_n.x.array[:] = u_.x.array[:]
> p_n.x.array[:] = p_.x.array[:]
> 
> # Write solutions to file
> xdmf.write_function(u_n, t)
> xdmf.write_function(p_n, t)
> 
> # Compute error at current time-step
> error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
> error_max = mesh.comm.allreduce(np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX)
> # Print error only every 20th step and at the last step
> # if (i % 20 == 0) or (i == num_steps - 1):
> # print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")
> # Close xmdf file
> xdmf.close()
> 
> ```
> 
> _ **To plot the velocity field and to make sure it works fine:** _
> 
> ```auto
> # pyvista.set_jupyter_backend("pythreejs")
> # pyvista.set_jupyter_backend("ipygany")
> 
> pyvista.set_jupyter_backend("ipyvtklink")
> pyvista.start_xvfb()
> 
> topology, cell_types, geometry = create_vtk_mesh(V)
> values = np.zeros((geometry.shape[0], 3), dtype=np.float64)
> values[:, :len(u_n)] = u_n.x.array.real.reshape((geometry.shape[0], len(u_n)))
> 
> # Create a point cloud of glyphs
> function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
> function_grid["u"] = values
> glyphs = function_grid.glyph(orient="u", factor=0.2)
> 
> # Create a pyvista-grid for the mesh
> grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh, mesh.topology.dim))
> 
> # Create plotter
> plotter = pyvista.Plotter()
> # plotter.add_mesh(grid, style="wireframe", color="k")
> plotter.add_mesh(glyphs)
> plotter.view_xy()
> if not pyvista.OFF_SCREEN:
> plotter.show()
> else:
> pyvista.start_xvfb()
> fig_as_array = plotter.screenshot("glyphs.png")
> 
> ```
> 
> **The code that I have so far, for the shear stress calculation:**
> 
> ```auto
> n = FacetNormal(mesh)
> T = -sigma(u, p)*n
> 
> Tn = inner(T, n) # Inner product: https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals_code.html?highlight=inner%20product#expressing-inner-products
> Tt = T - Tn*n
> 
> #I'm changing the piecewise constant to, order 1 and order 2
> p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
> v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
> scalar = FunctionSpace(mesh, p_cg1)
> vector = FunctionSpace(mesh, v_cg2)
> 
> v = TestFunction(scalar)
> w = TestFunction(vector)
> 
> normal_stress = Function(scalar)
> shear_stress = Function(vector)
> 
> Ln = (1/FacetArea(mesh))* v*Tn * ds
> Lt = (1/FacetArea(mesh))* inner(w,Tt) * ds
> 
> assemble_vector(normal_stress.vector, form(Ln))
> assemble_vector(normal_stress.vector, form(Lt))
> 
> ```

---

<div class="post-metadata">

### Author: ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 21, 2022, 7:35am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/2 "2022-11-21T07:35:20Z")

</div>

> [@charithjeewantha](#):
>
> ```auto
> Ln = (1/FacetArea(mesh))* v*Tn * dx
> Ln = (1/FacetArea(mesh))* inner(w,Tt) * dx
> 
> ```

why are you using `dx` and not `ds` as done in the old code? The assemble function should be replaced with `dolfinx.fem.petsc.assemble_vector(normal_stresses.vector, Ln)`

---

<div class="post-metadata">

### Author: ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)
#### Post date: [November 21, 2022, 7:57am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/3 "2022-11-21T07:57:27Z")

</div>

I’m sorry. That was a mistake from my side. I corrected it with `ds` (As we are using a surface integral)

Thank you for the instruction for the assembly. But when I used it, (` assemble_vector(normal_stress.vector, Ln)`), I’m getting the following error:

> ```auto
> ---------------------------------------------------------------------------
> TypeError Traceback (most recent call last)
> Input In [10], in <cell line: 26>()
> 21 Ln = (1/FacetArea(mesh))* v*Tn * ds
> 22 Ln = (1/FacetArea(mesh))* inner(w,Tt) * ds
> ---> 26 assemble_vector(normal_stress.vector, Ln)
> 
> File /usr/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
> 885 if not args:
> 886 raise TypeError(f'{funcname} requires at least '
> 887 '1 positional argument')
> --> 889 return dispatch(args[0]. __class__ )(*args, **kw)
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:206, in _assemble_vector_vec(b, L, constants, coeffs)
> 190 """Assemble linear form into an existing PETSc vector.
> 191 
> 192 Note:
> (...)
> 203 
> 204 """
> 205 with b.localForm() as b_local:
> --> 206 assemble._assemble_vector_array(b_local.array_w, L, constants, coeffs)
> 207 return b
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/assemble.py:198, in _assemble_vector_array(b, L, constants, coeffs)
> 172 @assemble_vector.register(np.ndarray)
> 173 def _assemble_vector_array(b: np.ndarray, L: FormMetaClass, constants=None, coeffs=None):
> 174 """Assemble linear form into a new Vector.
> 175 
> 176 Args:
> (...)
> 195 
> 196 """
> --> 198 constants = _pack_constants(L) if constants is None else constants
> 199 coeffs = _pack_coefficients(L) if coeffs is None else coeffs
> 200 _cpp.fem.assemble_vector(b, L, constants, coeffs)
> 
> TypeError: pack_constants(): incompatible function arguments. The following argument types are supported:
> 1. (form: dolfinx::fem::Form<double>) -> numpy.ndarray[numpy.float64]
> 2. (e: dolfinx::fem::Expression<double>) -> numpy.ndarray[numpy.float64]
> 3. (form: dolfinx::fem::Form<float>) -> numpy.ndarray[numpy.float32]
> 4. (e: dolfinx::fem::Expression<float>) -> numpy.ndarray[numpy.float32]
> 5. (form: dolfinx::fem::Form<std::complex<double> >) -> numpy.ndarray[numpy.complex128]
> 6. (e: dolfinx::fem::Expression<std::complex<double> >) -> numpy.ndarray[numpy.complex128]
> 7. (form: dolfinx::fem::Form<std::complex<float> >) -> numpy.ndarray[numpy.complex64]
> 8. (e: dolfinx::fem::Expression<std::complex<float> >) -> numpy.ndarray[numpy.complex64]
> 
> Invoked with: Form([Integral(Product(Division(IntValue(1), FacetArea(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0))), Inner(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2)), 0, None), Sum(ComponentTensor(Product(IntValue(-1), Indexed(ComponentTensor(Product(Indexed(FacetNormal(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0)), MultiIndex((Index(90),))), Conj(Inner(FacetNormal(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0)), ComponentTensor(IndexSum(Product(Indexed(FacetNormal(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0)), MultiIndex((Index(89),))), Indexed(ComponentTensor(Product(IntValue(-1), Indexed(Sum(ComponentTensor(Product(IntValue(-1), Indexed(ComponentTensor(Product(Indexed(Identity(2), MultiIndex((Index(82), Index(83)))), Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None)), MultiIndex((Index(82), Index(83)))), MultiIndex((Index(84), Index(85))))), MultiIndex((Index(84), Index(85)))), ComponentTensor(Product(Indexed(Sym(NablaGrad(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2)), 1, None))), MultiIndex((Index(80), Index(81)))), Product(IntValue(2), Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 2))), MultiIndex((Index(80), Index(81))))), MultiIndex((Index(86), Index(87))))), MultiIndex((Index(86), Index(87)))), MultiIndex((Index(88), Index(89))))), MultiIndex((Index(89),))), MultiIndex((Index(88),)))))), MultiIndex((Index(90),))), MultiIndex((Index(91),)))), MultiIndex((Index(91),))), ComponentTensor(IndexSum(Product(Indexed(FacetNormal(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0)), MultiIndex((Index(89),))), Indexed(ComponentTensor(Product(IntValue(-1), Indexed(Sum(ComponentTensor(Product(IntValue(-1), Indexed(ComponentTensor(Product(Indexed(Identity(2), MultiIndex((Index(82), Index(83)))), Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None)), MultiIndex((Index(82), Index(83)))), MultiIndex((Index(84), Index(85))))), MultiIndex((Index(84), Index(85)))), ComponentTensor(Product(Indexed(Sym(NablaGrad(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2)), 1, None))), MultiIndex((Index(80), Index(81)))), Product(IntValue(2), Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 2))), MultiIndex((Index(80), Index(81))))), MultiIndex((Index(86), Index(87))))), MultiIndex((Index(86), Index(87)))), MultiIndex((Index(88), Index(89))))), MultiIndex((Index(89),))), MultiIndex((Index(88),)))))), 'exterior_facet', Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), 'everywhere', {}, None)])
> 
> Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
> <pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
> conversions are optional and require extra headers to be included
> when compiling your pybind11 module.
> 
> ```

---

<div class="post-metadata">

### Author: ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 21, 2022, 8:09am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/4 "2022-11-21T08:09:10Z")

</div>

I typed this on my phone, and hence missed that you need to wrap `Ln`

> [@dokken](#):
>
> dolfinx.fem.petsc.assemble\_vector(normal\_stresses.vector, Ln)

should be

```python
dolfinx.fem.petsc.assemble_vector(normal_stresses.vector, dolfinx.fem.form(Ln))

```

---

<div class="post-metadata">

### Author: ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)
#### Post date: [November 21, 2022, 8:16am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/5 "2022-11-21T08:16:41Z")

</div>

mh… Sorry for the trouble but now I’m getting the following:

> ```auto
> ---------------------------------------------------------------------------
> UFLException Traceback (most recent call last)
> Input In [17], in <cell line: 1>()
> ----> 1 assemble_vector(normal_stress.vector, form(Ln))
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:166, in form(form, dtype, form_compiler_params, jit_params)
> 163 return list(map(lambda sub_form: _create_form(sub_form), form))
> 164 return form
> --> 166 return _create_form(form)
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:161, in form.<locals>._create_form(form)
> 158 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
> 159 return form argument"""
> 160 if isinstance(form, ufl.Form):
> --> 161 return _form(form)
> 162 elif isinstance(form, collections.abc.Iterable):
> 163 return list(map(lambda sub_form: _create_form(sub_form), form))
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py:135, in form.<locals>._form(form)
> 132 if mesh is None:
> 133 raise RuntimeError("Expecting to find a Mesh in the form.")
> --> 135 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
> 136 form_compiler_params=form_compiler_params,
> 137 jit_params=jit_params)
> 139 # For each argument in form extract its function space
> 140 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
> 51 @functools.wraps(local_jit)
> 52 def mpi_jit(comm, *args, **kwargs):
> 53 
> 54 # Just call JIT compiler when running in serial
> 55 if comm.size == 1:
> ---> 56 return local_jit(*args, **kwargs)
> 58 # Default status (0 == ok, 1 == fail)
> 59 status = 0
> 
> File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_params, jit_params)
> 202 # Switch on type and compile, returning cffi object
> 203 if isinstance(ufl_object, ufl.Form):
> --> 204 r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
> 205 elif isinstance(ufl_object, ufl.FiniteElementBase):
> 206 r = ffcx.codegeneration.jit.compile_elements([ufl_object], parameters=p_ffcx, **p_jit)
> 
> File /usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py:147, in compile_forms(forms, parameters, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
> 143 p = ffcx.parameters.get_parameters(parameters)
> 145 # Get a signature for these forms
> 146 module_name = 'libffcx_forms_' + \
> --> 147 ffcx.naming.compute_signature(forms, _compute_parameter_signature(p)
> 148 + str(cffi_extra_compile_args) + str(cffi_debug))
> 150 form_names = [ffcx.naming.form_name(form, i, module_name) for i, form in enumerate(forms)]
> 152 if cache_dir is not None:
> 
> File /usr/local/lib/python3.10/dist-packages/ffcx/naming.py:32, in compute_signature(ufl_objects, tag)
> 30 if isinstance(ufl_object, ufl.Form):
> 31 kind = "form"
> ---> 32 object_signature += ufl_object.signature()
> 33 elif isinstance(ufl_object, ufl.FiniteElementBase):
> 34 object_signature += repr(convert_element(ufl_object))
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/form.py:252, in Form.signature(self)
> 250 "Signature for use with jit cache (independent of incidental numbering of indices etc.)"
> 251 if self._signature is None:
> --> 252 self._compute_signature()
> 253 return self._signature
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/form.py:506, in Form._compute_signature(self)
> 503 def _compute_signature(self):
> 504 from ufl.algorithms.signature import compute_form_signature
> 505 self._signature = compute_form_signature(self,
> --> 506 self._compute_renumbering())
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/form.py:469, in Form._compute_renumbering(self)
> 466 def _compute_renumbering(self):
> 467 # Include integration domains and coefficients in renumbering
> 468 dn = self.domain_numbering()
> --> 469 cn = self.coefficient_numbering()
> 470 cnstn = self.constant_numbering()
> 471 renumbering = {}
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/form.py:238, in Form.coefficient_numbering(self)
> 235 """Return a contiguous numbering of coefficients in a mapping
> 236 ``{coefficient:number}``."""
> 237 if self._coefficient_numbering is None:
> --> 238 self._analyze_form_arguments()
> 239 return self._coefficient_numbering
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/form.py:456, in Form._analyze_form_arguments(self)
> 454 "Analyze which Argument and Coefficient objects can be found in the form."
> 455 from ufl.algorithms.analysis import extract_arguments_and_coefficients
> --> 456 arguments, coefficients = extract_arguments_and_coefficients(self)
> 458 # Define canonical numbering of arguments and coefficients
> 459 self._arguments = tuple(
> 460 sorted(set(arguments), key=lambda x: x.number()))
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/analysis.py:127, in extract_arguments_and_coefficients(a)
> 122 if len(bfnp) != len(set(bfnp.values())):
> 123 msg = """\
> 124 Found different Arguments with same number and part.
> 125 Did you combine test or trial functions from different spaces?
> 126 The Arguments found are:\n%s""" % "\n".join(" %s" % f for f in arguments)
> --> 127 error(msg)
> 129 # Build count: instance mappings, should be one to one
> 130 fcounts = dict((f, f.count()) for f in coefficients)
> 
> File /usr/local/lib/python3.10/dist-packages/ufl/log.py:158, in Logger.error(self, *message)
> 156 "Write error message and raise an exception."
> 157 self._log.error(*message)
> --> 158 raise self._exception_type(self._format_raw(*message))
> 
> UFLException: Found different Arguments with same number and part.
> Did you combine test or trial functions from different spaces?
> The Arguments found are:
> v_0
> v_1
> v_1
> 
> ```

Please note that I used:

```auto
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
scalar = FunctionSpace(mesh, p_cg1)
vector = FunctionSpace(mesh, v_cg2)

```

instead of the zero order that was in the FEniCS code.

---

<div class="post-metadata">

### Author: ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 21, 2022, 8:39am UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/6 "2022-11-21T08:39:00Z")

</div>

Please make a minimal reproducible code. I cannot spend time copy-pasting snippets of code from your previous posts guessing how you have added the changes. Please use a built-in mesh and make sure that one can run the code with DOLFINx by simply copy pasting it.

---

<div class="post-metadata">

### Author: ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)
#### Post date: [November 21, 2022, 7:35pm UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/7 "2022-11-21T19:35:51Z")

</div>

Sure. I edited my post with an executable example for the Navier Stokes problem.  
And I was trying to calculate the shear stresses for that vector field

---

<div class="post-metadata">

### Author: ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [November 21, 2022, 8:44pm UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/8 "2022-11-21T20:44:56Z")

</div>

> [@charithjeewantha](#):
>
> ```auto
> Ln = (1/FacetArea(mesh))* v*Tn * ds
> Lt = (1/FacetArea(mesh))* inner(w,Tt) * ds
> 
> assemble_vector(normal_stress.vector, form(Ln))
> assemble_vector(normal_stress.vector, form(Lt))
> 
> ```

`v` is part of the scalar space, while `w` is part of the vector function space.  
Are you sure that you want to call

> [@charithjeewantha](#):
>
> `assemble_vector(normal_stress.vector, form(Lt))`

There is also so much code here (You need to remove parts that are not needed to reproduce the error), it is really hard to see whats going on, your code

> [@charithjeewantha](#):
>
> ```auto
> n = FacetNormal(mesh)
> T = -sigma(u, p)*n
> 
> Tn = inner(T, n) # Inner product: https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals_code.html?highlight=inner%20product#expressing-inner-products
> Tt = T - Tn*n
> 
> ```

This shouldnt be `u` and `p`, but `u_` and `p_` as far as I can tell.

---

<div class="post-metadata">

### Author: ![charithjeewantha](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/charithjeewantha/32/3948_2.png) [@charithjeewantha](https://fenicsproject.discourse.group/u/charithjeewantha)
#### Post date: [November 21, 2022, 8:53pm UTC](https://fenicsproject.discourse.group/t/shear-stress-calculation-converting-a-fenics-code-to-fenicsx/9803/9 "2022-11-21T20:53:29Z")

</div>

@dokken Thank You!  
yes it shouldn’t be `u` and `p`.  
It should be `u_n` and `p_n`.