apply constrain after solution

Question

Hi, I am new to FEniCS.

I have solved a differential equation. Now i have to apply constraint (1) if sigma+ < sigma- then calculated solution on that node = 0.

How can i apply after calculation?

Thanks in advance.

Comments

What does $\sigma^\pm$ mean? Maybe post a minimal working example of what you're after.

---

\sigma +_ is the positive and negative part of the energy density function. That is based on the displacement.

---

Is $\sigma$ the solution of the differential equation? Or do you have it as an analytical expression or as a dolfin.Function?

---

yes sigma is the solution of the differential equation

Answer 1

Are you familiar with conditional statements in Expressions?

I.e.

expr = Expression("sigma1 < sigma2 ? 0 : sol", sigma1=sigma_plus, sigma2=sigma_minus, sol=solution)
solution = interpolate(expr,solution.function_space())

?

Comments

Hi, Thanks for the help. But after applying this conditional expression I am getting
Error: TypeError: expected default arguments for member variables to be scalars or a scalar GenericFunctions.

I have attached my code:

from dolfin import *
mesh = UnitSquareMesh(4,4)

plot(mesh, title = 'mesh', interactive = True)

def eps(v):
    return sym(grad(v))
def eps_dev(du):
    return dev(eps(du))
def sigma(u):
    return 2.0*mu*eps(u) + lmbda*tr(eps(u))*Identity(2)

V = VectorFunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
du = Function(V)

mu    = 80.77e3
lmbda = 121.15e3
kn = lmbda + mu
EDef = inner(sigma(u),grad(v))*dx


class top(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-10
        return abs(x[1]-1.0) < tol and on_boundary

class bottom(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-10
        return abs(x[1]) < tol and on_boundary
Top = top()
Bottom = bottom()
bclx= DirichletBC(V.sub(0), Constant(0.0), Bottom)
bcly = DirichletBC(V.sub(1), Constant(0.0), Bottom)
bcty = DirichletBC(V.sub(1), Constant(1e-4), Top)
bc_u = [bcly,bclx,  bcty]
solve(lhs(EDef)==rhs(EDef),du,bc_u)
#---------------------
# Definition for energy
#--------------------
def en_dens(u):
    str_ele = 0.5*(grad(u) + grad(u).T)
    IC = tr(str_ele)
    ICC = tr(str_ele * str_ele)
    return (0.5*lmbda*IC**2) + mu*ICC

WW = FunctionSpace(mesh, 'CG',1)
Histold = Function(WW)
HistNew = Function(WW)
# previous step solution
uold = Function(V)
# check enegy at the current step and previous step
expr = Expression('sigma1 < sigma2 ? 0: sol', sigma1 = en_dens(uold), sigma2 = en_dens(du), sol = HistNew)
HistNew = interpolate(expr, HistNew.function_space())

Answer 2

I assume that sigma is the solution that you computed and that sigmin is a given scalar value. Then I would proceed like:

sigarray = sigma.vector().array()  
# make it a numpy array

sigarray[sigarray < sigmin] = 0  
# set all values to zero that are < sigmin

adjsigma = sigma.vector().set_local(sigarray)  
# reassign the values of the dolfin function

Comments

Whenever I have applied

parray = p.vector().array()
parray[parray < 0] = 0
adjparray = p.vector().set_local(parray)

Here, p if temperature solution
I have got error:
parray = p.vector().array()
AttributeError: 'Indexed' object has no attribute 'vector'

---

p should be a dolfin.Function...

---

i have not got what do u mean by dolfin.Function p is the solution vector for temperature.

---

It should be a Function object in fenics. Like you would get it with this example

import dolfin
mesh = dolfin.UnitSquareMesh(10,10)
V = dolfin.FunctionSpace(mesh, 'CG', 1)
p = dolfin.Function(V)

or if you solve a problem via

# Compute solution
...
u = Function(V)
solve(a == L, u, bc)

cf., e.g., http://fenics.readthedocs.io/projects/dolfin/en/latest/demos/demo_poisson.py.html#demo-poisson-equation