Applying DirichletBC to a single node

[kothah]

Hi all,

I have been struggling with enforcing DirichletBC to nodes. I am trying to solve a problem that satisfies boundary conditions on a few nodes. My mesh has marked nodesets but from what I gathered petsc does not read the nodeset information from the exodus meshes. So, I am trying to mark the nodes directly in Firedrake.

Here is a snippet of the code which I have written to extract the node (0.,0.)

mesh = UnitSquareMesh(8,8)
x = SpatialCoordinate(mesh)
sd = SubDomainData(And(real(x[0]) < 1e-4, real(x[1]) < 1e-4))

V = VectorFunctionSpace(mesh, "CG", 1)
bc_bottom_left = DirichletBC(V, Constant(10.), sd)
print(bc_bottom_left.nodes)

This code does not work, the error message says "BC construction got invalid markers."

Can someone please let me know how can this be done?

Best,
Hardik

[pefarrell]

Here's an old (i.e. used to work, interfaces might have changed mildly) snippet of code that colours facets of a mesh:

from firedrake import *
from firedrake.cython import dmcommon

def mark_mesh(mesh):
    dm = mesh.topology_dm
    dim = mesh.geometric_dimension()
    sec = dm.getCoordinateSection()
    coords = dm.getCoordinatesLocal()
    dm.removeLabel(dmcommon.FACE_SETS_LABEL)
    dm.createLabel(dmcommon.FACE_SETS_LABEL)

    nface = dm.getStratumSize("exterior_facets", 1)
    if nface == 0:
        return mesh
    faces = dm.getStratumIS("exterior_facets", 1).indices

    for face in faces:
        vertices = dm.vecGetClosure(sec, coords, face).reshape(dim, dim)  # assumes simplices and geometric dim == topological dim

        if numpy.allclose(vertices[:, 0], 0):
            # all x-coordinates are near 0
            dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 1)
        elif numpy.allclose(vertices[:, 0], 10.0):
            # all x-coordinates are near 10
            dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 2)
        else:
            # otherwise
            dm.setLabelValue(dmcommon.FACE_SETS_LABEL, face, 3)
    return mesh

[kothah]

Hi Patrick,

Thank you for your reply. This code still works for marking faces.

Although, I did not figure out how to mark only a single node.
I would like to set up a simulation, where the load is applied only at a single point.
Do you have any recommendations for how to do so?

Best,
Hardik

[pefarrell]

Hi Hardik,

Yes, here is a subclass of DirichletBC that fixes all dofs at one vertex, in this snippet the origin:

class FixOriginBC(DirichletBC):
    def __init__(self, V, val, subdomain, method="topological"):
        super().__init__(V, val, subdomain, method)
        sec = V.dm.getDefaultSection()
        dm = V.mesh()._plex

        coordsSection = dm.getCoordinateSection()
        coordsDM = dm.getCoordinateDM()
        dim = dm.getCoordinateDim()
        coordsVec = dm.getCoordinatesLocal()

        (vStart, vEnd) = dm.getDepthStratum(0)
        indices = []
        for pt in range(vStart, vEnd):
            x = dm.getVecClosure(coordsSection, coordsVec, pt).reshape(-1, dim).mean(axis=0)
            if x.dot(x) == 0.0: # fix [0, 0] in original mesh coordinates (bottom left corner)
                if dm.getLabelValue("pyop2_ghost", pt) == -1:
                    indices = [pt]
                break

        nodes = []
        for i in indices:
            if sec.getDof(i) > 0:
                nodes.append(sec.getOffset(i))
        self.nodes = numpy.asarray(nodes, dtype=IntType)

        if len(self.nodes) > 0:
            print("Fixing nodes %s" % self.nodes, flush=True)
        else:
            print("Not fixing any nodes", flush=True)

If you want to fix a different vertex change the if x.dot(x) == 0.0.

This uses the coordinates associated to the DM, not the Firedrake coordinates (if you haven't intervened then these will agree).

Use it like bcs = [DirichletBC(V, 0, 1), FixOriginBC(V, 0, "on_boundary")] (the last "on_boundary" is ignored).

[kothah]

Hi Patrick,

This indeed worked! Thank you very much for your help.

Best,
Hardik

[nabw]

Given that this is still useful (at least in my case), here is the same code but working with more recent versions:

class FixOriginBC(DirichletBC):
    def __init__(self, V, val, subdomain, point): # Point must be a numpy array
        super().__init__(V, val, subdomain)
        sec = V.dm.getDefaultSection()
        dm = V.mesh().topological_dm
        coordsSection = dm.getCoordinateSection()
        coordsDM = dm.getCoordinateDM()
        dim = dm.getCoordinateDim()
        coordsVec = dm.getCoordinatesLocal()
        (vStart, vEnd) = dm.getDepthStratum(0)
        indices = []
        for pt in range(vStart, vEnd):
            x = dm.getVecClosure(coordsSection, coordsVec, pt).reshape(-1, dim).mean(axis=0)
            if numpy.linalg.norm(x-point) == 0.0:
                if dm.getLabelValue("pyop2_ghost", pt) == -1:
                    indices = [pt]
                break
        nodes = []
        for i in indices:
            if sec.getDof(i) > 0:
                nodes.append(sec.getOffset(i))
        self.nodes = numpy.asarray(nodes, dtype=int)
        if len(self.nodes) > 0:
            print("Fixing nodes %s" % self.nodes, flush=True)
        else:
            print("Not fixing any nodes", flush=True)

Modifications:

• Removed 'method' due to deprecation
• IntType -> int, the other one didn't exist
• Added 'point' to make it work on an arbitrary point
• _plex -> topological_dm

Thanks again Patrick!

[kothah]

Thanks, Nicolas!

[ABaierReinio]

I've found that the FixOriginBC code above doesn't work when using DG function spaces. In case it's helpful to anyone, here's a different piece of code (written thanks to the help of @pbrubeck) that has the same effect as FixOriginBC, but also works for DG function spaces (at least on triangular meshes).

class FixAtPointBC(DirichletBC):
    def __init__(self, V, g, sub_domain, bc_point):
        super().__init__(V, g, sub_domain)
        self.bc_point = bc_point

    @firedrake_utils.cached_property
    def nodes(self):
        V = self.function_space()

        x_sc = SpatialCoordinate(V.mesh())
        point_dist = interpolate(sqrt(dot(x_sc - self.bc_point, x_sc - self.bc_point)), V)

        with point_dist.dat.vec as v:
            min_index, min_value = v.min()
            assert min_value < 1e-14

            v.zeroEntries()
            v.setValue(min_index, 1.0)

        nodes, = np.nonzero(point_dist.dat.data_ro_with_halos)
        assert len(nodes) <= 1
        nodes = np.asarray([nodes[0]] if (len(nodes) > 0) else [], dtype=int)

        nodes_lgmap = V.dof_dset.lgmap.applyInverse([min_index])
        nodes_lgmap = nodes_lgmap[nodes_lgmap >= 0]
        assert len(nodes_lgmap) <= 1

        assert set(nodes) == set(nodes_lgmap)
        return nodes

[CABA1998]

I found this utility really useful. I'm wondering if it can be implemented with Gmsh and the Physical Point entity.