Trying to get the points that lie only on the boundary surface of a solidBody

[ZAARAOUI999]

Hello Mr. Andreas, I were going to calculate the distance (let's say the minimal gap) between two solidbodies (s1 and s2) (I'm trying to get points where a contact/collision can happen). So, I choosed to use the points of two boundaries regions with an activated only_surface flag, but there was sth that stopped me and it's that the points of those regions look like they are not only on the surface!!! Is there sth wrong in the definition of the boundaries regions?
The code:

import felupe as fem
import numpy as np
import matplotlib.pylab as plt
import matplotlib as mpl

mpl.rcParams['pgf.texsystem'] = 'pdflatex'
mpl.rcParams.update({'font.family': 'cambria', 'font.size': 13,
    'axes.labelsize': 15,'axes.titlesize': 15, 'figure.titlesize' : 15})


#Defining meshes
m1 = fem.Rectangle(a=(0,0), b=(10,10), n=10)
m2 = fem.Rectangle(a=(2.5,10), b=(7.5,15), n=5)

#Defining regions
r1 = fem.RegionQuad(m1)
r2 = fem.RegionQuad(m2)

#Defining regions boundaries
rb1 = fem.RegionQuadBoundary(m1, only_surface=True, mask=None)
rb2 = fem.RegionQuadBoundary(m2, only_surface=True, mask=None)

#Defining fields
f1 = fem.FieldsMixed(r1, n=1)
f2 = fem.FieldsMixed(r2, n=1)

#Defining materials
mat1 = fem.NeoHooke(mu=1, bulk=1)
mat2 = fem.LinearElastic(E=206000, nu=0.3)

#Defining solidbodies
s1 = fem.SolidBody(mat1, f1)
s2 = fem.SolidBody(mat2, f2)

#Defining masses
g1 = fem.SolidBodyGravity(s1.field, [9810, 0, 0], 1.522e-9)
g2 = fem.SolidBodyGravity(s2.field, [9810, 0, 0], 7.85e-9)

#Getting points from regions boundaries rb1 and rb2
points1 = rb1.mesh.points
points2 = rb2.mesh.points

#Plotting regions boundaries
plt.title('Regions Boundaries')
plt.xlabel(r'x$\longrightarrow$')
plt.ylabel(r'y$\longrightarrow$')
plt.scatter(points1[...,0], points1[...,1], c='r', label=r'rb1')
plt.scatter(points2[...,0], points2[...,1], c='b', label=r'rb2')
plt.legend()
plt.tight_layout()

The plot:

[ZAARAOUI999]

I have taked a look at the RegionBoundary code and I think that for the 2d mesh, the work should be with edges not cells. I mean that we should find out all edges that have a number of appearance equal to 1, that means that there is no two cells that share one of those edges at the same time. Is my idea clear?

[adtzlr]

Yes, that is how it works. Here is the relevant code:

felupe/felupe/region/_boundary.py

Lines 185 to 198 in 5b77a8f

	if self.only_surface:
	# sort faces, get indices of unique faces and counts
	cells_faces_sorted = np.sort(cells_faces, axis=1)
	cells_faces_unique, index, counts = np.unique(
	cells_faces_sorted, True, False, True, 0
	)
	self._index = index
	self._mask = counts == 1
	else:
	self._index = np.arange(len(cells_faces))
	self._mask = np.ones_like(self._index, dtype=bool)
	self._selection = self._index[self._mask]

[adtzlr]

cells_faces are in fact edges for 2d cells.

[ZAARAOUI999]

Oh, okey, I just mixed things up 😁.

[adtzlr]

Well, I know, the naming is a bit misleading. Sorry.

[adtzlr]

Thanks for investigating into that. If I remember correctly, I fixed something concerning the boundary regions. Please, could you have a look which version you are using? Thanks!

[adtzlr]

In 5.3.0 I made a small change; probably I mixed something up here. Will have a look...

[ZAARAOUI999]

I keep updating FElupe every time a new version is published, just this time I forgot, my vesion is __version__ = "5.3.0".I'm going to update it right now , thanks.

[ZAARAOUI999]

I give it a try and I have sth beautiful, it needs an enhancement in case if it will be helpful, there is a devil for-loop there, this is a part of the code:

def boundary_points(region):
    """Extract points that lie on the boundary of a given region"""
    mesh = region.mesh
    points = mesh.points
    cells = mesh.cells
    ncells = mesh.ncells
    cells = cells.reshape(ncells,-1, 2)
    edges = cells.reshape(-1, 2)
    edges=np.sort(edges,axis=0)
    counter = np.arange(edges.shape[0])
    for i, edge in enumerate(edges):
       mask = np.all(edges==edge, axis=1)
       appear = edges[mask]
       counter[i]=len(appear)

    mask = np.where(counter==1)
    points_mask = edges[mask]
    boundary_points = points[points_mask]
    boundary_points.reshape(-1,2)
    return boundary_points

