Tangent matrix and a residual using conditional expression

variational formulation
1

Hey,

I like to print the tangent matrix and the residual for validation. Here is the code:

from dolfin import *
import numpy as np

L = 3
mesh = IntervalMesh(31, 0, L)

Es, Ds, As, phi = Constant(2e8), Constant(12e-3), Constant(pi * (Ds/2)**2), Constant(pi * Ds)
Ec, Dc, Ac = Constant(3e7), Constant(90e-3), Constant(pi * (Dc/2)**2)
s1, s2, s3, s4 = Constant(3e-5), Constant(1e-3), Constant(3e-3), Constant(1e-2)
t1, t2, t4 = Constant(6e3), Constant(12e3), Constant(5e3)

Ue = FiniteElement("CG", mesh.ufl_cell(), 1)
U = FunctionSpace(mesh, MixedElement([Ue, Ue]))
Uc = U.sub(0).collapse()
Us = U.sub(1).collapse()

u = Function(U)
uc, us = split(u)
vc, vs = TestFunctions(U)
w = TrialFunction(U)

def tau(uc, us):
    s = us - uc
    ta = conditional(le(s, s1), (t1/s1) * s,\
        conditional(le(s, s2),(t1 - t2) / (s1 - s2) * s + (s1*t2 - s2*t1) / (s1 - s2),\
            conditional(le(s, s3), t2,\
                conditional(le(s, s4), (t2 - t4) / (s3 - s4) * s + (s3*t4 - s4*t2) / (s3 - s4), t4))))
    return ta

R = (Ec*Ac*uc.dx(0)*vc.dx(0) + Es*As*us.dx(0)*vs.dx(0) + phi * tau(uc, us)*(vs - vc)) * dx
dR = derivative(R, u, w)

uc = interpolate(Expression("(1.0e-4)*x[0]", element = Uc.ufl_element() ), Uc)
us = interpolate(Expression("(2.0e-4)*x[0]", element = Us.ufl_element() ), Us)

for i, c in enumerate(cells(mesh)):
    if i == 0:
        print("vertices element " + str(i))
        for v in vertices(c):
            print (v.point().x())
        A = assemble_local(dR, c)
        b = assemble_local(R, c)
        np.savetxt("K_t element 1.txt", A)
        np.savetxt("R element 1.txt", b)
        print("K_t: \n", A)
        print("R: \n", b)

For these uc, us I get the matrix I am expecting but I receive a vector of zeros instead of the residual I’m expecting.
If I use

uc = interpolate(Expression("(1.0e-4)*x[0]", element = Uc.ufl_element() ), Uc)
us = interpolate(Expression("(2.0e-4)*x[0] + (1.0e-4)", element = Us.ufl_element() ), Us)

which suppose to bring the first element to the second condition of tau I don’t even receive the correct matrix any more (and the residual is still zero).

I guess I’m doing something wrong with the assembly, because if I run the entire code and compare the displacement field I get good agreement.

Any help will be appreciated