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