Mesh refinement: adapt vs refine

1

Hi,

Hello,

I want to use mesh refinement in FEniCS.
While searching on the internet I found two commands: adapt and refine.

To try to understand how these commands work I tried to program a refinement loop that refines the mesh around a given point.
I have two problems:
when using adapt (commented line in the code) I have an error while I added the lines from the post “Equivalent commmand of adapt(space) in fenics 2019”

ETraceback (most recent call last):

File “example_04.py”, line 52, in
mesh = adapt(mesh,markers)
TypeError: adapt(): incompatible function arguments. The following argument types are supported:
1. (mesh_function: dolfin.cpp.mesh.MeshFunctionSizet, adapted_mesh: dolfin.cpp.mesh.Mesh) → dolfin.cpp.mesh.MeshFunctionSizet
2. (bc: dolfin.cpp.fem.DirichletBC, adapted_mesh: dolfin.cpp.mesh.Mesh, S: dolfin.cpp.function.FunctionSpace) → dolfin.cpp.fem.DirichletBC
3. (space: dolfin.cpp.function.FunctionSpace, adapted_mesh: dolfin.cpp.mesh.Mesh) → dolfin.cpp.function.FunctionSpace
4. (form: dolfin.cpp.fem.Form, adapted_mesh: dolfin.cpp.mesh.Mesh, adapt_coefficients: bool = True) → dolfin.cpp.fem.Form

Invoked with: <dolfin.cpp.generation.UnitSquareMesh object at 0x7f808b9cd360>, <dolfin.cpp.mesh.MeshFunctionBool object at 0x7f808a745f70>
rror message:

If I use refine, I don’t understand why it is not the marked cell that is refined (In fact, it works for the first refinement, not for the following ones)

Here is the script I am using.

from fenics import *
from matplotlib import pyplot

cpp_code = “”"
#include<pybind11/pybind11.h>
#include<dolfin/adaptivity/adapt.h>
#include<dolfin/mesh/Mesh.h>
#include<dolfin/mesh/MeshFunction.h>
#include<dolfin/fem/DirichletBC.h>
#include<dolfin/function/FunctionSpace.h>
#include<dolfin/fem/Form.h>

namespace py = pybind11;

