Vertices order

Question

Hi,

I cannot figure out what is the vertice order.

mesh = UnitSquareMesh(4, 4)
mesh.translate(Point(-.5,-.5))

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

c = TrialFunction(W)
u = TrialFunction(V)

nv  = mesh.num_vertices()
arr = numpy.zeros((nv*2,))
crr = numpy.zeros((nv*1,))

coor = mesh.coordinates()
for i in range(nv) :
    x = coor[i][0]
    y = coor[i][1]
    r = sqrt( x**2 + y**2 )     
    arr[i]     = x 
    arr[i*2]  = y 
    crr[i]     = r

u.vector().set_local(arr)
c.vector().set_local(crr)
plot(c)

vector c and vector u aren't what I'm trying to set up.
coor[i] doesn't correspond to what is in c[i].

I'm trying to avoid the use of Expression, because I'll have very complicated expression.

Thanks for any help.

Comments

The creation of a FunctionSpace will reorder the dofs, so that they don't necessarily correspond to the vertex-ordering of the mesh.

How would doing it like this be better than a complicated Expression? It will not save you any evaluations. If your expression is costly to evaluate, consider writing a C++-snippet for it (ref the documentation).

---

Thank you, I didn't know about C++-snippet and I didn't succeed to use them.
But in help(Expression), we can find 3 ways of doing it. I choose python overloading of class Expression.

An exemple code below :

from dolfin import * 
from math import *

def sign(x):
    if x > 0 :
        return 1
    elif  x < 0 :
        return -1
    return 0

# Overload Expression
class MyExpression0(Expression):
    def eval(self, value, x):
        X = x[0] - 0.5
        Y = x[1] - 0.5
        v = 0.1
        r = sqrt( X**2 + Y**2)
        if X != 0 :
            theta = atan(Y/X) 
        else :
            theta = pi / 2

        value[1] = -r*sign(X)*cos(theta)
        value[0] =  r*sign(X)*sin(theta)
    def value_shape(self):
        return (2,)

# define mesh
mesh = UnitSquare(10, 10)

# define function space
V = VectorFunctionSpace(mesh, "Lagrange", 2)


# Evaluate and plot  MyExpression0
f0 = MyExpression0()
U = interpolate(f0, V)
plot(U)

# Hold plot
interactive()

This exemple code defines a ortho-radial vector field.

---

I still don't understand why you would want to set the vector explicitly. As long as you're using python, I cannot see any significant speed-ups.

---

Well, it is not about speed-up. It's about ease of use.

I'd like to perform it via C++ snippet as well because speed-up could become an issue.