#getting points from regions r1 and r2
points1 = r1.mesh.points
points2 = r2.mesh.points

#getting boundaries points bp
bp1 = boundary_points(r1)
bp2 = boundary_points(r2)

#plotting regions boundaries
plt.title('Regions (in black) + Regions Boundaries (in red)')
plt.xlabel(r'x$\longrightarrow$')
plt.ylabel(r'y$\longrightarrow$')
plt.scatter(bp1[...,0], bp1[...,1], s=80, marker='s', c='r', label=r'rb1')
plt.scatter(bp2[...,0], bp2[...,1], s=80, marker='o', c='r', label=r'rb2')
plt.scatter(points1[...,0], points1[...,1], s=20, marker='s', c='k', label=r'r1')
plt.scatter(points2[...,0], points2[...,1], s=20, marker='o', c='k', label=r'r2')
plt.legend()
plt.tight_layout()

The plot:

That seems like what I were looking for.

[adtzlr]

The solution is:

A RegionBoundary shares the same points as a default inner Region. This is necessary to get equal-shaped assembled vectors and matrices out of the linear/bilinear integral (weak) forms. If you'd like to extract the points on the boundary, perform np.isin(...) as shown in the code below.

Script

import felupe as fem
import numpy as np
import matplotlib.pylab as plt
import matplotlib as mpl

mpl.rcParams['pgf.texsystem'] = 'pdflatex'
mpl.rcParams.update({'font.family': 'cambria', 'font.size': 13,
    'axes.labelsize': 15,'axes.titlesize': 15, 'figure.titlesize' : 15})


#Defining meshes
m1 = fem.Rectangle(a=(0,0), b=(10,10), n=10)
m2 = fem.Rectangle(a=(2.5,10), b=(7.5,15), n=5)

#Defining regions
r1 = fem.RegionQuad(m1)
r2 = fem.RegionQuad(m2)

#Defining regions boundaries
rb1 = fem.RegionQuadBoundary(m1, only_surface=True)
rb2 = fem.RegionQuadBoundary(m2, only_surface=True)

#Defining fields
f1 = fem.FieldsMixed(r1, n=1)
f2 = fem.FieldsMixed(r2, n=1)

#Defining materials
mat1 = fem.NeoHooke(mu=1, bulk=1)
mat2 = fem.LinearElastic(E=206000, nu=0.3)

#Defining solidbodies
s1 = fem.SolidBody(mat1, f1)
s2 = fem.SolidBody(mat2, f2)

#Defining masses
g1 = fem.SolidBodyGravity(s1.field, [9810, 0, 0], 1.522e-9)
g2 = fem.SolidBodyGravity(s2.field, [9810, 0, 0], 7.85e-9)

#Getting points from regions boundaries rb1 and rb2
points1 = rb1.mesh.points[np.isin(np.arange(rb1.mesh.npoints), rb1.mesh.cells_faces)]
points2 = rb2.mesh.points[np.isin(np.arange(rb2.mesh.npoints), rb2.mesh.cells_faces)]

#Plotting regions boundaries
plt.title('Regions Boundaries')
plt.xlabel(r'x$\longrightarrow$')
plt.ylabel(r'y$\longrightarrow$')
plt.scatter(points1[...,0], points1[...,1], c='r', label=r'rb1')
plt.scatter(points2[...,0], points2[...,1], c='b', label=r'rb2')
plt.legend()
plt.tight_layout()

A RegionBoundary operates on cells, not on faces. This is a decision by design in order to obtain a full-dimensional displacement or deformation gradient on the boundary. However, the quadrature points are of course located on the faces (3d) / edges (2d).

FYI, the boundary mesh looks like this:

fem.mesh.concatenate([rb1.mesh, rb2.mesh]).save("mesh_cells.vtk")

You may construct the 2d-edge-mesh with line-elements by yourself:

fem.mesh.concatenate([
    fem.Mesh(rb1.mesh.points, rb1.mesh.cells_faces, "line"),
    fem.Mesh(rb2.mesh.points, rb2.mesh.cells_faces, "line"),
]).save("mesh_faces.vtk")

As an idea for the future: The last one could be worth to be implemented as a method of RegionBoundary, something like region_boundary.mesh_boundary().

Hope that helps!

[ZAARAOUI999]

That's really what I were looking for a short and a clear solution. Thanks a lot!!!

[adtzlr]

Issue #333 was created based on this discussion. RegionBoundary.mesh_faces() is now implemented, see #334.