How to perform non-smooth operations on variables in finite element function space

variational formulation
1

When I tried to explore the coupling between the cahnhilliard equation and the ns equation, I found that I did not know how to perform nonsmooth operations on variables in the finite element space, such as when the relationship between density and concentration is a piecewise linear function.
The following is my code:

import random
from numpy.ma import tanh
import numpy as np
from fenics import *
from matplotlib import pyplot as plt


# piecewise function 
def den(c_):
    rho = np.piecewise(
        c_, [0 <= c_ < c_ox, c_ox <= c_ < c_met, c_met <= c_ <= 1],
        (lambda c_: (rho_eq_ox - rho_int_ox) / c_ox * c_ + rho_int_ox,
         lambda c_: (rho_eq_met - rho_eq_ox) / (c_met - c_ox) * (c_ - c_ox) + rho_eq_ox,
         lambda c_: (rho_eq_met - rho_int_met) / (c_met - 1) * (c_ - c_met) + rho_eq_met
         )
    )
    return rho


class InitialConditions(UserExpression):
    def eval(self, values, x):
        x_ = x[0]
        y_ = x[1]
        random.seed(x_)
        rx = 2 * random.random() - 1
        values[0] = 0.5 * (1 + tanh(2 / wint * (y_ - L - 0.5e-4 * L * rx)))
        values[1] = 0.0

    def value_shape(self):
        return 2,


# Sub domain for Periodic boundary condition
class PeriodicBoundaryX(SubDomain):
    def inside(self, x, on_boundary):
        return DOLFIN_EPS > x[0] > (- DOLFIN_EPS) and on_boundary

    def map(self, x, y):
        y[0] = x[0] - L
        y[1] = x[1]


# 
L = 0.176
wint = L / 16
rho_int_ox = 8.24783e3
rho_int_met = 6.96642e3
rho_eq_ox = 7.97994e3
rho_eq_met = 9.01478e3
c_ox = 1.87512e-3
c_met = 6.16085e-1
nx = 192
ny = 384


mesh = RectangleMesh(Point(0, 0), Point(L, 2 * L), nx, ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1 * P1)
du = TrialFunction(ME)
q, v = TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)

u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)

# 
u_0 = u.split(deepcopy=True)[0]
rho_0 = den(u_0.vector().get_local()) #How will I modify the line code?

The following error message is:

Traceback (most recent call last):
  File "/home/wangyao/PycharmProjects/hello_word/main.py", line 88, in <module>
    rho_0 = den(u_0.vector().get_local())
  File "/home/wangyao/PycharmProjects/hello_word/main.py", line 11, in den
    c_, [0 <= c_ < c_ox, c_ox <= c_ < c_met, c_met <= c_ <= 1],
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I hope to calculate the density distribution based on the c distribution.
Here is the relationship between density and concentration：

image
image1463×701 41.9 KB

Thank you for your help!

2

Hi, you can use conditional expressions in UFL https://fenics.readthedocs.io/projects/ufl/en/latest/manual/form_language.html#conditional-operators

1 Like
3

I tried ufl conditional expression, I’m not sure how to write it, here’s my code：

import random
from numpy.ma import tanh
from matplotlib import pyplot as plt
from ufl import *
from fenics import *  



# def den functions

def left(c_):
    return (rho_eq_ox - rho_int_ox) / c_ox * c_ + rho_int_ox


def mid(c_):
    return (rho_eq_met - rho_eq_ox) / (c_met - c_ox) * (c_ - c_ox) + rho_eq_ox


def right(c_):
    return (rho_eq_met - rho_int_met) / (c_met - 1) * (c_ - c_met) + rho_eq_met


class InitialConditions(UserExpression):
    def eval(self, values, x):
        x_ = x[0]
        y_ = x[1]
        random.seed(x_)
        rx = 2 * random.random() - 1
        values[0] = 0.5 * (1 + tanh(2 / wint * (y_ - L - 0.5e-4 * L * rx)))
        values[1] = 0.0

    def value_shape(self):
        return 2,


# Sub domain for Periodic boundary condition
class PeriodicBoundaryX(SubDomain):
    def inside(self, x, on_boundary):
        return DOLFIN_EPS > x[0] > (- DOLFIN_EPS) and on_boundary

    def map(self, x, y):
        y[0] = x[0] - L
        y[1] = x[1]



L = 0.176
wint = L / 16
rho_int_ox = 8.24783e3
rho_int_met = 6.96642e3
rho_eq_ox = 7.97994e3
rho_eq_met = 9.01478e3
c_ox = 1.87512e-3
c_met = 6.16085e-1
nx = 192
ny = 384
Vm = 1e-5
R = 8.31446
T = 3000
w_ref = 1e-9
DU = 1.55E-12
DO = 2.99e-9
DZ = DU
DF = DU
M_ref = Vm / (R * T) * (DU + DZ + DO + DF)
M = wint / w_ref * M_ref

mesh = RectangleMesh(Point(0, 0), Point(L, 2 * L), nx, ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1 * P1)
du = TrialFunction(ME)
q, v = TestFunctions(ME)
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)

u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)


u_0 = u.split(deepcopy=True)[0]
plot(u_0, mode='color', vmin=min(u_0.vector().get_local()), vmax=max(u_0.vector().get_local()))
plt.show()

rho_0 = conditional([0 <= u_0 < c_ox, c_ox <= u_0 < c_met, c_met <= u_0 <= 1],
                    [left(u_0), mid(u_0), right(u_0)], [left(u_0), mid(u_0), right(u_0)])

This is an error message：

Traceback (most recent call last):
  File "/home/wangyao/PycharmProjects/hello_word/main.py", line 102, in <module>
    rho_0 = conditional([0 <= u_0 < c_ox, c_ox <= u_0 < c_met, c_met <= u_0 <= 1],
  File "/usr/lib/python3/dist-packages/ufl/conditional.py", line 34, in __bool__
    error("UFL conditions cannot be evaluated as bool in a Python context.")
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.

Can you try to modify it on my basis?
Thank you for your help!

4

see here https://fenicsproject.org/qa/8175/conditional-expression-ufl-condition-as-boolean/

5

I noticed that this example is a typical 2-segment piecewise function. How should I modify my 3-segment piecewise function or more? Thank you for your help!

6

If there are not too many segments, it’s reasonable to nest conditionals, e.g.,

f = conditional(lt(x,x1),f1,
                         conditional(lt(x,x2),f2,
                                              f3))

to implement

f = \left\{\begin{array}{lcr}f_1&,&x<x_1\\ f_2&,&x_1\leq x<x_2\\ f_3&,&x\geq x_2\end{array}\right.
3 Likes
7

Thank you bleyerj and kamensky ,and the problem has been successfully solved!
image