# Initial Conditions and Function update on time in mixed Element

Version: Dolfinx 0.6.0 Python
Windows Subsystem for Linux (WSL) in Visual Studio Code

Hello guys, I’ve been trying to solve an Eletro-chemo-mechanical coupling problem, and I’ve been faced with some problems.
I’m trying to update these codes to the current version FENICSX, however, I’m having a lot of difficulties.
Due to the fact that the code is extremely extensive, I chose not to write it in full here, just to present the fundamental points

Some questions are:

1. How to determine an initial concentration the of mobile ions c0= 2740 , in time t =0 for each subdomain ( hydrogel and elastomer) , for a mixed elements.

2. How to update scalar and vector solution functions correctly in newmark functions (displacement, speed and acceleration)

In the previous version,
This was done in the following way

``````from dolfin import *
import numpy as np

# Define function space, scalar
U2 = VectorElement("Lagrange", mesh.ufl_cell(), 1)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
D0 = FiniteElement("DG", mesh.ufl_cell(), 0)
D1 = FiniteElement("DG", mesh.ufl_cell(), 1)

# DOFs
TH = MixedElement([U2, P1, P1, P1])
ME = FunctionSpace(mesh, TH) # Total space for all DOFs

W = FunctionSpace(mesh,P1)   # Scalar space for visualization later
W2 = FunctionSpace(mesh,U2)   # Vector space for visualization later
W3 = FunctionSpace(mesh,D0)   # DG space for visualization later
W4 = FunctionSpace(mesh,D1)   # DG space for visualization later

# Define test functions in weak form
dw = TrialFunction(ME)
(u_test, omgPos_test, phi_test, omgNeg_test)  = TestFunctions(ME)    # Test function

# Define actual functions with the required DOFs
w = Function(ME)
(u, omgPos, phi, omgNeg) = split(w)    # current DOFs

# A copy of functions to store values in last step for time-stepping.
w_old = Function(ME)
(u_old, omgPos_old, phi_old, omgNeg_old) = split(w_old)   # old DOFs

v_old = Function(W2)
a_old = Function(W2)

# Initial electro-chemical potentials
init_omgPos = Expression('abs(x[1])>=int2-tol?std::log((cPos0)):std::log((cNum))', int2=int2, tol = tol, cPos0 = cPos0, cNum=DOLFIN_EPS, degree=0)

omgPos_init = interpolate(init_omgPos,ME.sub(1).collapse())
assign(w_old.sub(1),omgPos_init)
assign(w.sub(1),omgPos_init)

# Update formula for acceleration
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dk
beta_ = beta
else:
dt_ = float(dk)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old

# Update formula for velocity
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dk
gamma_ = gamma
else:
dt_ = float(dk)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)

def update_fields(u_proj, u_proj_old, v_old, a_old):
"""Update fields at the end of each time step."""

# Get vectors (references)
u_vec, u0_vec  = u_proj.vector(), u_proj_old.vector()
v0_vec, a0_vec = v_old.vector(), a_old.vector()

# use update functions using vector arguments
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)

# Update (u_old <- u)
v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
#u_old.vector()[:] = u_proj.vector()

def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new

# Update fields for next step
u_proj = project(u, W2)
u_proj_old = project(u_old, W2)
update_fields(u_proj, u_proj_old, v_old, a_old)
w_old.vector()[:] = w.vector()

``````

I made some attempts to solve the problems, however I am not making any progress in terms of updating the vector fields and applying the initial condition, due to the fact that I must determine the initial concentration in different subdominions, for a geometry built in gmsh.
I made two attempts to determine the concentrations, one from the coordinates and the other through physical entities built in gmsh itself.Below is the summarized code.

From the coordinates:

``````from mpi4py import MPI
from dolfinx import fem, io, nls, log, mesh, plot
import numpy as np
import pyvista
from ufl import VectorElement,FiniteElement,MixedElement,TestFunction,TrialFunction,split,grad,tr,Identity,inner,dot
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx.io import gmshio

domain, cell_tags, facet_tags = gmshio.read_from_msh("capacitor_malha.msh", MPI.COMM_SELF,0, gdim=2)

"""Function Space"""
# Define function space, scalar
U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1)  # Displacent
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)  # concentrações + - , electric potential
D0 = ufl.FiniteElement("DG", domain.ufl_cell(), 0)
D1 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)

# DOFs
TH = ufl.MixedElement([U2, P1, P1, P1])
ME = fem.FunctionSpace(domain, TH)  # Total space for all DOFs

"""Extraindo os sub espaços do elemento misto e os mapas contendo os graus de liberdade """
num_subs = ME.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
space_i, map_i = ME.sub(i).collapse()
spaces.append(space_i)
maps.append(map_i)

(u_test, omgPos_test, phi_test, omgNeg_test)  = ufl.TestFunctions(ME)    # Test function

# Define test functions in weak form
dw = ufl.TrialFunction(ME)

# Define actual functions with the required DOFs
w = fem.Function(ME)
(u, omgPos, phi, omgNeg) = ufl.split(w)  # current DOFs

