Point source for nonlinear problem
I then defined a delta function using the limit definition. Obviously, this function is not well approximated by CG P1 elements. Would this work for DG?
My next idea is to modify the linearized system to add the source term, but I don't know how to do that.
Has anyone ever implemented a point source for a nonlinear problem? Thank you.
for use of PointSource(), I have the manual Newton's method solved with a point source:
from fenics import * mesh = UnitSquareMesh(50,50) Q = FunctionSpace(mesh, 'CG', 1) u = Function(Q) du = TrialFunction(Q) v = TestFunction(Q) p = Point(0.75,0.75) P = PointSource(Q, p, 100.0) f = Constant(100.0) def boundary(x, on_boundary): return on_boundary bc = DirichletBC(Q, 50.0, boundary) bcs = [bc] R = + inner(grad(u), grad(v)) * dx \ - f * v * dx J = derivative(R, u, du) # apply the boundary condition to the vector once, # the search direction is set to zero on the boundary # with homogenize() below, and therefore 'u' will not # change along the Dirichlet boundaries : for bc in bcs: bc.apply(u.vector()) atol = 1e-7 rtol = 1e-10 relaxation_param = 1.0 max_iter = 25 method = 'mumps' preconditioner = 'default' converged = False lmbda = relaxation_param # relaxation parameter nIter = 0 # number of iterations # need to homogenize the boundary, as the residual is always zero over # essential boundaries : bcs_u =  for bc in bcs: bc = DirichletBC(bc) bc.homogenize() bcs_u.append(bc) # the direction of decent : d = Function(u.function_space()) while not converged and nIter < max_iter: # assemble system : A, b = assemble_system(J, -R, bcs_u) # apply the point source : P.apply(b) # determine step direction : solve(A, d.vector(), b, method, preconditioner) # calculate residual : residual = b.norm('l2') # set initial residual : if nIter == 0: residual_0 = residual # the relative residual : rel_res = residual / (residual_0 + DOLFIN_EPS) # check for convergence : converged = residual < atol or rel_res < rtol # move U down the gradient : u.vector()[:] += lmbda*d.vector() # increment counter : nIter += 1 # print info to screen : if MPI.rank(mpi_comm_world()) == 0: string = "Newton iteration %d: r (abs) = %.3e (tol = %.3e) " \ +"r (rel) = %.3e (tol = %.3e)" print string % (nIter, residual, atol, rel_res, rtol) File('u.pvd') << u
edit: correctly applied the boundary condition to the initial guess to fix issue described by below comment.
tol = 3E-5 class Wells(Expression): def eval(self,value,x): Q = 1E-2 x0 = -25 y0 = 0 x1 = 25 y1 = 0 r0 = sqrt((x - x0)**2+(x - y0)**2) r1 = sqrt((x - x1) ** 2 + (x - y1) ** 2) if r0 < tol: value = Q elif r1 < tol: value = -Q else: value = 0
The tolerance gets a little tricky because, most likely, there is not a node in the position you want, a little bit of debugging is needed to check which one is the closest node and then adjust the tolerance.
To use it in your formulation you just do this, where q is your TestFunction
f=Wells(degree=3) F = (w**3/(12*mu))*inner(grad(q),grad(p))*dx-f*q*dx
Hope this helps!!