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)
```