Hello everyone,
I am new to the FEniCS. I am trying to compute a field (f_a_2) that is based on the line integral of displacement in a particular section of the mesh. The field (f_a_2) also depends on its value from previous iteration. To save the field, I project the field on functionspace and use assign command to use in next iteration. I noticed that during projection there are some negative values, which I supposed that are not possible (due to the use conditional command). Please, help me to understand the reason of these negative values. I attached the minimal code.
Thank you. 
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np
from ufl import replace
import sys, os, shutil, math
from numpy.linalg import eig
from dolfin import *
from numpy import array
import time
from ufl import eq
from ufl import le
from ufl import ge, ne
#Compiler parameter
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
# Define Data
mu = 80770.
lmbda = 121150.
gamma_val = 0.002194
beta_val = 23.649
del_f = 0.01
del_c = 0.0036
# Mesh Data
mesh = Mesh('mesh_2.xml')
ndim = mesh.topology().dim()
remain = CompiledSubDomain("x[1] >= 1.0 - tol",tol = 1E-14)
boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
remain.mark(boundaries, 1)
ds = Measure("ds",subdomain_data=boundaries)
V_u = VectorFunctionSpace(mesh, 'CG', 1)
V_alpha = FunctionSpace(mesh, 'CG', 1)
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)
u_0 = Function(V_u)
jump_0 = Function(V_alpha)
f_a = Function(V_alpha)
f_a_2 = Function(V_alpha)
f_a_0 = Function(V_alpha)
zero_1 = Constant(0.0)
f_a_0 = project(zero_1,V_alpha)
u.rename('displacement','displacement')
f_a_T1 = Constant(2.2e-70)
f_a = project(f_a_T1, V_alpha)
fdf_p = Function(V_alpha)
tol = 1E-10
tol_v = 1e-18
def boundary_D_l(x, on_boundary):
return x[1] <= 0.0 + tol
def boundary_D_o(x, on_boundary):
return x[0] <= 0.0 + tol and x[1] <= 0.0 + tol
def boundary_D_t(x, on_boundary):
return x[1] >= 1.0 - tol
bcl = DirichletBC(V_u.sub(1), Constant(0.0), boundary_D_l, method='pointwise')
bco = DirichletBC(V_u, Constant((0.0 ,0.0)), boundary_D_o, method='pointwise')
# Calculation of positive part of strain
def eps(u):
return sym(grad(u))
def sigma(u):
return 2.0*mu*(eps(u)) + lmbda*tr(eps(u))*Identity(ndim)
def line_integral(field, a,b,n):
assert n > 0
A = np.array([a, b-0.006])
B = np.array([a, b+0.006])
mesh_points = [A + t*(B-A) for t in np.linspace(0, 1, n+1)]
tdim, gdim = 1, len(A)
mesh_l = Mesh()
editor = MeshEditor()
editor.open(mesh_l, 'interval', tdim, gdim)
editor.init_vertices(n+1)
editor.init_cells(n)
for vi, v in enumerate(mesh_points):
editor.add_vertex(vi, v)
for ci in range(n): editor.add_cell(ci, np.array([ci, ci+1], dtype='uintp'))
editor.close()
V_l = FunctionSpace(mesh_l, 'Lagrange',1)
v = interpolate(field, V_l)
return assemble(v*dx)
energy = 0.5 * inner(sigma(u), eps(u)) * dx
E_u = derivative(energy,u,v)
Jd = derivative(E_u, u, du)
num_steps = 25
num_of_large_load_step = 25
total_disp = 0.01
u_R = Expression(('disp_A + disp_app*(n+1-B)'),disp_A = 0.0, disp_app = total_disp/(num_of_large_load_step), n=0., B=0., degree=0)
bct = DirichletBC(V_u.sub(1), u_R, boundary_D_t, method='pointwise')
bc_disp = [bcl, bco, bct]
problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
prm = solver_u.parameters
prm["newton_solver"]["relative_tolerance"] = 5E-3
prm["newton_solver"]["absolute_tolerance"] = 8E-3
prm["newton_solver"]["convergence_criterion"] = "residual"
prm["newton_solver"]["error_on_nonconvergence"] = True
prm["newton_solver"]["linear_solver"] = 'mumps'
prm["newton_solver"]["lu_solver"]["symmetric"] = False
prm["newton_solver"]["maximum_iterations"] = 5000
prm["newton_solver"]["relaxation_parameter"] = 1.0
prm["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True
class MyFunctionExpression(UserExpression):
def __init__(self, n,u_sol, **kwargs):
super().__init__(**kwargs)
self.n = n
self.u_sol = u_sol
def eval(self, value, x):
x_c = x[0]
y_c = x[1]
if y_c <= 0.501 and y_c >= 0.497:
q = line_integral(self.u_sol,x_c,y_c,self.n)
value[0]= q
else:
value[0]= 0.
def value_shape(self):
return (1,)
def alternate_minimization(u):
print('Solution for displacement')
solver_u.solve()
jump = MyFunctionExpression(n=10,u_sol=project(sym(grad(u))[1,1],V_alpha),element = V_alpha.ufl_element())
jump_sign = conditional(ge(jump,0.),1.,0.)
jump_rate = jump-MyFunctionExpression(n=10,u_sol=project(sym(grad(u_0))[1,1],V_alpha),element = V_alpha.ufl_element())
jump_rate_sign = conditional(ge(jump_rate,0.),1.,0.)
jump_star = conditional(ge(jump,del_c),(del_f - del_c)*f_a_0 + del_c,del_c)
alpha_fun = (f_a_0 + gamma_val)*pow((abs(jump_rate)/jump_star),beta_val)
alpha_fun_1 = conditional(lt(alpha_fun,0.),0.,alpha_fun)
f_a_2 = f_a_0 + alpha_fun_1
f_a = conditional(lt(f_a_2, f_a_T1), f_a_T1, f_a_2)
fdf = (2*f_a_T1)/(f_a + f_a_T1)
f_a_0.assign(project(f_a_2,V_alpha))
fdf_p = project(fdf, V_alpha)
file_f_a << (f_a_0,n)
file_fdf << (fdf_p,n)
u_0.assign(u)
return ()
savedir = "results/"
if os.path.isdir(savedir):
shutil.rmtree(savedir)
file_u = File(savedir+"/u.pvd")
file_f_a = File(savedir+"/f_a.pvd")
file_fdf = File(savedir+"/fdf.pvd")
def postprocessing():
file_u << (u,n)
for n in range(num_steps):
u_R.n = n
alternate_minimization(u)
postprocessing()