Conversation
tupek2
left a comment
tupek2
left a comment
There was a problem hiding this comment.
Please double check some slight additional logic. Make sure the exterior and interior areas are as expected. Is this run in parallel?
| if (!mesh.GetFaceInformation(f).IsInterior()) { | ||
| switch (geom) { | ||
| case mfem::Geometry::SEGMENT: | ||
| edge_id++; |
There was a problem hiding this comment.
Why are we incrementing the counter on exterior edges that we are not included in the interior edges set?
There was a problem hiding this comment.
It's because in the case that you have both an exterior and an interior boundary domain, and you don't increment values for the one you are excluding, you end up with domains of the same type that have repeated ID values in them. If you then try to do a domain operation with these two, like a union or difference, the repeated values cause an error. The incrementing on the edges we are skipping over makes sure that the two boundary domains have unique ID values.
There was a problem hiding this comment.
Referring to my comment above, I think that the Type should be different for interior and exterior faces, which resolves the problem of handling set operations (by forbidding them). If we agree on that, then the old logic should go back in: you won't count the exterior edges/faces in the internal boundary sets, and vice versa.
| { | ||
| assert(mesh.SpaceDimension() == d); | ||
|
|
||
| Domain output{mesh, d - 1, Domain::Type::BoundaryElements}; |
There was a problem hiding this comment.
Should the type be Domain::Type::InteriorFaces? I think that's the intent of that type, to differentiate between boundary and inter-element faces. @lihanyu97, you've worked on this more recently, can you confirm? My thinking is that Domain is meant for integrating over. We sort the domains into types of like elements so that Functional objects can iterate over all the elements and repeat the same calculation efficiently, not have to test each element to choose the right logic and data layout sizing.
There was a problem hiding this comment.
Yes, Domain::Type::BoundaryElements only refer to those faces/edges on the boundary of the global domain. Domain::Type::InteriorFaces refer to those faces/edges in the interior of the domain. So if a face/edge is in the interior of the global domain but on the local processor boundary, it is also considered a InteriorFace
| if (!mesh.GetFaceInformation(f).IsInterior()) { | ||
| switch (geom) { | ||
| case mfem::Geometry::SEGMENT: | ||
| edge_id++; |
There was a problem hiding this comment.
Referring to my comment above, I think that the Type should be different for interior and exterior faces, which resolves the problem of handling set operations (by forbidding them). If we agree on that, then the old logic should go back in: you won't count the exterior edges/faces in the internal boundary sets, and vice versa.
| #include "smith/smith_config.hpp" | ||
| #include "smith/mesh_utils/mesh_utils.hpp" | ||
|
|
||
| #include "smith/smith_config.hpp" |
There was a problem hiding this comment.
| #include "smith/smith_config.hpp" |
| // faces that satisfy the predicate are added to the domain | ||
| for (int f = 0; f < mesh.GetNumFaces(); f++) { | ||
| auto geom = mesh.GetFaceGeometry(f); | ||
|
|
There was a problem hiding this comment.
Calling the predicate to determine whether to add a interior face or not can be done here. Instead of calling GetBdrAttribute, you can directly point to the face/edge and call GetElementAttribute, assuming that in cubit you have a unique attribute assigned to the set of interior faces that you want to integrate on.
if one of your interior face is on the processor boundary, then it's going to get counted twice by two different processors and your QoI is integrated twice on that interior face. That's why I implemented insert_shared_interior_face_list which should always be invoked and present if you have a domain of interior faces. I have configured functional and functional qoi accordingly such that the integration is scaled by 0.5 so that when they get integrated twice, the global sum equal to one. This behavior is not related to whether you are on H1 or L2 spaces.
Finally, function makes element restriction to the domain you pass in. If geometry_dim + 1 = space_dim, indicating it's a face integral, ElementRestriction will call GetFaceDofs. The old version of the function is not deleted yet so it might confuse you. We are now calling the newer version here:
I implemented the DG part to grab DoFs on both side of the face. But you can see the search is handled differently. If the space is H1 or HCurl (line 607), it directly calls the native mfem function and you don't get dofs on both side.
| /////////////////////////////////////////////////////////////////////////////////////// | ||
|
|
||
| template <int d> | ||
| static Domain domain_of_boundary_elems(const mesh_t& mesh, |
There was a problem hiding this comment.
I would leave this method unchanged.
5 participants
Adding the ability to create a domain out of an interior boundary, as the previous method to create a domain of boundary elements excluded all interior elements.