Only mark the longest edge in a given element for refinement

long_refine.cpp

#include <apf.h>
#include <apfMesh2.h>
#include <apfMDS.h>
#include <apfNumbering.h>
#include <apfShape.h>
#include <gmi_mesh.h>
#include <ma.h>
#include <PCU.h>
#include <pcu_util.h>

static apf::MeshEntity* find_max_edge(apf::Mesh* m, apf::MeshEntity* elem) {
  apf::Downward edges;
  int nedges = m->getDownward(elem, 1, edges);
  apf::MeshEntity* max = 0;
  for (int i = 0; i < nedges-1; ++i) {
    double hi = apf::measure(m, edges[i]);
    double hip = apf::measure(m, edges[i+1]);
    if (hi > hip) max = edges[i];
    else max = edges[i+1];
  }
  return max;
}

struct LongRefiner : public ma::IdentitySizeField {
  LongRefiner(apf::Mesh2* m) : IdentitySizeField(m) {
    mesh = m;
    dim = mesh->getDimension();
    mark_edges();
  }
  bool shouldSplit(apf::MeshEntity* edge) {
    if (stored.count(edge)) return true;
    else return false;
  }
  void indicate_edges() {
    apf::MeshEntity* elem;
    apf::MeshIterator* elems = mesh->begin(dim);
    while ((elem = mesh->iterate(elems))) {
      apf::MeshEntity* edge = find_max_edge(mesh, elem);
      apf::setScalar(tags, edge, 0, 1.0);
    }
    mesh->end(elems);
    apf::accumulate(tags);
  }
  void gather_edges() {
    apf::MeshEntity* edge;
    apf::MeshIterator* edges = mesh->begin(1);
    while ((edge = mesh->iterate(edges)))
      if (apf::getScalar(tags, edge, 0) > 0)
        stored.insert(edge);
    mesh->end(edges);
      }
  void mark_edges() {
    apf::FieldShape* s = apf::getConstant(1);
    tags = apf::createField(mesh, "mark", apf::SCALAR, s);
    apf::zeroField(tags);
    indicate_edges();
    gather_edges();
    apf::destroyField(tags);
  }
  int dim;
  apf::Mesh2* mesh;
  apf::Field* tags;
  std::set<apf::MeshEntity*> stored;
};

int main(int argc, char** argv)
{
  PCU_ALWAYS_ASSERT(argc==3);
  MPI_Init(&argc,&argv);
  PCU_Comm_Init();
  gmi_register_mesh();
  apf::Mesh2* m = apf::loadMdsMesh(argv[1],argv[2]);
  ma::SizeField* size = new LongRefiner(m);
  ma::Input* in = ma::configureIdentity(m, size);
  in->shouldFixShape = false;
  in->shouldSnap = false;
  in->maximumIterations = 1;
  ma::adapt(in);
  delete size;
  m->verify();
  apf::writeVtkFiles("adapted", m);
  m->destroyNative();
  apf::destroyMesh(m);
  PCU_Comm_Free();
  MPI_Finalize();
}