How to modify the value of a Function object?

Linear Algebra
1

Hello everyone,

I define a Function object of c0/c1 to store the solution after each finite element step. After I obtain the solution c0/c1, I then use the conditional function (enclosed by # at the end of the code below) to make sure that all values in c0/c1 are between 0~3.75. However, when I print out the minimum and maximum value in c0/c1, it is way off the expected range.

Sorry for the long code. I tried my best to make it as short as possible. It is 2 linear FE processes, diffusion and elastic deformation, which are weakly coupled.

from __future__ import print_function
from fenics import *
from mshr import *
from ufl import nabla_div
import numpy as np

#Define material constants
T = 0.2           # final time
num_steps = 20     # number of time steps
dt = T / num_steps # time step size
D0 = Constant(100.0)         #diffusivity
omega = Constant(406.10315)     #diffusivity_2
nu = Constant(0.26)
alpha = Constant(0.157)   #volume expansion ratio due to concentration increase

# Create mesh and define boundaries
mesh = RectangleMesh(Point(0,0), Point(20,20), 50, 50)
def top_boundary(x):
    return near(x[1], 20)
def bottom_boundary(x):
    return near(x[1], 0)
def left_boundary(x):
    return near(x[0], 0)
#define function space
V  = FunctionSpace(mesh, "CG", 1)
Vu = VectorFunctionSpace(mesh, "CG", 3)
####################################################Diffusion Section
c = TrialFunction(V)
v = TestFunction(V)
c1 = Function(V, name="Concentration") #step n+1
c0 = interpolate(Constant(0.0), V) #step n
# Define boundary condition
c_bc = Constant(3.75)
bc = [DirichletBC(V, c_bc, top_boundary)]
#################################################Deformation Section
u = TrialFunction(Vu)
vu = TestFunction(Vu)
U = Function(Vu, name="Total Displacement")
# Define boundary condition
bcu_left = DirichletBC(Vu.sub(0), 0, left_boundary) #roller BC
bcu_bottom = DirichletBC(Vu.sub(1), 0, bottom_boundary) #roller BC
bcu = [bcu_left, bcu_bottom]

#define functions and variational form
def lame_Si(c):
    E = 100.0-20.0*c
    lambda_ = E*nu/(1.0+nu)/(1.0-2.0*nu)
    mu = E/2.0/(1.0+nu)
    kappa = E/3.0/(1.0-2.0*nu)
    return lambda_, mu, kappa

def sigma(u,c):
    eps = sym(grad(u))
    lambda_, mu, kappa = lame_Si(c)
    return lambda_*tr(eps)*1.5*Identity(2) + 2*mu*eps

def sigma_thermo(c):
    lambda_, mu, kappa = lame_Si(c)
    return alpha*c*(3.0*kappa)*Identity(2)

def grad_mean_stress(u,c):
    s = sigma(u,c)
    ms = tr(s)/3.
    return grad(ms)

#########################################################
#Diffusion
a = v*c*dx + dt*D0*dot(grad(v), grad(c))*dx - dt*omega*dot(grad(v), grad_mean_stress(U,c0))*c*dx
L = v*c0*dx
##########################################################
#Deformation
au = inner(grad(vu), sigma(u,c1))*dx
Lu = inner(grad(vu), sigma_thermo(c1))*dx - inner(grad(vu), sigma(U,c1))*dx
##########################################################

# Time-stepping
t = 0.0;
for n in range(num_steps):
    print("Step:", str(n))

    # Compute solution for diffusion
    solve(a == L, c1, bc)
############################################################################
    c_ = conditional(lt(c1,0.0), 0.0, c1) #make sure c>=0                  #
    c1.assign(project(c_,V))                                               #
    c_ = conditional(gt(c1,3.75), 3.75, c1) #make sure c<=3.75             #
    c1.assign(project(c_,V))                                               #
############################################################################
    # Compute solution for deformation
    solve(au == Lu, U, bcu)
    c0.assign(c1)
    #print results
    print("    min c:", c0.vector().get_local().min())
    print("    max c:", c0.vector().get_local().max())

My question is what should I do if I want to make 0<c<3.75 after I get the FE solution?

Thanks!

2

To me, it seems like your solution (c1) before trying to use the conditional, is highly unstable (spurious behavior). So this is something you have to fix regardless of how you would like to bound your function.
The following alternative works for bounding the variables:

ff = File("c1.pvd")
for n in range(num_steps):
    print("Step:", str(n))

    # Compute solution for diffusion
    solve(a == L, c1, bc)
    ff << c1
    indices = np.argwhere(c1.vector()[:] > 3.75)
    c1.vector().vec().setValuesLocal(np.array(indices,dtype=np.int32), np.full(len(indices), 3.75))
    indices = np.argwhere(c1.vector()[:] < 0)
    c1.vector().vec().setValuesLocal(np.array(indices,dtype=np.int32), np.full(len(indices),0))
    c1.vector().apply("insert")

    solve(au == Lu, U, bcu)
    c0.assign(c1)
    #print results
    print("    min c:", c1.vector().get_local().min())
    print("    max c:", c1.vector().get_local().max())
    ff << c1

And as you can observe by inspecting c1.pvd, the output is highly spurious

3

Thank you. I will check my equations and code.