solving Eikonal equation

On Wed, 01 May 2013 20:16:21 -0000
Pary Parsian <email address hidden> wrote:
> New question #227969 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227969
>
> I am a new Fenics user, and I want to solve eikonal equation, but I
> don't know what is is the problem with my code: from dolfin import *
> import numpy
> import sys
>
> mesh = UnitInterval(51, )
>
> class RefractiveIndex(Expression):
> def eval(self, values,x):
> if x[0]<0.5:
> values[0] = 1.0
> else:
> values[0] = 1.0
>
> def light_exposed_area(x, on_boundary):
> return near(x[0],0.0) and on_boundary
> light_boundary = AutoSubDomain(light_exposed_area)
>
> V = FunctionSpace(mesh,'CG',1)
> u = Function(V)
> du = TrialFunction(V)
> v = TestFunction(V)
> eps = Constant(0.01)
> n = RefractiveIndex()
>
>
> F1 = inner(grad(u), grad(v))*dx - n**2*v*dx

This is only rank 1 form. It does not contain trial function so
lhs(F1) is empty. Use u = TrialFunction(V)

> a1, L1 = lhs(F1), rhs(F1)
> A1, b1 = assemble_system(a1, L1)
>
> for bc in bcs: bc.apply(A1, b1)

You did not define bcs.

> solve(A1, u.vector(), b1)
>
>
> F = inner(grad(u), grad(u))*v*dx - n**2*v*dx +
> eps*inner(grad(u),grad(v))*dx J = derivative(F,u,du)
>
> bc = DirichletBC(V,Constant(0.0),light_boundary)
>
> problem = NonlinearVariationalProblem(-F, u, bcs=[bc,], J=J))

You provide Jacobian of F but tell DOLFIN that it is Jacobian of -F. Do
it consistently.

> solver = NonlinearVariationalSolver(problem)
> solver.solve()
> plot(u)
>

Oh, I made some mistakes! I have boundary condition just in one point (0,1)!
My code is not converging even thought I am using initial guess!
Thanks

from dolfin import *
import numpy
import sys

# Define parameters
nx = 4
ny = 4

# Create mesh and define function space
mesh=Rectangle(-1, -1, 1, 1, nx, ny,"left/right")
#Definig boundary condition
V = FunctionSpace(mesh,'CG',1)
u0 = Constant(0.0)
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14 # tolerance for coordinate comparisons
        return on_boundary and \
               (abs(x[0]) < tol and abs(x[1] - 1) < tol)

u0_boundary = DirichletBoundary()
bc = DirichletBC(V, u0, u0_boundary)

u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
eps = Constant(0.01)

f = Constant(1.0)
# Initialization problem to get good initial guess for nonlinear problem:
F1 = inner(grad(du), grad(v))*dx - f*v*dx
a1, L1 = lhs(F1), rhs(F1)
A1, b1 = assemble_system(a1, L1)
# Apply boundary conditions: DirichletBC = 0 on the boundary
bc.apply(A1, b1)
solve(A1, u.vector(), b1)

F = inner(grad(u), grad(u))*v*dx - f*v*dx + eps*inner(grad(u),grad(v))*dx
J = derivative(F,u,du)

problem = NonlinearVariationalProblem(-F, u, bc, J=J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
plot(u)

Thanks Johan Hake, that solved my question.

Dear all,
I used the implementation described here to solve the eikonal equation and it works just fine.
Does any of you has any (citable) reference to a paper in which this method is discribed (and maybe analyzed)?

The only reference I found is which is kind of close:
Differential equation-based wall distance computation for
DES and RANS by P.G. Tucker, Journal of Computational Physics 190 (2003) 229–248.

Thanks and best regards,
Jan