Duplicated coordinates with P2 elements

Question

I don't understand why the tabulated coordinates are duplicated when using P2 elements.

# Mesh details 
uplim = 1
size = 50
# Create Mesh    
domain = Rectangle(Point(0., 0.), Point(uplim, uplim))
mesh = generate_mesh(domain,size)
# Define Function Space
V = VectorFunctionSpace(mesh, 'CG', pdim, dim=2)

# Cartesian coordinates:
dim = V.dim()
print dim
N = mesh.geometry().dim()
coor = V.tabulate_dof_coordinates().reshape(dim,N)

print coor[0:10,:]
print coor.size

This shows that each 2D coordinate is doubled. When using p2 elements I would expect a higher number of points (not necessarily doubled with respect to p1 elements), some of which are nodes and some others are edges middle points.

Moreover, what can I use in python to obtain the size of an object as an ndarray? size gives the number of all elements, while I'm looking for an output of type (dim,N).

Thanks

Answer 1

Hi,

to answer your two questions:
1) This is exactly the expected behavior, since you are initializing V as a 2D vector function space (since dim=2). This means that at each dof position, you actually have two components/dofs (or even more, for higher dimensional VectorFunctionSpaces and TensorFunctionSpaces ). So when listing the coordinates of the dofs as in your example, you will indeed observe that different dofs share the same coordinates.

2) To get the shape of an ndarray, simply use:

coor.shape

returning a tuple of the ndarray shape.

Comments

Hi,

Thank you, my bad: I confused pdim that indicates the element type (linear/quadratic), with dim that gives the dimension of the system.

The output variable is a 2D vector u=[u_x,u_y]^T, and u_x and u_y share the same coordinate; in which case this would not be the case? I don't understand the interest of having double dofs/components.

My goal is to set a condition over the coordinates and then save the indices of the dofs that correspond to this request. To do so I need to skip all the even entries of coord.

I've done:

xs1 = 0.495
xs2 = 0.505
indices = np.where(np.logical_and(coor_x > xs1, coor_x < xs2))[0]
xs = coor[indices,:] 

vals = u1_proj.vector()[indices] #

source_vec = np.zeros((len(xs)/2,4))
count = 0
for i in range(0,len(xs)-1,2):
    source_vec[count,:] = [ xs[i,0], xs[i,1], vals[i], vals[i+1] ]
    count = count+1

where u1_proj is defined as u1_proj.vector()[:] = project(u1, V).vector() and u1 is a function given by the sum of the solution and other not relevant stuff for this question.

Is there a better way to do so???