Regional initial condition - a trial w/ if statement

[Kaa]

Failed Matter: Setting different initial condition for a region.

For a rectangular domain on (0, 0) (10, 10), I want to set different initial conditions for a variable’s (u1) initial condition. For a rectangle of (0, 0) (5, 5), IC should be 1 and 0.5 for the rest of the region.

Failed Logic:

• IC is set as an expression in C++ syntax.

• Add a spatial if expression for an IC.

MWE:

from fenics import *

# Time features
T = 50 #-> 300
num_steps = 100 #-> 600
dt = num_steps / T

mesh = RectangleMesh(Point(0.0, 0.0), Point(10.0, 10.0), 100, 100)
P1 = FiniteElement('CG', triangle, 2)
W = FunctionSpace(mesh, MixedElement([P1,P1]))


# Boundary
def boundary(x):
  return near(x[0], 0.0, DOLFIN_EPS) or near(x[1], 0.0, DOLFIN_EPS)\
    or near(x[0], 10.0, DOLFIN_EPS) or near(x[1], 10.0, DOLFIN_EPS)

# Conditions
u10 = Constant(0.5)
u20 = Constant(2)
bc1 = DirichletBC(W.sub(0), u10, boundary)
bc2 = DirichletBC(W.sub(1), u20, boundary)
bcs = [bc1, bc2]

# Initial Condition 
u_0 = Expression(('((x[0] <= (5.0 + DOLFIN_EPS)) and x[1] <= (5.0 + DOLFIN_EPS))) ? 1 : 0.5', '2'), degree=2)

# Test - trial functions
(v1, v2) = TestFunctions(W)
w = Function(W)
u1, u2 = split(w)
u_n = project(u_0, W)
u_n1, u_n2 = split(u_n) 

# Eq.s
F1 = u1*v1*dx - u_n1*v1*dx + 0.5*(dt*dot(grad(u1), grad(v1))*dx + dt*dot(grad(u_n1), grad(v1))*dx)
F2 = u2*v2*dx - u_n2*v2*dx + 0.5*(dt*dot(grad(u2), grad(v2))*dx + dt*dot(grad(u_n2), grad(v2))*dx)
F = F1 + F2 

# Save
file_u1 = File("u1.pvd")
file_u2 = File("u2.pvd")


#Time-stepping
t = 0 

for n in range(num_steps): 

    #Update 

    t += dt 
    # Solve
    solve(F == 0, w, bcs)
    # Save sol
    _u_1, _u_2 = w.split() 
    _u_1.rename('u_1', 'x')
    _u_2.rename('u_2', 'y')
    file_u1 << (_u_1, t)
    file_u2 << (_u_2, t)
    u_n.assign(w)

Does anyone know how this can be done as simple as possible? Thank you in advance.

[dokken]

How about:

u_0 = Expression(("x[0]<=(5.0 + DOLFIN_EPS) && x[1] <= (5.0 + DOLFIN_EPS) ? 1:0.5", "2"), degree=2)

as and is not C++ syntax

[Kaa]

What a reasonable suggestion! Thanks.

[VINICIUS_DANIEL_PERE]

DOLFINX VERSION .
Hello, I need these same condition models, y>= int2-tol, I have a mesh, divided into hydrogels and elastomers, where I need to define the concentrations of cations and anions, which are different in each region, could anyone help me?

"""Concentrações iniciais"""

x = ufl.SpatialCoordinate(domain)
DOlfin_Eps= 0.3e-17
tol = 1e-12
scaleX = 2.0e4   # 2 cm
scaleY = 900.e0  # 0.9 mm
# N number of elements in y-direction    
int1 = 200.e0 # 
int2 = scaleY/2.-int1
int3 = int1/2.0 + int2

class InitOmg():
    def __init__(self):
        self.t= 0.0
        self.Tramp =T2_tot
        self.scalex = scaleX
        self.cPos0= cPos0
        self.cNeg0= cNeg0
        self.cNum= DOlfin_Eps
    
    def eval(self,x):
        if abs(x[1])>= int2-tol:
            return ufl.ln((self.cNeg0))
        else:
            return ufl.ln((self.cNum))
       
Init_POS= InitOmg()
POs= fem.Function(spaces[1])
POs.interpolate(Init_POS.eval)
w_old.sub(1).interpolate(POs)

I have already tried different ways, including using this format, however, the concentrations found do not satisfy the desired solution.

domain, cell_tags, facet_tags = gmshio.read_from_msh("malhas/capacitor_geometria_pronta_g.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)

"""Concentrações iniciais"""
x = ufl.SpatialCoordinate(domain)
DOlfin_Eps= 0.3e-17

# 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]
# 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)
w_old.sub(1).interpolate(Init_CPos)

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(3).interpolate(Init_CNeg)

here is the model that implemented it correctly, however it is in the version where the Expression argument is obsolete

DOLFIN VERSION
# 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)

below we have the expected behavior.

and subsequently the incorrect one, it is worth highlighting that the division is hydrogel, elastomer, hydrogel

I can open another call and make the complete code available if necessary for better interpretation, thank you in advance for any effort to help me

[dokken]

Please open a new post with a reproducible minimal example.
The code you have supplied is not executable for anyone, making it hard to give you guidance.