# A copy of functions to store values in last step for time-stepping.
w_old = fem.Function(ME)
(u_old, omgPos_old, phi_old, omgNeg_old) = ufl.split(w_old)  # old DOFs

"""Concentrações iniciais"""
x = ufl.SpatialCoordinate(domain)

DOlfin_Eps= 0.3e-16
scaleY = 900.e0 #Problema aqui
tol = 1e-12
cPos0 = 0.274
int1 = 200.e0 #
int2 = scaleY/2.-int1

class init_Omg:
def __init__(self):
self.scaley= scaleY
self.int2= int2
self.tol= tol
self.c= 0
self.CNum= DOlfin_Eps

def init_Omg_eval(self,x):
if x[1]>= 700e-6 + tol:
#return ufl.ln(self.c) #
return np.log(self.c, dtype=np.int32)
if x[1]<= 200e-6 - tol:
return np.log(self.c, dtype=np.int32)
else:
return np.log(self.CNum, dtype=np.int32)

def init_Omg_(self,x):
return ufl.conditional(ufl.ge(x[1],200e+6 + tol), ufl.ln(self.c), ufl.ln(self.CNum))

my_bc_value = init_Omg()
my_bc_value.c=cPos0
bc_expr = fem.Expression(my_bc_value.init_Omg_eval, spaces[1].element.interpolation_points())
w_old.sub(1).interpolate(bc_expr)

``````

by physical entities:

``````
""" Implementation of concentrations in each region"""
def tag(n_tag):
return cell_tags.find(n_tag)

Q = fem.FunctionSpace(domain, D0)

eletrodo_sup = tag(26)
eletrodo_inf = tag(27)
gel_p = tag(28)

# A ordem de entrada das propriedades na lista deve ser equivalente ao espaço no qual o material ocupa dentro do dominio
material_ = [gel_p, eletrodo_inf, eletrodo_sup]

Init_CPos= fem.Function(spaces[1])
Init_CPos.interpolate(my_bc_value.init_Omg_)

w_old.x.array[gel_p]= Init_CPos.x.array
``````

both attempts don’t work.

I guess you are referring to this?

In the code below, you are mixing a whole lot of `ufl` and `numpy`

Let us consider a simple MWE, on a unit square:

``````from IPython import embed
from mpi4py import MPI
from dolfinx import fem, io, nls, log, mesh, plot
import numpy as np
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx.io import gmshio

domain = mesh.create_rectangle(MPI.COMM_WORLD, [[2, 2], [3, 3]], [10, 10])
"""Function Space"""
# Define function space, scalar
U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1)  # Displacent
# concentrações + - , electric potential
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
D0 = ufl.FiniteElement("DG", domain.ufl_cell(), 0)
D1 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)

# DOFs
TH = ufl.MixedElement([U2, P1, P1, P1])
ME = fem.FunctionSpace(domain, TH)  # Total space for all DOFs
w = fem.Function(ME)

a = 2.5
cpos0 = 1.1
cNum = 2

def init_omega(x):
condition = x[1] > a
condition_inverse = np.invert(condition)
return condition*np.log(cpos0) + condition_inverse * np.log(cNum)
w.sub(1).interpolate(init_omega)
``````

Here I’ve created a simply numpy function with the same functionality as

where you replace `a` with `int2-tol` etc.

As long as your example is kept as partial snippets, it is quite hard to give you any further guidance (and it is aso the reason why it has taken so long to write up a reply). I would strongly suggest you try to simplify your code, and have a look at: PETSc NonlinearProblem running slow - #4 by dokken
as it seems rather relevant for the second question.

Thank you for your answer friend, however, the model you implemented didn’t work for me, however, I got a solution to my problem, which apparently satisfied the behavior I expected, below is the model I used as a reference, In case anyone found it useful.

``````
from mpi4py import MPI
from dolfinx import fem,nls, log
import numpy as np

from ufl import VectorElement,FiniteElement,MixedElement,TestFunction,TrialFunction,split,grad,tr,Identity,inner,dot
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx.io import gmshio
from petsc4py import PETSc

domain, cell_tags, facet_tags = gmshio.read_from_msh("malhas/capacitor_geometria__triangulo_13.msh", MPI.COMM_WORLD,0, gdim=2)

"""Function Space"""
# Define function space, scalar
U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1)  # Displacent
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)  # concentrações + - , electric potential
D0 = ufl.FiniteElement("DG", domain.ufl_cell(), 0)
D1 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)

# DOFs
TH = ufl.MixedElement([U2, P1, P1, P1])
ME = fem.FunctionSpace(domain, TH)  # Total space for all DOFs

W = fem.FunctionSpace(domain,P1)   # Scalar space for visualization later

"""Extraindo os sub espaços do elemento misto e os mapas contendo os graus de liberdade """
num_subs = ME.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
space_i, map_i = ME.sub(i).collapse()
spaces.append(space_i)
maps.append(map_i)

"""Aplicando as propriedades dos materiais"""
"""
1 28 "Voltagem"
1 27 "Aterramento"
1 22 "Engaste_esquerda"
1 23 "Engaste_direita"
2 24 "Hidrogel_superior"
2 25 "Hidrogel_inferior"
2 26 "elastomero"
2 27 "dominio"

"""
"""Extrair as tags da malha"""

