Find collinear sides of tets in a fracture

[rebpoli]

Hi

I have two sidesets of tets with nodes split, forming a fracture.
I need to pair the sides from one set to the sides of the other set.
Basically, build a map [elem,side] => [elem,side].

What is the best way to make this search?
Does libmesh have some smart tool for that? PointLocator?
Or should I build my own B-tree or similar?

Regards,
Renato

[jwpeterson]

Do the sides share the same nodes, or are there two separate sets of coincident nodes defining the face on either side?

[rebpoli]

The side nodes share the same coordinates, but are different nodes so that the DoFs are split from the mesh geometry.
So I would say we need a geometry search right? The sidesets are assigned different boundary_ids ("frac_p" and "frac_n").

[roystgnr]

Definitely need a geometry search to start in that case. So PointLocator would be the place to begin. And if the performance is enough, that's the place to end too. If the runtime isn't good enough, though, you might have better luck by, at each boundary element after the first, initially checking the point neighbors of the last cross-boundary element you found (O(1), when you're lucky) before falling back to the PointLocator (O(log N)).

[rebpoli]

I believe it is enough... I was afraid of rounding issues. I can get the centroid of a face and look for faces at that point. Can't we get into rounding issues? I am considering if the right search would be "what is the nearest face to this centroid"... Would PointLocator give me that?

[jwpeterson]

So I would say we need a geometry search right?

Yes, I would think so. You could use a PointLocator for this... there's a PointLocator::operator() API that fills a std::set of "candidate" elements for a given search point. It is intended for searching for Points that lie on the face between two neighboring elements. So I would loop over the set of elements in either frac_p or frac_n and do a PL search for the vertex_average() of each face, then look in the set of candidate elements for the one that is not the current Elem.

[rebpoli]

I see. The set of candidates can be the opposite sideset? That would solve the issue I guess:

foreach side in frac_p
found = pointlocator( candidate=frac_n )

[jwpeterson]

The candidate_elements set is an output of the search. You can restrict the returned set by subdomain (in case the elements on either side of the fracture happen to be in different subdomains) but in practice I don't think it's necessary as ideally you will only get two candidates back from each search: the element you are searching from and its (geometric) face neighbor.

[rebpoli]

I see. I will test and check how things go.
Thanks.