Assigning Dirichlet BCs on Specific Nodes

1

Hello,

I apologize in advance that I won’t be able to include an example reproducible code at the moment, but I will try to explain the problem as much as I could and hopefully I could get some clues on how to debug my problem.

I am currently running a 2D simulation of a large mesh (generated by Gmsh) of ~30000 nodes, of which specific nodes (~4000 nodes) on the boundaries and inside the domain are to be assigned a zero Dirichlet BC. I obtained a list of these nodes from Gmsh and if I directly plot these nodes on top of the mesh on any visualization software they match fine. However, when the list is passed to fenics to solve the problem, many of the nodes fail to get assigned with the BCs. Following is an example code of how I incorporate the list of nodes into fenics:

BC_nodes = np.loadtxt('list.txt',delimiter=',')
def specified_nodes(x, on_boundary):
    Status = False
    for i in range(BC_nodes.shape[0]):
        Status = near(x[0], BC_nodes[i,0])  and  near(x[1], BC_nodes[i,1] )
        if Status:
            break
    return Status

The .txt file contains coordinates of all the specified nodes, the first column contains the x coordinates and the second column contains the y coordinates. The above function is being passed to Dirichlet BC with this:

bcs  = DirichletBC(V, Constant(0.0), specified_nodes)

More information on this code, I use Newton solver:

parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

solver = NewtonSolver()
solver.parameters["linear_solver"] = "mumps"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
solver.parameters["report"] = False

And assemble the nonlinear problem with the class:

class PDEEquation(NonlinearProblem):
    def __init__(self, L, a, bcs):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bcs = bcs

    def F(self, b, x):
        assemble(self.L, tensor=b)
        self.bcs.apply(b,x)

    def J(self, A, x):
        assemble(self.a, tensor=A)
        self.bcs.apply(A)

Furthermore, the mesh generated by Gmsh is saved in .msh format and is converted to .xml using dolfin-convert.

Any insight or suggestion on where to debug is greatly appreciated, I am also thinking if there’s a more efficient way to tell fenics these specific nodes instead of looping through the whole list every time the solver is called.

Thank you!

2

Consider the following minimal example:

import fenics as fn
import numpy as np
mesh = fn.UnitSquareMesh(10, 10)
V = fn.FunctionSpace(mesh, "CG", 3)

u = fn.Function(V)

nodes = np.array([[0.2, 0.3], [0.5, 0.5], [0.7, 0.8]])


def specified_nodes(x, on_boundary):
    return any([np.allclose(x, row, rtol=1e-14) for row in nodes])


bc = fn.DirichletBC(V, fn.Constant(5), specified_nodes, "pointwise")

bc.apply(u.vector())
fn.File("u.pvd") << u
2 Likes
4

Thank you! Indeed I just found out one major bug that the nodes read from the .txt file has precision issues. I am fixing that and also going to try your method, I will let you know! Thank you!