Classifying an existing mesh to mark boundaries uniquely

[a-jp]

Hi,

I'm calling the following code in 3D on an existing mesh. I generated the code mainly from trial and error using the icesheet.cc and classifyThenAdapt.cc examples. The following "works" in that it doesn't fail when I call

  apf::deriveMdsModel(mesh);
  apf::verify(mesh, true);

afterwards. Beyond this I don't know if it's right, and I have a few questions. Can someone take a look and see if what I'm doing is correct?

  const auto mdim = mesh->getDimension();
  const int INTERIOR_CELLS = 0;
  const int INTERIOR_FACES = 1;
  const int BOUNDARY_FACES = 2;
  const int INTERIOR_EDGES = 3;
  const int INTERIOR_VERTS = 4;
  {
    apf::ModelEntity *interior = mesh->findModelEntity(mdim, INTERIOR_CELLS);
    apf::MeshIterator *it = mesh->begin(mdim);
    apf::MeshEntity *ent;
    while ((ent = mesh->iterate(it)))
      mesh->setModelEntity(ent, interior);
    mesh->end(it);
  }
  {
    apf::ModelEntity *interior = mesh->findModelEntity(0, INTERIOR_VERTS);
    apf::MeshIterator *it = mesh->begin(0);
    apf::MeshEntity *ent;
    while ((ent = mesh->iterate(it)))
      mesh->setModelEntity(ent, interior);
    mesh->end(it);
  }
  {
    apf::MeshIterator *it = mesh->begin(mdim-1);
    apf::MeshEntity *ent;
    while ((ent = mesh->iterate(it)))
    {
      const auto numCellNeighbours = mesh->countUpward(ent);
      if (numCellNeighbours == 1)
      {
        mesh->setModelEntity(ent, mesh->findModelEntity(mdim-1, BOUNDARY_FACES));
      }
      else
      {
        mesh->setModelEntity(ent, mesh->findModelEntity(mdim, INTERIOR_FACES));
      }      
    }
    mesh->end(it);
  }
  if(mdim == 3)
  {
    apf::ModelEntity *interior = mesh->findModelEntity(1, INTERIOR_EDGES); // Why one in the volume?? 
    apf::MeshIterator *it = mesh->begin(1);
    apf::MeshEntity *ent;
    while ((ent = mesh->iterate(it)))
      mesh->setModelEntity(ent, interior);
    mesh->end(it);
  }

If this is correct, then I think what it's doing is saying interior_faces are within the volume, so they're on model_dim=3, where as boundary_faces are on the boundary so they're on model_dim=2. What I can't understand is why I've not had to define the edges to be on model_dim=3 in the interior volume, and on model_dim=2 on the boundary, indeed they're marked as model_dim=1 which was required to get it to work. Likewise for the vertices? Could someone help clarify my misunderstanding?

[mortezah]

@a-jp, the classification of mesh entities follows a specific rule which is defined in some of the papers published by SCOREC (for example look at section 2.2 of https://scorec.rpi.edu/reports/view_report.php?id=666 ). But in short, each mesh entity is classified on the lowest order model entity. For example, if an edge is on a model edge it will have model dimension 1, and if it is classified on a model face it will have model dimension 2 and so on.

There is another problem in your code (not related to classification) if you run it in parallel. If I recall correctly, countUpward called on a mesh face that is on the part boundary will give you the answer 1, even for an interior mesh face.

[a-jp]

There is another problem in your code (not related to classification) if you run it in parallel. If I recall correctly, countUpward called on a mesh face that is on the part boundary will give you the answer 1, even for an interior mesh face.

That's frustrating I wasn't aware of that. I had expected that one of the face-cells would have been the ghost cell, and the other the interior cell. I had expected that countUpward always returned 2 for a parallel mesh face on a part boundary. Exactly, how am I meant to index ghost cells then if not from the face on the partition boundaries?

As to your other point, thanks for the help. Yes I'd read that page just don't understand it. It's all very well from a simple 2D drawing defining that. But for a very complex mesh in 3D that already exists, I'm struggling to work out how to define the edges/vertices properly, as that amounts to identifying feature edges in the mesh (to define the edges properly with a model dimension of 1) and corners (to define the vertices with a model dimension of 0). If I read a mesh in from say vtk none of that information is there. I'm trying to apply the classification after the fact without CAD to fall back on. In that instance it's not clear how to define the edges or vertices on the boundary. Do you have sample code for this that I could integrate with the above posted code?

The logical conclusion of your comment is that for edges and vertices fully contained in the interior they should be defined on a model dimension of 3. Is that correct?

[mortezah]

The problem with countUpward is relatively easy to rectify. The only thing you need to do to figure out if a face is interior or not is to change your check to the following:

If countUpward returns 2 OR if the face has a remote copy, then it must be interior. (To get remote copies, check the api getRemotes.)

The struggle that you are facing to classify mesh entities manually is exactly the reason that at SCOREC we recommend to start the mesh generation process with a CAD model and use a tool like Simmetrix to generate a mesh that is already "properly" classified on the model. (Note that this approach will always work no matter the complexity of the model).

Assuming you don't have access to the above workflow you can do either of the following:

1. if you know model entities (like verts, edges, faces, and regions) and also know which mesh entities each one includes, you can do the classification manually, following the same definition in the papers that I sent you earlier. Note that even though things are described in the context of a 2D example in that paper, they work exactly the same for a 3D mesh, and once again the complexity of the domain is not an issue.
2. If your intention is to pass the verify test, you can use the api apf::deriveMdsModel or apf::deriveMdlFromManifold. The first one classifies everything on the boundary to a dummy model face and everything interior to a dummy model region. The second one uses some extra information to come up with a model that is a bit more sophisticated (but note that it only works for very simple domains). Please refer to the api documentation here for more details https://www.scorec.rpi.edu/pumi/doxygen/index.html