how to manually build a mesh and use it to solve a pde

Question

I'm trying to build a mesh manually because I am linking Fenics into another project I'm working on where a mesh structure is already defined. It makes more sense to build the mesh in code rather than import two separate files.

Further, I would like to be able to assign the vertex and cell indices so that they correspond with the other mesh data structure that I'm already using.

When I try to do this, I am getting the error

*** Error: Unable to set given (local) rows to identity matrix.
*** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where: This error was encountered inside PETScMatrix.cpp.
*** Process: 0

---

*** DOLFIN version: 2016.1.0
*** Git changeset: unknown

Here is a minimum working example

#include <iostream>
/*---------------------------------------------------------------------------*/
#include <dolfin.h>
#include "pde.h"
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public dolfin::SubDomain
{
  bool inside(const dolfin::Array<double>& x, bool on_boundary) const
  { return on_boundary; }
};

// a function that gives the representation vector at every point on the boundary
class RepresentationVectorFunction : public dolfin::Expression
{
  public:
    RepresentationVectorFunction(std::size_t dim):dolfin::Expression(dim){}

  private:
    void eval(dolfin::Array<double> &values, const dolfin::Array<double> &x) const
    {
      values[0] = 0;
      values[1] = 1;
    }
};


/*---------------------------------------------------------------------------*/
int main()
{
  //Create the mesh and open for editing
  auto mesh = std::make_shared<dolfin::Mesh>();
  dolfin::MeshEditor editor;
  editor.open(*mesh, 2, 2);

  std::vector<double> v = {0,0   //0
                          ,0,1   //1
                          ,0,2   //2
                          ,1,0   //3
                          ,1,1   //4
                          ,1,2   //5
                          ,2,0   //6
                          ,2,1   //7
                          ,2,2}; //8

  editor.init_vertices(v.size());

  int vertex_id = 0;
  for (unsigned int i = 0; i < v.size(); i += 2)
  {
    editor.add_vertex(vertex_id++,v[i],v[i+1]);
  }

  int num_triangles = 8;
  editor.init_cells(num_triangles);

  editor.add_cell(0, 0,1,4);
  editor.add_cell(1, 0,4,3);
  editor.add_cell(2, 1,2,5);
  editor.add_cell(3, 1,5,4);
  editor.add_cell(4, 3,4,7);
  editor.add_cell(5, 3,7,6);
  editor.add_cell(6, 4,5,8);
  editor.add_cell(7, 4,8,7);


  // Close the mesh for editing
  editor.close();

  dolfin::plot(mesh);
  dolfin::interactive();

  // Create function space
  auto V = std::make_shared<pde::FunctionSpace>(mesh);

  auto f = std::make_shared<dolfin::Constant>(0,0);

  // Define boundary condition
  auto u0 = std::make_shared<RepresentationVectorFunction>(2);
  auto boundary = std::make_shared<DirichletBoundary>();
  dolfin::DirichletBC bc(V, u0, boundary);

  // Define pde variational problem
  pde::BilinearForm a(V, V);
  pde::LinearForm L(V);
  L.f = f;

  // Compute pde solution
  auto u = std::make_shared<dolfin::Function>(V);
  dolfin::solve(a == L, *u, bc);

  return 0;
}

The UFL file is

# Elements
V = VectorElement("Lagrange", triangle, 2)

# Trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
f = Coefficient(V)

# Bilinear form
a = inner(grad(u), grad(v))*dx

# Linear form
L = dot(v,f)*dx

Answer 1

Hi, init_vertices should take the number of vertices, in your case v.size()/2.

Comments

Thanks for your answer MiroK. That did in fact fix the code I posted, however in my original code I know for sure that am passing the right number to init_vertices, but I'm still getting the same error.

Here is my original code (which won't be able to run without linking to the rest of the project). It should be fairly straightforward to read. The variable m_mesh is the other mesh data structure that I am using to build the mesh in fenics. Would you mind looking through it to see if you notice anything that would be causing this error?

 //Create the mesh and open for editing
  auto mesh = std::make_shared<dolfin::Mesh>();
  dolfin::MeshEditor editor;
  editor.open(*mesh, 2, 2);

  int num_nodes = m_mesh->getNbNodes();
  editor.init_vertices(num_nodes);

  std::map<gmds::TCellID,std::size_t> node_map;

  int vertex_id = 0;
  gmds::IGMesh::node_iterator node_it = m_mesh->nodes_begin();
  for (; !node_it.isDone(); node_it.next())
  {
    gmds::Node n = node_it.value();
    node_map[n.getID()] = vertex_id;
    editor.add_vertex(vertex_id++,n.X(),n.Y());
  }

  int num_triangles = m_mesh->getNbFaces();
  editor.init_cells(num_triangles);

  int cell_id = 0;
  gmds::IGMesh::face_iterator f_it = m_mesh->faces_begin();
  for(; !f_it.isDone(); f_it.next())
  {
    gmds::Face f = f_it.value();
    std::vector<gmds::Node> nodes = f.get<gmds::Node>();

    std::size_t v0,v1,v2;
    v0 = node_map[nodes[0].getID()];
    v1 = node_map[nodes[1].getID()];
    v2 = node_map[nodes[2].getID()];
    editor.add_cell(cell_id++,v0,v1,v2);
  }


  // Close the mesh for editing
  editor.close();

  dolfin::plot(mesh);
  dolfin::interactive();


  // Create function space
  auto V = std::make_shared<initial::FunctionSpace>(mesh);

  auto f = std::make_shared<dolfin::Constant>(0,0);

  // Define boundary condition
  auto u0 = std::make_shared<RepresentationVectorFunction>(2,*mesh);
  auto boundary = std::make_shared<DirichletBoundary>();
  dolfin::DirichletBC bc(V, u0, boundary);

  // Define initial variational problem
  initial::BilinearForm a(V, V);
  initial::LinearForm L(V);
  L.f = f;

  // Compute initial solution
  auto u = std::make_shared<dolfin::Function>(V);
  dolfin::solve(a == L, *u, bc);

---

I suggest you inspect the mesh when you are done editing. For example, is the number of vertices/cell (in python mesh.num_vertices(), mesh.num_cells()) as you expect?

---

They both are what I expect. Also when I plot the mesh it looks exactly as I expect.

---

Well it looked okay before too, didn't it. Are there any duplicate coordinates? Anyways, consider also passing in keep_diagonal=true to solve (via parameters)

---

I tried keep_diagonal=true in the parameters but the solver didn't recognize it. I think that it might be a PETSc specific parameter.

I also built my own mesh file and when I run it with that the code runs to completion and everything works as expected. The problem must be in the file that I've been trying to use. Anyway, thanks for the help! Your suggestions helped me narrow down the problem.