Multiphysics: Solving PDEs on subdomains

[rsn766]

HI!

I have been messing around with this page Multiphysics: Solving PDEs on subdomains — FEniCS Workshop

And when I copied and pasted the second part where it graphs the boundaries, I do not get the same graph as yours. wondering if it has to do with my mac book?

from mpi4py import MPI
import dolfinx
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 40, 40)

def Omega1(x, tol=1e-14):
    return x[0] <= 0.7 + tol

def Omega0(x, tol=1e-14):
    return 0.7-tol <= x[0]

tdim = mesh.topology.dim
cell_map = mesh.topology.index_map(tdim)
num_cells_local = cell_map.size_local + cell_map.num_ghosts
marker = np.empty(num_cells_local, dtype=np.int32)
heat_marker = 1
stokes_marker = 3

marker[dolfinx.mesh.locate_entities(mesh, tdim, Omega0)] = heat_marker
marker[dolfinx.mesh.locate_entities(mesh, tdim, Omega1)] = stokes_marker

cell_tags = dolfinx.mesh.meshtags(mesh, mesh.topology.dim,
                                    np.arange(num_cells_local, dtype=np.int32),
                                    marker)

import sys, os

import pyvista

if sys.platform == "linux" and (os.getenv("CI") or pyvista.OFF_SCREEN):
    pyvista.start_xvfb(0.05)


def plot_mesh(mesh: dolfinx.mesh.Mesh, tags: dolfinx.mesh.MeshTags=None):
    plotter = pyvista.Plotter()
    mesh.topology.create_connectivity(tdim-1, tdim)
    if tags is None:
        ugrid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
    else:
        # Exclude indices marked zero
        exclude_entities = tags.find(0)
        marker = np.full_like(tags.values, True, dtype=np.bool_)
        marker[exclude_entities] = False
        ugrid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh, tags.dim, tags.indices[marker]))
        print(tags.indices[marker], tags.values[marker])
        ugrid.cell_data["Marker"] = tags.values[marker]

    plotter.add_mesh(ugrid, show_edges=True, line_width=3)
    plotter.show_axes()
    plotter.view_xy()
    plotter.show()

plot_mesh(mesh, cell_tags)

mesh.topology.create_connectivity(tdim-1, tdim)
facet_map = mesh.topology.index_map(tdim-1)
num_facets_local = facet_map.size_local + facet_map.num_ghosts
facet_values = np.zeros(num_facets_local, dtype=np.int32)
outer_marker = 4
facet_values[dolfinx.mesh.exterior_facet_indices(mesh.topology)] = 4

def inlet(x, tol=1e-14):
    return np.isclose(x[0], 0.0) & ((x[1] >= 0.4-tol) & (x[1]<=0.6+tol))

inlet_marker = 1
facet_values[dolfinx.mesh.locate_entities_boundary(mesh, tdim-1, inlet)] = inlet_marker

def outlets(x, tol=1e-14):
    return ((np.isclose(x[1], 0.0) &  ((0.4 - tol <= x[0]) & (x[0]<=0.6 + tol))) |
            (np.isclose(x[1], 1.0) &  ((0.2 - tol <= x[0]) & (x[0] <= 0.35 + tol))))

outlet_marker = 2
facet_values[dolfinx.mesh.locate_entities_boundary(mesh, tdim-1, outlets)] = outlet_marker

interface_marker = 3
facet_values[dolfinx.mesh.locate_entities(mesh, tdim-1, lambda x: np.isclose(x[0], 0.7))] = interface_marker

facet_tags = dolfinx.mesh.meshtags(mesh, tdim-1, np.arange(num_facets_local, dtype=np.int32), facet_values)
plot_mesh(mesh, facet_tags)

with dolfinx.io.XDMFFile(mesh.comm, "tags.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags, mesh.geometry) 
```

[dokken]

What version of FEniCS are you running (and how did you install it)?
Could you additionally provide the outputs of

python3 -c "import dolfinx; print(dolfinx.__version__, dolfinx.git_commit_hash)"

I think that what you are seeing could be related to: [BUG]: vtk_mesh on subset of cells based of function space · Issue #3607 · FEniCS/dolfinx · GitHub

[rsn766]

@dokken

python3 -c "import dolfinx; print(dolfinx.__version__, dolfinx.git_commit_hash)"

Dolfinx version 0.9.0

install by conda.
using:
conda create -n fenicsx-env
conda activate fenicsx-env
conda install -c conda-forge fenics-dolfinx mpich pyvista

have a mac book pro m3

This is the readout from program:

[dokken]

Very strange. I cannot reproduce this on my system (linux, ubuntu 24.04).
Could you try opening the facet_tags in Paraview (tags.xdmf) and see if it makes sense there?

[rsn766]

Well that’s weird.

[dokken]

Very weird…
I’m not sure how I can be of any further help at the moment, as the tags seems fine in XDMF, and fine on my system.

[rsn766]

You were actually a lot of help. If it worked on your system and not mine, has to be MAC related problem. Think I tracked it down to this:

if sys.platform == "linux" and (os.getenv("CI") or pyvista.OFF_SCREEN):           pyvista.start_xvfb(0.05)

I’ll try a workaround. I was mainly wondering if my program wasn’t recognizing boundaries, but you allayed my worries!

[rsn766]

Figured out the Paraview problem. I needed to do some post processing to get rid of mesh that also graphed with the boundaries. No big deal.

Still never figured out why my machine wouldn’t graph pyvista. It has worked on other Macs.

[dokken]

For later versions of pyvista (>0.45) one no longer needs the start_xvfb on any system. This has been taken into account in the demo that you are running now (see the webpage for the new code).