[C++] Mesh refinement produces no result

Question

Hello,

In the Poisson demo I tried to refine the mesh in a small circle around the center of the mesh, using the example I found here.

First of all, the example left untouched won't compile, I had to change

for (CellIterator cell(mesh); !cell.end(); ++cell)
--> for (CellIterator cell(*mesh); !cell.end(); ++cell)

and

mesh = refine(mesh, cell_markers);
--> refine(*mesh, cell_markers);

to make it compile.

After that, it turns out the mesh is left totally unchanged, although the terminal in which I launched the simulation reads

Number of cells increased from 200 to 312 (56.0% increase).

My modifications probably messed up something important so that there is no effect on the mesh.

What is then the correct way to refine a mesh?

Here is the modified main.cpp of the Poisson demo I used.

#include <dolfin.h>
#include "Poisson.h"
#include <mshr.h>

using namespace dolfin;

// Source term (right-hand side)
class Source : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    double dx = x[0] - 0.5;
    double dy = x[1] - 0.5;
    values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
  }
};

// Normal derivative (Neumann boundary condition)
class dUdN : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = sin(5*x[0]);
  }
};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS;
  }
};

int main()
{
    auto mesh = std::make_shared<UnitSquareMesh>(10, 10);


        CellFunction<bool> cell_markers(mesh);
        cell_markers.set_all(false);
        Point origin(0.5, 0.5, 0.0);
        for (CellIterator cell(*mesh); !cell.end(); ++cell)
        {
            Point p = cell->midpoint();
            if (p.distance(origin) < 0.2)
        cell_markers[*cell] = true;
        }
        refine(*mesh, cell_markers);

  auto V = std::make_shared<Poisson::FunctionSpace>(mesh);

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

  // Define variational forms
  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);

  auto f = std::make_shared<Source>();
  auto g = std::make_shared<dUdN>();
  L.f = f;
  L.g = g;

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

  // Plot solution
  plot(u);
  interactive();

  return 0;
}

Answer 1

When using shared pointers, consider:

auto mesh2 = std::make_shared<Mesh>();
refine(*mesh2, *mesh, *cell_markers, false);
mesh = mesh2;

Comments

Thank you very much!

Your suggestion works, but I had to change

refine(*mesh2, *mesh, *cell_markers, false);
--> refine(*mesh2, *mesh, cell_markers, false);