def tag(n_tag):
return cell_tags.find(n_tag)

"Determinação das propriedades de um material "

def mat_features(function_descontinuo, material, constanste):
space_Mat = fem.Function(function_descontinuo)
for i in range(len(material)):
space_Mat.x.array[material[i]] = np.full_like(material[i], constanste[i], dtype=ScalarType)
return space_Mat

""" Implementação das propriedades em cada regiao"""
Q = fem.FunctionSpace(domain, D0)

hidrogel_sup = tag(24)
hidrogel_inf = tag(25)
elastomero = tag(26)

# A ordem de entrada das propriedades na lista deve ser equivalente ao espaço no qual o material ocupa dentro do dominio
material_ = [elastomero, hidrogel_inf, hidrogel_sup]
Gshear = mat_features(Q, material_, [0.034e-6, 0.003e-6, 0.003e-6])
Kbulk = mat_features(Q, material_, [2000*0.034e-6,2000*0.003e-6, 2000*0.003e-6])
#intInd = mat_features(Q, material_, [1, 0, 0])
matInd = mat_features(Q, material_, [0.0, 1, 1])
Gshear0 = 100.0e-6  # uso na formução fraca e no cconstrutor
Im_gent = mat_features(Q, material_, [90, 300.0, 300.0])

"""Constantes"""
D = 1.0e-2  # 1.0e0                 # Diffusivity [L2T-1]
RT = 8.3145e-3 * (273.0 + 20.0)  # Gas constant*Temp [ML2T-2#-1]
Farad = 96485.0e-6  # Faraday constant [Q#-1]

""" Initial concentrations"""
cPos0 = 0.274  # Initial concentration [#L-3]
cNeg0 = cPos0  # Initial concentration [#L-3]
cMax = 10000 * 1e-9  # 0.00001*1e-9 # 10000e0*1.e-9

"""Eletrical permittivity """
vareps0 = fem.Constant(domain, 8.85e-12 * 1e-6)  #
vareps_num = mat_features(Q, material_, [1.0 ,1.0e3, 1.0e3])  # Permissividade do gel igual a 10000
vareps_r = mat_features(Q, material_, [6.5, 80., 80.])  # eletrical permittivity for material
vareps = vareps0 * vareps_r * vareps_num

# Mass density
rho = fem.Constant(domain, 1e-9)  # 1e3 kg/m^3 = 1e-9 mg/um^3,

# Rayleigh damping coefficients
eta_m = 0.00000  # Constant(0.000005)
eta_k = 0.00000  # Constant(0.000005)

""" Generalized-alpha method parameters """
alpha = 0.0
gamma = 0.5 + alpha
beta = ((gamma + 0.5) ** 2) / 4

""" Simulation time related params (reminder: in microseconds)"""
ttd = 0.01
# Step in time
t = 0.0  # initialization of time
T_tot = 0e6 * 1.0e6  # 200.0 #0.11        # total simulation time

study = "Electrochemical only"
#study = "Capacitive strain sensing"

if study == "Electrochemical only":
T2_tot = 30.0e6
elif study == "Capacitive strain sensing":
T2_tot = 60.0e6

dt = T2_tot / 20  # incrementos de tempo

phi_norm = RT / Farad  # "Thermal" Volt

"""unknown functions, testing functions and step functions"""
(u_test, omgPos_test, phi_test, omgNeg_test)  = ufl.TestFunctions(ME)    # Test function

# Define test functions in weak form
dw = ufl.TrialFunction(ME)

# Define actual functions with the required DOFs
w = fem.Function(ME)
(u, omgPos, phi, omgNeg) = ufl.split(w)  # current DOFs

# A copy of functions to store values in last step for time-stepping.
w_old = fem.Function(ME)
(u_old, omgPos_old, phi_old, omgNeg_old) = ufl.split(w_old)  # old DOFs

"""acceleration and speed"""
W2 = fem.FunctionSpace(domain, U2)
v_old = fem.Function(W2)
a_old = fem.Function(W2)

"""Implementação das Concentrações iniciais, cations e anions, initial concentrations"""
x = ufl.SpatialCoordinate(domain)
DOlfin_Eps= 0.3e-15

# Initial electro-chemical potentials
Init_CPos= fem.Function(Q)
Init_CPos.x.array[material_[0]]= ufl.ln(DOlfin_Eps)
Init_CPos.x.array[material_[1]]= ufl.ln(cPos0)
Init_CPos.x.array[material_[2]]= ufl.ln(cPos0)

Init_CNeg= fem.Function(Q)
Init_CNeg.x.array[material_[0]]= ufl.ln(DOlfin_Eps)
Init_CNeg.x.array[material_[1]]= ufl.ln(cNeg0)
Init_CNeg.x.array[material_[2]]= ufl.ln(cNeg0)

w_old.sub(1).interpolate(Init_CPos)
w_old.sub(3).interpolate(Init_CNeg)

w.sub(3).interpolate(Init_CNeg)
w.sub(1).interpolate(Init_CPos)

``````

Below is the link to the complete code and the mesh/geometry used