Boundary condition

variational formulation
1

hello,

im trying to work on a code which solves two things. the code is below for reference. the first part is a variational problem to solve for (u) and the second part of the code is basically from the tutorial. however, im trying to change it to suit the problem im solving. im trying to make phie (source term) = u (from my first part of the code) and i believe that the boundary condition (g_D) has to change because the source term = −(sum of the second derivation of x and y). i tried placing (u) in the expression for the boundary condition and it cant work, i was wondering if theres a way for me to extract the values of (u) and then place it in the expression for the boundary condition. or if theres another way in solving it

currently for my boundary condition, i just tried ‘1’ and the code runs fine, i just know that its not mathematically correct because my g_D has no (u) in it.

ps: i believe its due to the fact the u is a varying function and thats why i cant just put u in the g_D?

heres the code:


from fenics import *
import math
import os
import numpy as np



end_time = 160.0  
dt = 0.1          
num_steps = 1600
eps = 0.1




mesh = RectangleMesh(Point(0,0), Point(30,30), 10, 10, 'crossed')
dimension = 2

#Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)

V = FunctionSpace(mesh, P1)

#Define test functions
v_1 = TestFunction(V)

#Define functions for velocity and concentrations
u = Function(V)
u_n = Function(V)

output_frequency = 10

vtkfile_phi1  = File('newtrial3/first.pvd')
vtkfile_phi2  = File('newtrial3/second.pvd')

**FIRST PART**
#Define source terms
i_s1 = 0  # initial condition
stim = Expression(
    'x[0] <= 1 ? i_s1 : 0.0', degree = 0, i_s1 = i_s1)

#Define expressions used in variational forms
k = Constant(dt)
eps = Constant(eps)
c = Constant(8.0)
a1 = Constant(0.01)


#Define variational problem
F = ((u - u_n)/k)*v_1*dx \
  + dot(eps*grad(u), grad(v_1))*dx \
  + c*u *(u - a1)*(u- Constant(1.0))*v_1*dx \
  - stim*v_1*dx 
  

 **SECOND PART**
g_D = Expression('1 ', degree=2)
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, g_D, boundary) 
g= TrialFunction(V)
g1= TestFunction(V)
phie= (u)
a= (inner(grad(g), grad(g1))*dx) 
L= phie*g1*dx 


#Compute solution
g = Function(V)


**main loop**

t=0
for n_step in range(num_steps):
    t += dt 

    # The s1 stimulus
    if t <= 0.5:
        stim.i_s1 = 1.0
    else:
        stim.i_s1 = 0.0   
    
    
    # Solve the problems

    solve(F== 0, u)
    
    solve(a == L, g, bc)


    u_n.assign(u)

    vtkfile_phi1 << (u, t)
    vtkfile_phi2 << (g, t)
2

I do not know what your code is supposed to do, but you can put u in the declaration of the DirichletBC, no need to define an Expression:

bc = DirichletBC(V, u, 'on_boundary')