Hi, I’m quite noob in Fenics 2019.1.0 and I’m having serious problems in the implementation of a Picard iteration method for a nonlinear coupled mixed-primal problem that goes as follow:

Find (\sigma,u,\phi) \in \mathbb{X} \times \mathbf{M} \times H such that

where k is a constant vector, f is a vector function, g is a scalar function.

My idea is to iterate between the first two equations and the third one. That is, start with a given \phi_N, solve the first two equations obtaining \sigma_N and u_N. Then, replace u_N and \phi_N in the third equation (\phi_N just inside \gamma) and solve it for a new \phi_N. Then, repeat until some convergence criteria is met.

It seems that the solutions is not updating properly because the plots are empty.

The importants part of my unworking code are

```
# Load mesh from file
mesh = Mesh() #malla
hdf = HDF5File(mesh.mpi_comm(), "mallas/cuadrado.h5", "r")
hdf.read(mesh, "/mesh", False)
# Referencias de los dominios
domains = MeshFunction("size_t", mesh, dim) #CellFunction("size_t", mesh) #dominio
hdf.read(domains, "/domains")
#Referencias de los lados
ref_sides = MeshFunction("size_t", mesh, dim-1) #FacetFunction("size_t", mesh) #caras
hdf.read(ref_sides, "/bc_markers")
ds = Measure('ds', domain=mesh, subdomain_data=ref_sides)
n = FacetNormal(mesh)
# Problem data
def mu(s):
return pow(1 - c*s,-2)
def gamma(s):
return c*s*pow(1-c*s,2)
def theta(s):
return m1+m2*pow(1+s*s,0.5*m3-1)
#Auxiliar function
def norma(x):
return sqrt(pow(x[0],2)+pow(x[1],2))
def deviator(A):
return A - 0.5*tr(A)*Identity(2)
# Define function spaces
X = VectorElement("RT",mesh.ufl_cell(), degree)
M = VectorElement("DG",mesh.ufl_cell(), degree)
H = FiniteElement("Lagrange",mesh.ufl_cell(),degree + 1)
#Mixed spaces
WS = FunctionSpace(mesh,MixedElement([X,M]))
WStilde = FunctionSpace(mesh,H)
# Define trial and test functions
(sigma,u) = TrialFunctions(WS)
su = Function(WS)
phi = TrialFunction(WStilde)
(tau,v) = TestFunctions(WS)
psi = TestFunction(WStilde)
#Define iteration (previous step) functions
suN = Function(WS)
sigmaN,uN = split(suN)
phiN = Function(WStilde)
#Null boundary conditions for phi
bcphi = DirichletBC(WStilde,Constant(0), ref_sides, 0)
bcs = [bcphi]
for i in range(10):
#S(phi) = (sigma,u)
lhss1 = inner(deviator(sigma)/mu(phiN),deviator(tau))*dx + inner(u,div(tau))*dx + inner(v,div(sigma))*dx
rhss1 = inner(tau*n,u_exact)*ds(0) - inner(f,phiN*v)*dx
solve(lhss1 == rhss1,suN)
sigmaN,uN = split(suN)
#Stilde(phi,u) = phitilde
lhss2 = inner(theta(norma(grad(phiN)))*grad(phi),grad(psi))*dx - inner(phi*uN,grad(psi))*dx
rhss2 = inner(gamma(phiN)*kv,grad(psi))*dx + inner(g,psi)*dx
phiN1 = Function(WStilde)
solve(lhss2 == rhss2,phiN1,bcs)
phiN.assign(phiN1)
print(i)
import matplotlib.pyplot as plt
plot(sigmaN)
plt.show()
plot(uN)
plt.show()
plot(phiN)
plt.show()
```

What am I doing wrong?

I must notice that I’m not imposing any boundary condition in the first variational problem (although I should impose \int_\Omega \operatorname{tr} \sigma = 0, but I don’t know how to do that)