initializing mesh with user defined coordinates and connectivity

Question

If I specify

coord = numpy.array([ [ 0, 0], [1, 0], [0, 1], [1,1] ])
connectivity = numpy.array([ [0, 1, 3], [1, 2, 3]  ])
mesh = Mesh()
mesh.coordinates()[:] = coord
mesh.cells()[:] = connectivity
plot(mesh)

it gives an error as below:

Traceback (most recent call last):
File "2d_mesh_expts.py", line 141, in <module>
mesh.coordinates()[:] = coord
ValueError: operands could not be broadcast together with shapes (0,0)  (4, 2) 

Is there any way to do this ? Once I initialize the mesh, I plan to refine it to desired level. ( Of course this is a trivial MWE)

Answer 1

For this particular example you can just create a RectangleMesh in dolfin:

 mesh = RectangleMesh(0.0, 0.0, 1.0, 1.0, 1, 1, 'right')

(here 'right' specifies the direction of the diagonal).

The reason why your example doesn't work is because when you say mesh=Mesh() this creates an empty mesh, which is why its coordinates array has shape (0,0) so you can't just assign a bigger array of coordinates. However, you can add vertices and cells by using dolfin's MeshEditor class:

from dolfin import *
import numpy as np
editor = MeshEditor()
mesh = Mesh()
editor.open(mesh, 2, 2)  # top. and geom. dimension are both 2
editor.init_vertices(4)  # number of vertices
editor.init_cells(2)     # number of cells
editor.add_vertex(0, np.array([0.0, 0.0]))
editor.add_vertex(1, np.array([1.0, 0.0]))
editor.add_vertex(2, np.array([0.0, 1.0]))
editor.add_vertex(3, np.array([1.0, 1.0]))
editor.add_cell(0, np.array([0, 1, 3], dtype=np.uintp))
editor.add_cell(1, np.array([0, 2, 3], dtype=np.uintp))
editor.close()
plot(mesh)

Note that the vertex arrays must have dtype double and the cell arrays must have dtype uintp.

I think there is also an example demonstrating the MeshEditor class in the manual or the demos.

Comments

Thanks! That seems to work.I was aware of the RectangleMesh , but wanted to know ifthis could be done.
I have 2 followup questions though:

• After editor(close) , when I try refine(mesh), there is no refinement. What is the problem there ?

• Subjective question: How do I find out information , like the
  MeshEditor for instance ? I find am completely reliant on
  information given in this forum and the tutorial/demo programs. I
  dont remember seeing this in the tutorial or the demos for instance.

---

When you say refine(mesh) it doesn't refine it in place but returns the refined mesh, so you need to say:

mesh_fine = refine(mesh)

Regarding the second question, the MeshEditor is mentioned in the FEniCS manual (which can be downloaded in pdf form here: http://fenicsproject.org/documentation/). It's described in section 3.3.2 "Meshes" (p. 100). I recommend browsing at least the table of contents of the FEniCS manual and the FEniCS book (which can be obtained from here: http://fenicsproject.org/book/index.html#book by following the link "Download from Launchpad") and start by reading sections that seem to be interesting or relevant to your application.