Problem with Dirichlet BC on edges in a 3D mesh
Hi, I created a a couple of physical groups within GMSH with appropriate labels/tags. I then tried to use the same labels after importing the .msh file in firedrake. I got the following error: Please find below the relevant section in the .py file that generates the gmsh mesh and saves it in .msh file gmsh.option.setNumber("Mesh.MshFileVersion", 2) # Physical curves defined with appropriate tags starting from 1000 gmsh.model.addPhysicalGroup(1, b_o_c_arcs,1001) gmsh.model.addPhysicalGroup(1, b_i_c_arcs,1002) gmsh.model.addPhysicalGroup(1, list(t_o_c_arcs),1003) gmsh.model.addPhysicalGroup(1, list(t_i_c_arcs),1004) # Physical surfaces defined with appropriate tags starting from 2000 gmsh.model.addPhysicalGroup(2, circ_pl_surf,2001) gmsh.model.addPhysicalGroup(2, extsurfs,2002) # Physical volumes defined with appropriate tags starting from 3000 gmsh.model.addPhysicalGroup(3, voltags,3001) gmsh.model.mesh.generate(3) gmsh.write("Pipe3D.msh") Here is the snippet of the .py file executed within firedrake: mesh=Mesh('Pipe3D.msh') bc1=DirichletBC(V,zero,1003) When I execute the solve command I get this error: Error: {1003} are not a valid markers (not in {2001,2002}) I think I have a workaround that is not yet giving me errors wherein I define a class derived from DirichletBC and redefine the @utils.cached_property nodes function. However, It would be great if I could do it using the existing DirichletBC class. I would be very grateful if you could help me out here. Thanks -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran@netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:0311bcdf-830c-456c-a671-ff3e827d18f2]
Hi Abhishek, There seems to be a few missing pieces to naturally apply edge boundary conditions in Firedrake, for instance: * EDGE_SETS_LABEL, which is to map edge markers (in your case 1001, 1002, etc) to associated edge entities in the mesh, is not automatically constructed (I think) (Firedrake does construct FACE_SETS_LABEL, which is the facet counterpart), * `bc.nodes` eventually calls `dmcommon.facet_closure_nodes` (https://github.com/firedrakeproject/firedrake/blob/master/firedrake/cython/d...) with given subdomain markers to collect all facet nodes on which we would like to enforce values on, but we would need to generalise this function for edges. So I'm afraid subclassing DirichletBC probably is the best workaround to apply edge conditions at the moment. Regards, Koki ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Venketeswaran, Abhishek (FN) <Abhishek.Venketeswaran@netl.doe.gov> Sent: Thursday, July 15, 2021 7:12 PM To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Problem with Dirichlet BC on edges in a 3D mesh Hi, I created a a couple of physical groups within GMSH with appropriate labels/tags. I then tried to use the same labels after importing the .msh file in firedrake. I got the following error: Please find below the relevant section in the .py file that generates the gmsh mesh and saves it in .msh file gmsh.option.setNumber("Mesh.MshFileVersion", 2) # Physical curves defined with appropriate tags starting from 1000 gmsh.model.addPhysicalGroup(1, b_o_c_arcs,1001) gmsh.model.addPhysicalGroup(1, b_i_c_arcs,1002) gmsh.model.addPhysicalGroup(1, list(t_o_c_arcs),1003) gmsh.model.addPhysicalGroup(1, list(t_i_c_arcs),1004) # Physical surfaces defined with appropriate tags starting from 2000 gmsh.model.addPhysicalGroup(2, circ_pl_surf,2001) gmsh.model.addPhysicalGroup(2, extsurfs,2002) # Physical volumes defined with appropriate tags starting from 3000 gmsh.model.addPhysicalGroup(3, voltags,3001) gmsh.model.mesh.generate(3) gmsh.write("Pipe3D.msh") Here is the snippet of the .py file executed within firedrake: mesh=Mesh('Pipe3D.msh') bc1=DirichletBC(V,zero,1003) When I execute the solve command I get this error: Error: {1003} are not a valid markers (not in {2001,2002}) I think I have a workaround that is not yet giving me errors wherein I define a class derived from DirichletBC and redefine the @utils.cached_property nodes function. However, It would be great if I could do it using the existing DirichletBC class. I would be very grateful if you could help me out here. Thanks -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran@netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:0311bcdf-830c-456c-a671-ff3e827d18f2]
On 16 Jul 2021, at 17:30, Sagiyama, Koki <k.sagiyama@imperial.ac.uk> wrote:
Hi Abhishek,
There seems to be a few missing pieces to naturally apply edge boundary conditions in Firedrake, for instance: • EDGE_SETS_LABEL, which is to map edge markers (in your case 1001, 1002, etc) to associated edge entities in the mesh, is not automatically constructed (I think) (Firedrake does construct FACE_SETS_LABEL, which is the facet counterpart), • `bc.nodes` eventually calls `dmcommon.facet_closure_nodes` (https://github.com/firedrakeproject/firedrake/blob/master/firedrake/cython/d...) with given subdomain markers to collect all facet nodes on which we would like to enforce values on, but we would need to generalise this function for edges.
FWIW if you don't need variable layer extruded meshes, that code (in dmcommon) is completely generic and will work for marked edges you would need to use the EDGE_SETS_LABEL name and also in dmcommon.complete_facet_labels call DMPlexLabelComplete on the EDGE_SETS_LABEL as well. So I think would actually not be a lot of work to support edge-based dirichlet data. Although perhaps one would want to pick a new name. Lawrence
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 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. Thanks -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran@netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:a14b973b-8fc2-4e72-93d9-c72dff10cffc] ________________________________ From: Sagiyama, Koki <k.sagiyama@imperial.ac.uk> Sent: Friday, July 16, 2021 12:30 PM To: Venketeswaran, Abhishek (FN) <Abhishek.Venketeswaran@netl.doe.gov>; firedrake <firedrake@imperial.ac.uk> Subject: [EXTERNAL] Re: Problem with Dirichlet BC on edges in a 3D mesh Hi Abhishek, There seems to be a few missing pieces to naturally apply edge boundary conditions in Firedrake, for instance: * EDGE_SETS_LABEL, which is to map edge markers (in your case 1001, 1002, etc) to associated edge entities in the mesh, is not automatically constructed (I think) (Firedrake does construct FACE_SETS_LABEL, which is the facet counterpart), * `bc.nodes` eventually calls `dmcommon.facet_closure_nodes` (https://github.com/firedrakeproject/firedrake/blob/master/firedrake/cython/d...) with given subdomain markers to collect all facet nodes on which we would like to enforce values on, but we would need to generalise this function for edges. So I'm afraid subclassing DirichletBC probably is the best workaround to apply edge conditions at the moment. Regards, Koki ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Venketeswaran, Abhishek (FN) <Abhishek.Venketeswaran@netl.doe.gov> Sent: Thursday, July 15, 2021 7:12 PM To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Problem with Dirichlet BC on edges in a 3D mesh Hi, I created a a couple of physical groups within GMSH with appropriate labels/tags. I then tried to use the same labels after importing the .msh file in firedrake. I got the following error: Please find below the relevant section in the .py file that generates the gmsh mesh and saves it in .msh file gmsh.option.setNumber("Mesh.MshFileVersion", 2) # Physical curves defined with appropriate tags starting from 1000 gmsh.model.addPhysicalGroup(1, b_o_c_arcs,1001) gmsh.model.addPhysicalGroup(1, b_i_c_arcs,1002) gmsh.model.addPhysicalGroup(1, list(t_o_c_arcs),1003) gmsh.model.addPhysicalGroup(1, list(t_i_c_arcs),1004) # Physical surfaces defined with appropriate tags starting from 2000 gmsh.model.addPhysicalGroup(2, circ_pl_surf,2001) gmsh.model.addPhysicalGroup(2, extsurfs,2002) # Physical volumes defined with appropriate tags starting from 3000 gmsh.model.addPhysicalGroup(3, voltags,3001) gmsh.model.mesh.generate(3) gmsh.write("Pipe3D.msh") Here is the snippet of the .py file executed within firedrake: mesh=Mesh('Pipe3D.msh') bc1=DirichletBC(V,zero,1003) When I execute the solve command I get this error: Error: {1003} are not a valid markers (not in {2001,2002}) I think I have a workaround that is not yet giving me errors wherein I define a class derived from DirichletBC and redefine the @utils.cached_property nodes function. However, It would be great if I could do it using the existing DirichletBC class. I would be very grateful if you could help me out here. Thanks -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran@netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:0311bcdf-830c-456c-a671-ff3e827d18f2] ******************************************************************** This message does not originate from a known Department of Energy email system. Use caution if this message contains attachments, links or requests for information. ********************************************************************
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
participants (3)
-
Lawrence Mitchell
-
Sagiyama, Koki
-
Venketeswaran, Abhishek (FN)