On 16 Jul 2021, at 17:50, Venketeswaran, Abhishek (FN) <Abhishek.Venketeswaran@netl.doe.gov> wrote:
Dear Koki and Lawrence Thanks for the detailed response. I would really appreciate it if you could take a quick look at this and tell me if I am subclassing DirichletBC properly. Here, I have a cylinder with an outer radius Ro and axial length Lz. The CollarDirichletBC class finds the nodes of the edge corresponding to radius Ro and z=0.
class CollarDirichletBC(DirichletBC): @utils.cached_property def nodes(self): bcnodes=[] Ro=8.0 atol=1e-3 coords=self.function_space().mesh().coordinates.dat.data_ro rvec=np.linalg.norm(coords[:,0:2],axis=1) zvec=coords[:,2] for ind in range(coords.shape[0]): r=rvec[ind] z=zvec[ind] if np.isclose(z,0.0,atol=atol) and np.isclose(r,Ro,atol=atol): bcnodes.append(ind) return bcnodes
That looks plausible.
Also, Lawrence, Would you be able to give me an outline of the code that I would need to change? This is my first time working with open source software and would appreciate any help I could get.
In firedrake/cython/dmcommon.pyx there is a facet_closure_nodes function. You will want to modify that function to take an additional argument that is the label name (as well as the sub_domain). I recommend defining a "constant" like we have done for FACE_SETS_LABEL and CELL_SETS_LABEL. Then remove the line that does `label = FACE_SETS_LABEL` in the else branch. To get that change visible you'll need to run firedrake-update (to rebuild the cython files). Then you'll want to trace through the call stack from DirichletBC construction through the call to `boundary_nodes` to pass the relevant information through. Lawrence