Integrate values and set values to Neumann condition

dolfinx
1

Hi everyone, I’m new of FEnicsX, and I’ve tried to compute galvanic corrosion using FEniCSX.
To calculate that, following procedure will be conducted.

  1. setting initial potential to the field
  2. calculating current density at the boundaries according to the functions, two boundaries for galvanic corrosion
  3. adjusting total current density at the boundaries to zero
  4. setting adjusted current density to Neumann boundary conditions
  5. calculating PDE

Then I have two quenstions,

  1. In adjusting total current density, I have to obtain values of total current density at the boundary. Can I calculate integral of current density at the boundary? For example, can I calculate integral from value of “u” using function “fCDs” at the boundary “boundaryS”?
  2. After obtaining adjusting current density at the boundary, How do I set values to Neumann boundary condition?

Following is a program I wrote halfway.

Thank you in advance for any help.

Takahiro.

import numpy as np
from mpi4py import MPI
from dolfinx.mesh import create_rectangle

# parameters
lx , ly = 1e-2, 1.0e-3
nx, ny = 200, 20

pN = 100.0
F = 96500.0
R = 8.3144
T = 300.0
DO = 1.9e-9
delta = ly
Einit = -0.9
sigma = 0.2
Ecors = -0.3
Ecorz = -0.9
co = 0.25
ic = 0.35

# create mesh
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0.0, 0.0]), np.array([lx, ly])], [nx, ny])

# Functions
from dolfinx.fem import FunctionSpace, Function
from ufl import TestFunction, TrialFunction

V = FunctionSpace(mesh, ("CG", 1))
u = Function(V)
v = TestFunction(V)

# initial value of u
from dolfinx.fem import Constant, Expression
from petsc4py.PETSc import ScalarType

expr = Expression(Constant(mesh, ScalarType(Einit)), V.element.interpolation_points)
u.interpolate(expr)

# boundaries
def boundaryZ(x):
    return np.logical_and(np.isclose(x[1], 0), x[0] < 0.5 * lx)

def boundaryS(x):
    return np.logical_and(np.isclose(x[1], 0), x[0] >= 0.5 * lx)
    

# functions of current density
def fCDs(x):
	return ic * np.exp((x - Ecors) * 0.38 * F / (R * T))

def fCDz(x):
	return ic * np.exp((x - Ecorz) * 1.2 * F / (R * T))