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:
-
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.
-
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.