Hi, John, Thanks for your helpful suggestions. Having tested several cases, I found that the problem comes from assigning back the active side to the FEMContext c.side. Now things look good.
Many thanks again for your quick response and helps! Best regards, Eli ________________________________ From: John Peterson <jwpeter...@gmail.com> Sent: Friday, August 19, 2022 12:26 AM To: Yu Zhong <elijahzh...@live.com> Cc: libmesh-users@lists.sourceforge.net <libmesh-users@lists.sourceforge.net> Subject: Re: [Libmesh-users] imposing boundary conditions inside the domain On Thu, Aug 18, 2022 at 5:20 AM Yu Zhong <elijahzh...@live.com<mailto:elijahzh...@live.com>> wrote: Dear all, For a 2D problem, I'm trying to impose Dirichlet BC on a physical boundary connecting two subdomains in a full domain (containing more than 2 subdomains). With the penalty method, I tried with c.has_side_boundary_id(ID) to find this physical boundary, where ID is the physical curve ID defined in gmsh. However, it does not work. This method seems only working with outer boundaries instead of those within the full domain. I went through the examples, and found the Miscellaneous Example 9 provides a solution for this issue, where two different meshes are created in two subdomains. Then the boundary connecting these two subdomains can be found as those outer boundaries for a single domain. This method would ask for heavy works on the mesh generation side when there are many such interior B.C. needed within the full domain. I just wonder whether libmesh has any routine that can automatically find the interior edges with specific IDs, such that specific B.C. can be imposed onto these edges. Hi, Libmesh has support for interior sides, so there is most likely some bug, either on our part of yours, that prevents this from working. Off the top of my head, I see that the has_side_boundary_id() function only works for the FEMContext's currently active side, bool FEMContext::has_side_boundary_id(boundary_id_type id) const { return _boundary_info.has_boundary_id(&(this->get_elem()), side, id); } so you definitely need to be sure that you are calling this function in a loop over all elem sides (and not just exterior boundary sides). Another possibility is that the Gmsh reader is not properly storing your internal sides, so one debugging idea would be to call BoundaryInfo::print_info() after reading in your Mesh, and just confirm that all the (elem, side, id) tuples that you expect to be present are actually present. You will most likely need to share either the real code or (preferably) a simplified example which demonstrates the error in order for us to debug the problem further. -- John _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