PYBIND11_MODULE(SIGNATURE, m) {
   m.def("adapt", (std::shared_ptr<dolfin::MeshFunction<std::size_t>> (*)(const dolfin::MeshFunction<std::size_t>&,
      std::shared_ptr<const dolfin::Mesh>)) &dolfin::adapt,
      py::arg("mesh_function"), py::arg("adapted_mesh"));
   m.def("adapt", (std::shared_ptr<dolfin::DirichletBC> (*)(const dolfin::DirichletBC&,
   std::shared_ptr<const dolfin::Mesh>, const dolfin::FunctionSpace&)) &dolfin::adapt,
    py::arg("bc"), py::arg("adapted_mesh"), py::arg("S"));
   m.def("adapt", (std::shared_ptr<dolfin::FunctionSpace> (*)(const dolfin::FunctionSpace&,
   std::shared_ptr<const dolfin::Mesh>)) &dolfin::adapt,
    py::arg("space"), py::arg("adapted_mesh"));
   m.def("adapt", (std::shared_ptr<dolfin::Form> (*)(const dolfin::Form&,
   std::shared_ptr<const dolfin::Mesh>, bool)) &dolfin::adapt,
    py::arg("form"), py::arg("adapted_mesh"), py::arg("adapt_coefficients")=true);
}
"""

adapt = compile_cpp_code(cpp_code).adapt

n_elem = 1
n_ref = 5
mesh = UnitSquareMesh(n_elem,n_elem,“crossed”)

Praf = Point(0.20981,0.10987652,0.0)
for i_ref in range(n_ref):
print(“-------------------------------------”)
print(“Refinement “+str(i_ref))
print(”-------------------------------------”)
print(" “)
markers = MeshFunction(“bool”, mesh, False)
for cell in cells(mesh):
markers[cell] = cell.contains(Praf)
if markers[cell]:
print(cell.get_vertex_coordinates())
print(”")

plot(mesh)
pyplot.show()
mesh = refine(mesh,markers)

mesh = adapt(mesh,markers)

plot(mesh)
pyplot.show()

Thank you very much.

Xavier

2

I’m not sure if this completely answers your question:

As far as I’m aware, adapt() was introduced a long time ago as part of this work (and perhaps later increments). It was a “cute” interface for a problem which is hard to tame.

refine() is the standard interface for mesh refinement with old dolfin which implements one or two algorithms whose names I forget. I can’t get your code example to run on my machine, but I’m assuming that more than one cell is being refined. This is likely because of the necessity of old dolfin to only allow a conforming mesh structure. So neighbours of cells being refined will need to also be refined to support any floating vertices which are generated.

3

Thank you Nate for your quick answer.

I am aware that it is not possible to divide only one cell in a conforming mesh.
So it’s not a problem for me if refine cuts more than one cell.
But I do not understand why refine cuts cells close to the top left of the unit square while I mark (or try to mark) a cell close to the bottom left.
As it seems I did not copy properly my script in my previous mail (indentations disapear), I copy it a second time in this mail.

Xavier

from fenics import *
from matplotlib import pyplot

cpp_code = """
    #include<pybind11/pybind11.h>
    #include<dolfin/adaptivity/adapt.h>
    #include<dolfin/mesh/Mesh.h>
    #include<dolfin/mesh/MeshFunction.h>
    #include<dolfin/fem/DirichletBC.h>
    #include<dolfin/function/FunctionSpace.h>
    #include<dolfin/fem/Form.h>

    namespace py = pybind11;

    PYBIND11_MODULE(SIGNATURE, m) {
       m.def("adapt", (std::shared_ptr<dolfin::MeshFunction<std::size_t>> (*)(const dolfin::MeshFunction<std::size_t>&,
          std::shared_ptr<const dolfin::Mesh>)) &dolfin::adapt,
          py::arg("mesh_function"), py::arg("adapted_mesh"));
       m.def("adapt", (std::shared_ptr<dolfin::DirichletBC> (*)(const dolfin::DirichletBC&,
       std::shared_ptr<const dolfin::Mesh>, const dolfin::FunctionSpace&)) &dolfin::adapt,
        py::arg("bc"), py::arg("adapted_mesh"), py::arg("S"));
       m.def("adapt", (std::shared_ptr<dolfin::FunctionSpace> (*)(const dolfin::FunctionSpace&,
       std::shared_ptr<const dolfin::Mesh>)) &dolfin::adapt,
        py::arg("space"), py::arg("adapted_mesh"));
       m.def("adapt", (std::shared_ptr<dolfin::Form> (*)(const dolfin::Form&,
       std::shared_ptr<const dolfin::Mesh>, bool)) &dolfin::adapt,
        py::arg("form"), py::arg("adapted_mesh"), py::arg("adapt_coefficients")=true);
    }
    """
adapt = compile_cpp_code(cpp_code).adapt
#
n_elem = 1
n_ref =  5
mesh = UnitSquareMesh(n_elem,n_elem,"crossed")
#
Praf = Point(0.20981,0.10987652,0.0)
for i_ref in range(n_ref):
    print("-------------------------------------")
    print("Refinement "+str(i_ref))
    print("-------------------------------------")
    print(" ")
    markers = MeshFunction("bool", mesh, False)
    for cell in cells(mesh):
        markers[cell] = cell.contains(Praf)
        if markers[cell]:
            print(cell.get_vertex_coordinates()) 
    print("")
#
    plot(mesh)
    pyplot.show()
    mesh = refine(mesh,markers)
#    mesh = adapt(mesh,markers)
plot(mesh)
pyplot.show()
4

See my modifications of your code:

from dolfin import *
import matplotlib.pyplot as plt


n_elem = 1
n_ref =  2
mesh = UnitSquareMesh(n_elem,n_elem,"crossed")
Praf = Point(0.20981,0.10987652,0.0)


fig, axs = plt.subplots(nrows=n_ref+1, ncols=1, constrained_layout=True)

plt.sca(axs[0])
plot(mesh)
plt.scatter(Praf.x(), Praf.y())


for i_ref in range(n_ref):
    print("-------------------------------------")
    print("Refinement "+str(i_ref))
    print("-------------------------------------")
    print(" ")
    markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)
    for cell in cells(mesh):
        markers[cell] = cell.contains(Praf)
        if markers[cell]:
            print(cell.get_vertex_coordinates()) 
    print("")

    mesh = refine(mesh,markers)

    plt.sca(axs[i_ref+1])
    plt.scatter(Praf.x(), Praf.y())
    plot(mesh)


plt.savefig("meshes.png")

Which gives me this:

meshes

The key change here was

markers = MeshFunction("bool", mesh, mesh.topology().dim(), False)

A MeshFunction requires the topology dimension over which it’s defined. By providing False as the third argument, you were essentially creating a vertex function (False is interpreted as 0 which corresponds to vertex topology). This is one of the many reasons I recommend dolfinx over old dolfin since it employs explicit interfaces where possible.

1 Like
5

Thank you Nate.
It is exactly what I’m looking for.

Regards

Xavier

1 Like