# solving Eikonal equation

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.

Thanks,

Pary

from dolfin import *

import numpy

import sys

mesh = UnitInterval(51, )

class RefractiveIndex

def eval(self, values,x):

if x[0]<0.5:

else:

def light_exposed_

return near(x[0],0.0) and on_boundary

light_boundary = AutoSubDomain(

V = FunctionSpace(

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

a1, L1 = lhs(F1), rhs(F1)

A1, b1 = assemble_system(a1, L1)

for bc in bcs: bc.apply(A1, b1)

solve(A1, u.vector(), b1)

F = inner(grad(u), grad(u))*v*dx - n**2*v*dx + eps*inner(

J = derivative(F,u,du)

bc = DirichletBC(

problem = NonlinearVariat

solver = NonlinearVariat

solver.solve()

plot(u)

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Johan Hake

- Solved:
- 2013-05-08

- Last query:
- 2013-05-08

- Last reply:
- 2013-05-05

Jan Blechta (blechta) said : | #1 |

On Wed, 01 May 2013 20:16:21 -0000

Pary Parsian <email address hidden> wrote:

> New question #227969 on DOLFIN:

> https:/

>

> 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

> def eval(self, values,x):

> if x[0]<0.5:

> values[0] = 1.0

> else:

> values[0] = 1.0

>

> def light_exposed_

> return near(x[0],0.0) and on_boundary

> light_boundary = AutoSubDomain(

>

> V = FunctionSpace(

> 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(

>

> bc = DirichletBC(

>

> problem = NonlinearVariat

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

it consistently.

> solver = NonlinearVariat

> solver.solve()

> plot(u)

>

Pary Parsian (pary-karimi) said : | #2 |

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(

u0 = Constant(0.0)

class DirichletBounda

def inside(self, x, on_boundary):

tol = 1E-14 # tolerance for coordinate comparisons

return on_boundary and \

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(

J = derivative(F,u,du)

problem = NonlinearVariat

solver = NonlinearVariat

solver.solve()

plot(u)

Johan Hake (johan-hake) said : | #3 |

Try:

from dolfin import *

nx = 10

ny = 10

# Create mesh and define function space

mesh=RectangleM

#Definig boundary condition

V = FunctionSpace(

u0 = Constant(0.0)

class DirichletBounda

def inside(self, x, on_boundary):

return near(x[0], -1.0) and near(x[1], 1.0)

u0_boundary = DirichletBoundary()

bc = DirichletBC(V, u0, u0_boundary, "pointwise")

u = Function(V)

du = TrialFunction(V)

v = TestFunction(V)

eps = Constant(0.1)

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)

solve(a1==L1, u, bc)

plot(u, interactive=True)

F = sqrt(inner(grad(u), grad(u)))*v*dx - f*v*dx +

eps*inner(

solve(F==0, u, bc)

plot(u, interactive=True)

On 05/04/2013 10:16 PM, Pary Parsian wrote:

> Question #227969 on DOLFIN changed:

> https:/

>

> Status: Answered => Open

>

> Pary Parsian is still having a problem:

> 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(

> u0 = Constant(0.0)

> class DirichletBounda

> 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(

> J = derivative(F,u,du)

>

> problem = NonlinearVariat

> solver = NonlinearVariat

> solver.solve()

> plot(u)

>

Pary Parsian (pary-karimi) said : | #4 |

Thanks Johan Hake, that solved my question.

jfpxtal (jan-fp) said : | #5 |

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:

Diﬀerential 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