# Why is the normal stress non-zero at boundary where I impose zero normal stress?
**URL:** https://fenicsproject.discourse.group/t/why-is-the-normal-stress-non-zero-at-boundary-where-i-impose-zero-normal-stress/5193
**Category:** variational formulation
**Created:** [March 1, 2021, 4:44pm UTC](https://fenicsproject.discourse.group/t/why-is-the-normal-stress-non-zero-at-boundary-where-i-impose-zero-normal-stress/5193 "2021-03-01T16:44:57Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@kamilova-maths](https://fenicsproject.discourse.group/u/kamilova-maths)
#### Post date: [March 1, 2021, 4:44pm UTC](https://fenicsproject.discourse.group/t/why-is-the-normal-stress-non-zero-at-boundary-where-i-impose-zero-normal-stress/5193/1 "2021-03-01T16:44:57Z")
Hello,
I am solving Stokes flow for a rectangular domain. The top of the domain is split into two halves, where the boundary condition on the right half should be zero normal stress. Even though my solution looks ok when I plot it, the normal stress I calculate through that boundary during post-processing is not zero. How can this be?
This is my minimal working example
```auto
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Define mesh and geometry - We solve for half of the domain we need, and impose symmetry
mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)
n = FacetNormal(mesh)
# Define Taylor--Hood function space W
V = VectorElement("CG", triangle, 2)
Q = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))
# Define Function and TestFunction(s)
w = Function(W)
(u, p) = split(w)
(v, q) = split(TestFunction(W))
# Define the viscosity and bcs
mu = Constant(1.0)
# mu = Expression('exp(-a*pow(x[0],2))', degree=2, a=10)
u_in = Constant(-2.0)
u_c = Constant(-1.0)
inflow = 'near(x[1], 1.0) && x[0]<=0.5'
wall = 'near(x[0], 1)'
centre = 'near(x[0], 0.0)'
outflow = 'near(x[1], 0.0)'
bcu_inflow = DirichletBC(W.sub(0), (0.0, u_in), inflow)
bcu_wall = DirichletBC(W.sub(0), (0.0, u_c), wall)
bcu_outflow = DirichletBC(W.sub(0), (0.0, u_c), outflow)
bcu_symmetry = DirichletBC(W.sub(0).sub(0), Constant(0.0), centre)
bcs = [bcu_inflow, bcu_wall, bcu_outflow, bcu_symmetry]
# Define stress tensor
x = SpatialCoordinate(mesh)
def epsilon(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1), 0],
[v[1].dx(0), v[1].dx(1), 0],
[0, 0, 0]]))
# stress tensor
def sigma(v, p):
return 2*mu*epsilon(v)-Id(p)
def Id(p):
return as_tensor([[p, 0, 0],
[0, p, 0],
[0, 0, p]])
def cond(v):
return sym(as_tensor([[v[0].dx(0), v[0].dx(1)],
[v[1].dx(0), v[1].dx(1)]]))
def sigmabc(v, p):
return 2*mu*cond(v) - p*Identity(2)
# Define the variational form
f = Constant((0, -1))
colors = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
colors.set_all(0) # default to zero
# We match the colours to the defined sketch in the Fenics chapter
CompiledSubDomain("near(x[0], 0.0)").mark(colors, 4)
CompiledSubDomain("near(x[1], 1.0) && x[0]<=0.5").mark(colors, 0)
CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.5").mark(colors, 1)
CompiledSubDomain("near(x[0], 1.0)").mark(colors, 2) # wall
CompiledSubDomain("near(x[1], 0.0)").mark(colors, 3) # outflow
# Create the measure
ds = Measure("ds", domain=mesh, subdomain_data=colors)
a1 = (inner(sigma(u, p), epsilon(v))) * dx
a2 = (- div(u) * q - dot(f, v)) * dx
b1 = - dot(dot(sigmabc(u, p), v), n) * ds(0)
b3 = - dot(dot(sigmabc(u, p), v), n) * ds(2)
b4 = - dot(dot(sigmabc(u, p), v), n) * ds(3)
F = a1 + a2
# Solve problem
solve(F == 0, w, bcs)
# Plot solutions
(u, p) = w.split()
# We check stress at the end
stress = (dot(as_vector([sigmabc(u, p)[0, 0], sigmabc(u, p)[1, 1]]), n) ) * ds(1) # stress calculation at ds(1) - this normal stress should be zero!
flux = dot(u, n) * dot(u, n) * ds(1)
total_stress_new = assemble(stress)
total_flux_new = assemble(flux)
print("Total stress on boundary ds1 is", total_stress_new)
print("Total flux on boundary ds1 is", total_flux_new)
```
Does anyone know why my normal stress is not zero? Or maybe a better way to impose zero normal stress at ds(1)?
---
### Author:  [@kamilova-maths](https://fenicsproject.discourse.group/u/kamilova-maths)
#### Post date: [March 1, 2021, 5:35pm UTC](https://fenicsproject.discourse.group/t/why-is-the-normal-stress-non-zero-at-boundary-where-i-impose-zero-normal-stress/5193/2 "2021-03-01T17:35:50Z")
Furthermore, I can calculate the stress explicitly (as I am probably doing it wrong the way it is in the code), as
```auto
stress = (sigmabc(u, p)[0, 0]*n[0] + sigmabc(u, p)[1, 1]*n[1])*ds(1)
```
which is also non-zero. But
```auto
stress = (sigmabc(u, p)[0, 0]*n[0] )*ds(1)
```
is zero. Now I’m really confused… The normal vector at the boundary ds(1) should be as\_vector(n[0], 0), it’s just a vector pointing up, right? Then why does adding the component multiplied by n[1] make it non-zero? I know that for this particular example it doesn’t matter and I can just write it times n[0] as above, but my goal is to solve on a domain where the normal to the boundary does not align with the components of the velocity vectors (e.g. an inclined boundary on a Polygon domain), and so I can’t just make the normal vector (1,0). Any advice would be greatly appreciated. I have been stuck on this for weeks and have been unable to make progress on my actual problem, with the Polygon boundary.
---
### Author:  [@Francesca\_Ballatore](https://fenicsproject.discourse.group/u/Francesca_Ballatore)
#### Post date: [March 2, 2024, 12:34pm UTC](https://fenicsproject.discourse.group/t/why-is-the-normal-stress-non-zero-at-boundary-where-i-impose-zero-normal-stress/5193/3 "2024-03-02T12:34:30Z")
> [@kamilova-maths](#):
>
> `n = FacetNormal(mesh)`
Have you ever solved? I have the same problem!
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [March 19, 2024, 7:53pm UTC](https://fenicsproject.discourse.group/t/why-is-the-normal-stress-non-zero-at-boundary-where-i-impose-zero-normal-stress/5193/4 "2024-03-19T19:53:48Z")
The normal vector at ds(1) should be (0, 1), not (1, 0), as you are looking at the top boundary
> [@kamilova-maths](#):
>
> `mesh = RectangleMesh(Point(0, 0), Point(1, 1), 100, 100)`
> [@kamilova-maths](#):
>
> `CompiledSubDomain("near(x[1], 1.0) && x[0]>=0.5").mark(colors, 1)`
So it is expected that the first term is 0, and that the second term gives a non-zero contribution