Mesh Coordinates Index doesn't Match Function Index in FEniCS 1.4.0

Question

This question is related to my previous posted question
[http://fenicsproject.org/qa/5359/zero-dirichlet-boundary-condition-doesnt-work][1]. After some tests, I found the problem is from the mismatch between mesh coordinates index doesn't match function index in FEniCS 1.4.0

Here is the simple test code:

from dolfin import *
from math import *
import numpy as np
import scipy as sc
import scipy.linalg
import datetime
import sys

bl=1.4

mesh = BoxMesh(-50,-50,-50,50,50,50,2,2,2)
origin1 = Point(-bl,.0,.0)
origin2 = Point(bl,.0,.0)

V = FunctionSpace(mesh, "Lagrange", 1)

#relationship between vertex and element                        
e_v=mesh.cells()
noe=e_v.shape[0]

#coordinates of vertex                                          
coor = mesh.coordinates()
dof=coor.shape[0]

v_ext=Expression("1.0/(pow(pow(x[0]-bl,2)+pow(x[1]-0.0,2)+pow(x[2]-0.0,2),0.5))+1.0/(pow(pow(x[0]+bl,2)+pow(x[1]-0.0,2)+pow(x[2]-0.0,2),0.5))",bl=bl)

for i in range(dof):
    print coor[i],(interpolate(v_ext,V)).vector().array()[i]

Here are the results in FEniCS 1.4.0

[-50. -50. -50.] 0.0230940101543
[ 0. -50. -50.] 0.0282870413401
[ 50. -50. -50.] 0.0399843292138
[-50.   0. -50.] 0.0282787291596
[  0.   0. -50.] 0.0282870413401
[ 50.   0. -50.] 0.0230940101543
[-50.  50. -50.] 1.42857142857
[  0.  50. -50.] 0.0399843292138
[ 50.  50. -50.] 0.0282870413401
[-50. -50.   0.] 0.0230940101543
[  0. -50.   0.] 0.0400313846055
[ 50. -50.   0.] 0.0282787291596
[-50.   0.   0.] 0.0230940101543
[ 0.  0.  0.] 0.0230940101543
[ 50.   0.   0.] 0.0282787291596
[-50.  50.   0.] 0.0399843292138
[  0.  50.   0.] 0.0282870413401
[ 50.  50.   0.] 0.0399843292138
[-50. -50.  50.] 0.0282870413401
[  0. -50.  50.] 0.0400313846055
[ 50. -50.  50.] 0.0282870413401
[-50.   0.  50.] 0.0282870413401
[  0.   0.  50.] 0.0230940101543
[ 50.   0.  50.] 0.0282870413401
[-50.  50.  50.] 0.0230940101543
[  0.  50.  50.] 0.0282787291596
[ 50.  50.  50.] 0.0230940101543

Here are the results in FEniCS 1.1.0

[-50. -50. -50.] 0.0230940101543
[  0. -50. -50.] 0.0282787291596
[ 50. -50. -50.] 0.0230940101543
[-50.   0. -50.] 0.0282870413401
[  0.   0. -50.] 0.0399843292138
[ 50.   0. -50.] 0.0282870413401
[-50.  50. -50.] 0.0230940101543
[  0.  50. -50.] 0.0282787291596
[ 50.  50. -50.] 0.0230940101543
[-50. -50.   0.] 0.0282870413401
[  0. -50.   0.] 0.0399843292138
[ 50. -50.   0.] 0.0282870413401
[-50.   0.   0.] 0.0400313846055
[ 0.  0.  0.] 1.42857142857
[ 50.   0.   0.] 0.0400313846055
[-50.  50.   0.] 0.0282870413401
[  0.  50.   0.] 0.0399843292138
[ 50.  50.   0.] 0.0282870413401
[-50. -50.  50.] 0.0230940101543
[  0. -50.  50.] 0.0282787291596
[ 50. -50.  50.] 0.0230940101543
[-50.   0.  50.] 0.0282870413401
[  0.   0.  50.] 0.0399843292138
[ 50.   0.  50.] 0.0282870413401
[-50.  50.  50.] 0.0230940101543
[  0.  50.  50.] 0.0282787291596
[ 50.  50.  50.] 0.0230940101543

FEniCS 1.1.0 gives the right results, while FEniCS 1.4.0 gives different results. How can we get the relationship between vertex coordinates and function values at a certain vertex in FEniCS 1.4.0?

Answer 1

Hi, you need to work with DofMap object. You can then obtain coordinates of dofs with different tabulate methods, see here. Or you can use dof_to_vertex_map as done here. Note that the function is now free so the linked answer needs a slight tweak.

Comments

Thanks for the explanation~~

Answer 2

Hi,

Here is an example of how to access the correct coordinate. The script below assigns the value x+y+z to each vertex at (x,y,z) in a box mesh.

from dolfin import *
import numpy as np

mesh = BoxMesh(-1,-1,-1,1,1,1,20,20,20)
V = FunctionSpace(mesh, "Lagrange", 1)
func = Function(V)

n = V.dim()
d = mesh.geometry().dim()

coor = V.dofmap().tabulate_all_coordinates(mesh)
coor.resize((n,d))

vertex_values = np.zeros(mesh.num_vertices())
for vertex in vertices(mesh):
    x = vertex.x(0)
    y = vertex.x(1)
    z = vertex.x(2)

    vertex_values[vertex.index()] = x+y+z

func.vector()[:] = vertex_values[dof_to_vertex_map(V)]

plot(func,interactive=True)

Comments

Thank you very much~~