On Sat, 4 May 2019, Manav Bhatia wrote:
> I am working on immersed boundary problems where I cut a FE based on a > level-set function. On each element obtained by the cut operation, I > ask the FE to be initialized, either using a quadrature rule, or by specific > in the QP locations. > > While this forward operation of QP -> compute shape functions, > derivatives and normals (on sides) should work out fine, I'd be scared of the derivatives. If dxi/dx is going off to infinity (in some eigenvalues anyway) then our standard attempt to compute dphi/dx as dphi/dxi*dxi/dx sounds perilous. > the issue is showing > up in compute_face_map() ( > https://github.com/libMesh/libmesh/blob/c4c9fd5450489486fd5f7baf2ec9a8ac0d47fc99/src/fe/fe_boundary.C#L224 > ). > > It appears that compute_face_map() is using the xyz locations of the > quadrature points to figure out non-dimensional location of the > quadrature points here: > https://github.com/libMesh/libmesh/blob/c4c9fd5450489486fd5f7baf2ec9a8ac0d47fc99/src/fe/fe_boundary.C#L754 > > > So, every once in a while the level-set based intersection leads to > sliver cut-cells, which case this inverse-map method to fail. > > I am guessing there is a reason to do things this way, even though the > user may have explicitly provided the quadrature points to begin with. > Do you know if there is a way to bypass the inverse_map for FE reinit? I'm actually not aware of any reason to do things this way. It seems like a bad idea in general, not just in your case! There are cases where we do inverse_map() in AMR problems and for edge and face reinits, because we're too lazy to do direct translation function calculations for every single combination of parent<->child element configuration, and it's easy to just do inverse_map() on one of the xyz coordinates calculated from the other. And that's embarrassing, but at least we have an excuse. But here I don't actually see an excuse. Maybe we wrote the compute_face_map API first, then we added the ability to get proper tangents of 2D elements in 3D space later, and someone wanted to do that without changing the API? I'd say you should just change the APIs for compute_face_map and compute_edge_map. I'll bet nobody is directly using them; they're technically public but they're basically internal. Hmm.. that fix will work for face reinits, and for edge reinits with manually supplied quadrature points, but for edge reinits based on an quadrature rule we're *also* using inverse_map(xyz) for some reason, rather than using the refspace_nodes calculation we do for face reinits. That should probably be fixed one of these days. --- Roy _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
