Thanks, but this doesn't seem to work if I'm working with first-order (planar triangle) elements, but require higher order basis nodes. Do I have to create a dummy gmsh model with a higher-order mesh to get this information?

I am using the following function call in the python API to get Lagrange basis functions for a triangle:

`_, basis_functions, _ = gmsh.model.mesh.getBasisFunctions(2, integration_coords, f"Lagrange{order}")`

where `integration_coords`

is a set of points on which the basis functions are sampled, and `order`

is the desired basis order.

How can I can extract the coordinates of the *nodes on which the basis functions are defined*? I don't think there's an API function for this, but I'm happy to go digging for lower level functions if pointed in the right direction.

It would be even more helpful if I can also extract a mapping from the nodes on which the functions are defined to the functions themselves in the order in which they appear, so that I know, for example, which of those functions' nodes sit on the triangle's edges and which ones are internal nodes. I know that node ordering is given in the docs for higher order functions, but it's not clear to me if this ordering also corresponds to the Lagrange basis functions.

Thank you.

Thanks, I'll try that approach. However, I also need to make sure that the normal vectors of the mesh triangles always point outwards of an object. With the post-processing approach, for an object with the intersecting surface, the normals will be outward for one of the associated objects and inwards for the other. Therefore, I would still have to compute intersections between all object pairs beforehand, after generating the fragments, to keep track of the cases where normals have to be reoriented. Is there a simpler approach that I'm missing?

I realize that the boolean fragments operation can be used to achieve a conformal 2D mesh on the surfaces of objects which are in contact with each other. I am looking for a way to achieve a conformal 2D mesh on the intersecting surfaces, but where each object has *its own* version of that surface. I.e., there should be distinct mesh triangles associated with each object, and those triangles on the intersecting surfaces should still match up conformally. Nodes need not be distinct, only triangles.

Currently, I can think of a rather complicated way of achieving this by combining BooleanFragments, BooleanIntersection, and Periodic Surfaces, as in the .geo script copied below. I'm actually using the C++ api, but prototyping the ideas in .geo scripts.

My main challenge is automating the process of making a duplicate version of the intersecting surface, making it periodic w.r.t. the original, and then assigning this new copy to one of the objects (lines 40-53 of the below .geo file).

For the case of only two objects, automating this process in the api is relatively simple. However, when there are many objects, any of which may be in contact with each other (but guaranteed not to have a volumetric overlap), the only way I can think of is to compare every possible pair of objects one at a time. Otherwise, if I pass the surfaces of many objects into BooleanIntersections in one go, I don't have a way of knowing which objects are associated with the intersecting surface.

I thought that the 'outDimTagsMap' argument of occ::fragment() in the api might provide some information about the intersecting surfaces and their 'parent' objects, but I don't think that variable contains any surface-related info.

Is there an easier way to achieve this, than comparing pairs of objects one at a time? Alternatively, is there a better approach than using fragments? I thought of just using intersections + periodic surfaces, but I think I would still have the same issue when there are many objects.

Thanks for your time, and for the fantastic software.

.geo file copied below:

```
SetFactory("OpenCASCADE");
// Width, height, length, and corner position of first box
w1 = 20.0e-6;
h1 = 20e-6;
l1 = 400e-6;
x1 = 0.0;
y1 = 0.0;
z1 = 0.0;
// Width, height, length, and corner position of second box
w2 = 40.0e-6;
h2 = 20e-6;
l2 = 400e-6;
x2 = 0.0;
y2 = 200.0e-6;
z2 = 20e-6;
// Characteristic lengths for each box
lc1 = 0.25*w1;
lc2 = 0.5*w2;
// Make boxes
Box(1) = {x1, y1, z1, w1, l1, h1};
Box(2) = {x2, y2, z2, w2, l2, h2};
// Compute fragments
out() = BooleanFragments{ Volume{1,2}; Delete; }{};
// Get surfaces of the resulting volumes
S1() = Boundary{ Volume{1}; };
S2() = Boundary{ Volume{2}; };
Printf("all tags: ", S1());
Printf("all tags: ", S2());
// Obtain all intersecting surfaces, which must be duplicated
out2() = BooleanIntersection{ Surface{S1()}; }{ Surface{S2()}; };
// Duplicate the intersecting surfaces and reverse their normals
out3() = Translate {0, 0, 0} { Duplicata { Surface{out2()}; } };
// Reverse Surface {out3(0)}; // comment out for Gmsh >= 4.7.1
// Ensure that the duplicated surface will have the same mesh as the original one
Periodic Surface {out3(0)} = {out2(0)} Translate {0, 0, 0};
// For the volume for which the intersecting surface has wrongly directed normals, replace the original intersecting surfaces with the duplicated ones, which have oppositely oriented normals
// TODO: how do we automatically figure out which volume had the normals correctly oriented, and which one needs to have its intersecting surface replaced?
// TODO: how do we know which index the original surface is stored at?
// TODO: how do we do this efficiently when there are many objects whose intersections must be found?
S2(4) = out3(0); // for Gmsh >= 4.7.1
// S2(3) = out3(0); // for Gmsh 3.0.6
Printf("all tags: ", S1());
Printf("all tags: ", S2());
Printf("all tags: ", out2());
Characteristic Length { PointsOf {Volume{1};}} = lc1;
Characteristic Length { PointsOf {Volume{2};}} = lc2;
Physical Surface("object1") = { S1() };
Physical Surface("object2") = { S2() };
Mesh.Algorithm = 6;
Mesh 2;
Coherence Mesh;
```